862
Estimating Quebec provincial forest resources using ICESat/GLAS Ross Nelson, Jonathan Boudreau, Timothy G. Gregoire, Hank Margolis, Erik Næsset, Terje Gobakken, and Go¨ran Sta˚hl
Abstract: Ground plots, airborne profiling and space lidar (light detection and ranging) measurements of canopy height and crown closure, space radar topographic data, a Landsat cover type map, and a vegetation zone map were used in a model-assisted, two-phase sampling design to estimate the aboveground biomass and carbon resources of Quebec. It was determined that a simple random sampling estimator, with covariance terms added, could be used to quantify the variability of regional Geoscience Laser Altimeter System (GLAS) biomass estimates where interorbit distances are, on average, ‡15 km apart. Prediction error increased standard errors, on average, 24.4%, 4.6%, and 2.8% at the cover type, vegetation zone, and provincial levels, respectively. Inclusion of the covariance term in the calculation of grouped cover type variances increased the vegetation zone standard errors up to 3.7 times and the provincial standard errors 15.6 times. In the southern commercial forests of Quebec, GLAS underestimated ground-based biomass values by 7.3% (stratified linear model) and 10.2% (nonstratified linear model). Quebec forests support 2.57 ± 0.33 gigatonnes of carbon (nonstratified linear model). Approximately 25% of that carbon was found to be located in two southern vegetation zones (northern hardwood and mixedwood), another 25% in two northern vegetation zones (taiga and treed tundra), and the remaining 50% in the boreal zone. Re´sume´ : Des parcelles d’inventaire au sol, des mesures de hauteur de la canope´e et de la densite´ du couvert obtenues avec un lidar profileur ae´roporte´ et un lidar spatioporte´, des donne´es topographiques provenant d’un radar spatial, une carte du type de couvert Landsat et une carte des zones de ve´ge´tation sont utilise´es dans une me´thode d’e´chantillonnage en deux phases assiste´e par un mode`le pour estimer la biomasse ae´rienne et les ressources en carbone du Que´bec. Un estimateur par e´chantillonnage ale´atoire simple, avec l’ajout des termes de covariance, peut eˆtre utilise´ pour quantifier la variabilite´ des estimations re´gionales de biomasse GLAS avec une distance moyenne ‡15 km entre les orbites. L’erreur de pre´diction augmente les erreurs-types, en moyenne, de respectivement 24,4 %, 4,6 % et 2,8 % a` l’e´chelle du type de couvert, de la zone de ve´ge´tation et de la province. L’inclusion du terme de covariance dans le calcul de la variance des types de couvert regroupe´s augmente l’erreur-type de la zone de ve´ge´tation jusqu’a` 3,7 fois et celle de la province jusqu’a` 15,6 fois. Dans la portion sud des foreˆts commerciales du Que´bec, GLAS sous-estime les valeurs de biomasse mesure´es au sol par 7,3 % (mode`le line´aire stratifie´) et 10,2 % (mode`le line´aire non stratifie´). Les foreˆts du Que´bec supportent 2,57 ± 0,33 gigatonnes de carbone (mode`le line´aire non stratifie´). Approximativement 25 % de ce carbone est situe´ dans deux zones de ve´ge´tation du sud (foreˆt de´cidue et foreˆt me´lange´e), un autre 25 % se retrouve dans deux zones de ve´ge´tation du nord (taı¨ga et toundra forestie`re) et le 50 % qui reste est dans la zone bore´ale continue. [Traduit par la Re´daction]
Introduction Imprecision and large discrepancies in our knowledge of global carbon (C) pools, especially in the tropics (DeFries et al. 2002) and in inaccessible boreal forests, provide an impetus to develop robust regional, continental, and global
C measurement techniques using Earth remote sensing technology. Currently, the land and oceans are thought to be absorbing approximately 5.0 gigatonnes (Gt) C per year, or 55% of the 9.1 Gt C/year generated by land-use change (1.5 ± 0.5 Gt C/year) and by our consumption of coal, natural gas, and oil (7.6 Gt C/year). The remainder, 4.1 Gt C, is
Received 5 February 2008. Accepted 19 December 2008. Published on the NRC Research Press Web site at cjfr.nrc.ca on 16 April 2009. R. Nelson.1 Biospheric Sciences Branch, Code 614.4, NASA-Goddard Space Flight Center, Greenbelt, MD 20771, USA. J. Boudreau.2 Centre d’e´tude de la Foreˆt (CEF), De´partement des sciences du bois et de la foreˆt, Universite´ Laval, Que´bec, QC G1K 7P4, Canada. T.G. Gregoire. School of Forestry and Environmental Studies, Yale University, 360 Prospect Street, New Haven, CT 06511-2189, USA. H. Margolis. Faculte´ de foresterie et de ge´omatique, Centre d’e´tude de la Foreˆt (CEF), Universite´ Laval, Que´bec, QC G1K 7P4, Canada. E. Næsset and T. Gobakken. Department of Ecology and Natural Resource Management, Norwegian University of Life Sciences, ˚ s, Norway. P.O. Box 5003, NO-1432, A G. Sta˚hl. Department of Forest Resource Management and Geomatics, Swedish University of Agricultural Sciences, 90183 Umea˚, Sweden. 1Corresponding
author (e-mail:
[email protected]). address: Financie`re agricole du Que´bec – De´veloppement International, 1400, boulevard de la Rive-Sud, St-Romuald, QC G6W 8K7, Canada.
2Present
Can. J. For. Res. 39: 862–881 (2009)
doi:10.1139/X09-002
Published by NRC Research Press
Nelson et al.
added to the Earth’s atmosphere annually. Oceanic uptake is estimated to be approximately 2.2 ± 0.4 Gt C/year, with the remainder being absorbed by a ‘‘missing’’ terrestrial sink (2.8 ± 0.7 Gt C/year). (All estimates and standard errors (SEs) are from Canadell et al. (2007).) An appreciable amount of remote sensing work has been done to locate this sink, and much of this research has been conducted using optical satellite data available globally, e.g., Advanced Very High Resolution Radiometer and MODerate resolution Imaging Spectrometer (MODIS) (Myneni et al. 2001; DeFries et al. 2002). Researchers suspect that much of the C uptake occurs in the northern temperate forests (e.g., hardwood encroachment; Pacala et al. 2001) or in the northern Eurasian forests and North American temperate and boreal zones (Myneni et al. 2001). There is a great deal of uncertainty about the location, extent, and dynamics of this global C sink. The uncertainty is driven, in part, by limitations associated with the use of optical sensors to estimate standing carbon stocks and changes to those stocks. Optical sensors, such as the Thematic Mapper, the Advanced Very High Resolution Radiometer, or the MODIS, can precisely define the location and extent of forested areas that have been cleared or burned, i.e., obvious land cover change, but the spectral measurements cannot be used to accurately estimate the amount of biomass or C lost from these converted areas. Houghton et al. (2001) compared seven different estimates of biomass in the Brazilian Amazon: trans-Amazon estimates varied by a factor of *2.4 (from 39 to 93 Gt C) and were inconsistent with respect to locating areas of high and low biomass. The authors suggest that one way to improve estimation of subcontinental C stocks is to employ airborne or satellite observations that measure vertical canopy structure, e.g., lidar (light detection and ranging) technology. The overall purpose of this study was to assess the accuracy and to characterize the precision of satellite lidar-based estimates of regional biomass and C estimates. In addition, we sought to identify a variance estimator that would capture a portion of the total variance inherent in a sampling system that relies on biomass models within a design-based sampling framework driven by satellite lidar observations. We state from the outset that at this point, we cannot identify an optimal variance estimator because we have no information concerning the true population values of the cover type, vegetation zone, and provincial biomass variances. However, we argue for inclusion of certain error sources in a variance estimator based on the characteristics of the data. Further, we can compare estimators in the absence of ‘‘truth’’ to test for stability among estimators, i.e., to look for a tendency to consistently underestimate or overestimate relative to other variance estimators. We utilized the statistical framework of Boudreau et al. (2008) but with a different biomass model (an explicitly linear model) to address the following specific subobjectives: (1) Quantify the effects of including prediction error and covariance in biomass estimates at the cover type, vegetation zone, and provincial levels. In the context of this study, prediction error is that variation associated with the use of a predictive regression to estimate biomass as a function of laser measurements. (2) Compare three alternative estimators for the variance of
863
estimated regional biomass for a Geoscience Laser Altimeter System (GLAS)-based regional and (or) subcontinental, model-assisted, two-phase sampling procedure. (3) Compare nonstratified versus stratified linear model results to assess the impact of using different linear models to estimate regional biomass. (4) Compare GLAS estimates of biomass with that of ground-based estimates to assess the accuracy and precision of the satellite lidar survey.
Background GLAS as a regional sampling tool for forest resources assessment The Ice, Cloud, and land Elevation Satellite (ICESat)/ GLAS lidar provides near-global observations of forest canopy structure, and these structural measurements can be related to forest height (Lefsky et al. 2005; Sun et al. 2008) and biomass (Lefsky et al. 2005). GLAS is a large-footprint space lidar designed primarily for ice elevation measurements, but much effort is being expended to extend its utility to vegetation measurement. The elliptical area illuminated on the ground by a single GLAS pulse is about 110 m 70 m (Abshire et al. 2005), though the size and shape of the footprint has changed over the course of the ICESat mission. In areas with forests on steep slopes, the forest height measurements are convolved with terrain, and the larger the footprint, the greater the topographic smearing of the ground return. Lefsky et al. (2005) developed equations to predict field-measured maximum tree height as a function of GLAS waveform measurements on three sites in (1) Tennessee — eastern US hardwoods, (2) Oregon — temperate needle-leaf forests, and (3) Brazil — tropical broadleaf forests. The equations employed two GLAS variables to predict field-measured height: the total waveform extent (i.e., signal start to signal end, in metres) and the length of the leading edge of the waveform (in metres). A third variable was derived from digital elevation models to characterize the slope of the local terrain. These equations explained between 59% and 68% of the field height variability and root mean squared errors varied from 4.9 to 12.7 m. They used similar models to predict ground-measured biomass near Santarem, Para, Brazil; their equation explained 73% of the biomass variability, with an root mean squared error of 58 t/ha. Lefsky et al. (2007) suggest the use of a trailing edge metric to mitigate pulse spreading in forests growing on steep slopes, precluding the need for coincident digital elevation models. Nelson et al. (2009) used GLAS data in conjunction with a MODIS land cover map to estimate timber volume on an 811 414 km2 area in south-central Siberia. They trained a neural network to predict ground-measured timber volume using GLAS waveform measurements and used that net to estimate volume for each of the MODIS land cover classes. Considering only those pulses falling on slopes £108, GLAS/ MODIS volumes averaged 163.4 ± 11.8 m3/ha based on GLAS pulses on slopes £108. A comparable ground-based estimate is 146 m3/ha (Shepashenko et al. 1998), a difference of 11.9%. Total regional volume estimates more closely agree given that the lower ground-based per-hectare estimates are spread across a larger forest area, as reported Published by NRC Research Press
864
by Grasia (1990). The ground-based total volume estimate for all forested areas is 7.46 109 m3; the comparable GLAS/MODIS estimate is 7.38 109 m3, a difference of less than 1.1%. Boudreau et al. (2008), working with the same data sets utilized in this current study, employed GLAS as a regional sampling tool to estimate total aboveground dry biomass and aboveground standing C throughout the province of Quebec. They describe a model-assisted two-phase sampling procedure (their Appendix A), where they use an airborne lidar profiler to tie existing ground-plot observations of biomass to GLAS height measures. Their average per hectare provincial biomass estimate was within 4% of an independent estimate developed for the North American boreal forest by Botkin and Simpson (1990). Their GLAS-based northern temperate forest biomass estimate was within 17% of a laurentian highlands estimate and within 3% of a St. Lawrence lowland forest estimate reported by Botkin et al. (1993). Comparison with a ground-based estimate for commercial forests from the Ministry of Natural Resources Quebec (Ministe`re des ressources naturelles, de la Faune et des Parcs Que´bec (MRNFPQ)) suggests an 8.3% GLAS underestimate (GLAS: 75.1 ± 2.2 t/ha versus MRNFPQ: 81.9 ± 0.5 t/ha). Estimating biomass variance Boudreau et al. (2008) calculated estimates of cover type variance using a weighted version of a variance estimator (see eq. A-3 in Boudreau et al. (2008) and eq. 1 in Nelson et al. 2008) that is design-unbiased under simple random sampling (SRS). For convenience, we refer to this estimator as the SRS estimator. Treating a systematic sample as a random sample may lead to significant overestimation of the variance when there is a trend in the sampling frame with respect to the variable of interest (Osborne 1942; Nyysso¨nen et al. 1967, 1971). Areal sampling frames, such as landscapes or selected cover types within landscapes, may exhibit trends with respect to tree size, biomass, and volume, especially if one samples across extensive ecotones, e.g., forest to prairie, taiga to tundra, or across elevation, moisture, or temperature gradients. Of concern in our study is a situation where a biomass gradient runs perpendicular to the north–south GLAS orbits. Such a gradient is not obvious but could exist in Quebec, from the western clay belt near James Bay to the larger maritime forests of the Gaspe´ and St. Lawrence to the east. As early as 1926, Lindeberg (1926) suggested the use of a successive differences (SD) estimator to better estimate variance in a systematically acquired line sample. Subsequent ground-based studies, e.g., Wolter (1985, p. 283), have confirmed Lindeberg’s findings. In the context of a lidar sample, a SD estimator (Scheaffer et al. 1990, see eq. 7.27; Nelson et al. 2008, see eq. 2) estimates the variance of a systematically acquired sample by summing squared differences between spatially adjacent flight lines. The Newton’s Method (NM) estimator (Nelson et al. 2008, see eq. 3; Cochran 1977, see eq. 8.45 for unweighted formula) operates similarly but instead sums squared differences among three adjacent flight lines. Nelson et al. (2008), analyzing a systematically acquired airborne profiling lidar data set in Delaware, showed that a SD or NM estimator could be used in place of a SRS esti-
Can. J. For. Res. Vol. 39, 2009
mator to mitigate, but not completely negate, a tendency of the SRS to overestimate variability when flight lines were relatively closely spaced, i.e., 2–6 km apart. When larger interflight line distances up to 28 km apart were considered, the SRS estimator showed no consistent, conservative tendency, and of the three estimators evaluated, the SRS estimator tracked systematic sampling variation most closely. In the current study, we investigated, in the context of a satellite lidar data set where the average distance between adjacent orbits is 15.6 km, the three variance estimators tested by Nelson et al. (2008) to see if the SRS variance estimates are inflated relative to the SD and NM estimators. (Note: A statistical appendix to this report, similar to the Boudreau et al. (2008) appendix, but that includes formulae for the three variance estimators discussed above as well as for the calculation of prediction error, is available upon request from the primary author.)
Materials The study integrated six data sets, the same data sets described by Boudreau et al. (2008). These data sets included the Landsat Enhanced Thematic Mapper Plus (ETM+) land cover map of the province, a digital map of the six vegetation zones, ground plots, GLAS data, SRTM (Shuttle Radar Topographic Mission) topographic data, and airborne profiling lidar data. The Landsat land cover mosaic and the vegetation zone map were used to stratify the province. Data acquired by an airborne profiling lidar were used to tie ground-plot observations to GLAS height and SRTM topographic measurements. The GLAS instrument serves as the provincial sampling tool to attribute each of the land cover classes in each vegetation zone with estimates of biomass per hectare. The SRTM data were used to correct those GLAS measurements for the effects of topography. Each of these data sets is briefly described below. Landsat land cover Ninety-five Landsat scenes were mosaicked and classified using unsupervised clustering techniques to produce a 23class land cover map of the province for the area south of the tree line, i.e., the northernmost black line in Fig. 1. The Landsat imagery (1998–2003) was resampled to a 25 m grid. The 23 land cover classes are listed in the Fig. 1 legend; the areas associated with each of these cover types and vegetation zones are reported in Table A1. Vegetation zone map A vector map divides the province south of the tree line into the five vegetation zones listed in Table A1 and illustrated in Fig. 1. The vegetation zones are defined by climatic factors and elevation, both of which drive tree species distributions (MRNFPQ 2003). One of the five vegetation zones, the boreal vegetation zone, is bisected by the commercial forest line (Fig. 1), which runs roughly east–west just north of the city of Chibougamau. This partition increases the total number of vegetation zones to six. The boreal zone was divided into a commercial (south) and noncommercial (north) zone to facilitate an accuracy assessment of the GLAS biomass estimates. The ground plots used in the accuracy assessment are located only in the commerPublished by NRC Research Press
Nelson et al.
865
Fig. 1. A 23-class Landsat land cover map of Quebec, with the red line and the most northerly black line denoting the commercial forest line and tree line, respectively. The numbers on the map correspond to the seven areas where ground plots were overflown with the airborne profiling lidar: (1) Trois Rivieres, (2) Mont Laurier, (3) Riviere du Loup, (4) Lac St. Jean, (5) Baie Comeau, (6) Chibougamau, and (7) Radisson. Yellow lines mark locations of GLAS pulses overflown with the airborne profiling lidar. Five vegetation zones are identified on this map, from south to north: northern hardwood, mixedwood, boreal, taiga, and treed tundra. The boreal vegetation zone is bisected by the commercial forest line (red) into commercial boreal (south) and noncommercial boreal (north). The southern Arctic vegetation zone is north of the tree line. (Landsat land cover classification was provided by Natural Resources Canada, Canadian Forest Service, Laurentian Forestry Centre.)
cial boreal zone and in the two vegetation zones south of the boreal. Ground plots Three hundred ground plots located in the northern temperate, mixedwood, and southern (commercial) boreal vegetation zones were selected from a list of 12 000+ MRNFPQ temporary sample plots (TSPs) measured between 2000 and 2004. Each circular, fixed-area plot encompasses 400 m2
and has a plot radius of 11.28 m. The 300 TSPs were purposively selected to cover, as best we could, as wide a range of biomass conditions, tree species conditions, and geography as possible in each of the vegetation zones. The 300 were clustered around the seven airports illustrated in Fig. 1 to facilitate flight logistics. Of these 300 preselected sites, 295 were actually overflown. The remaining five ground plots were not flown because of the presence of conditions that made low-altitude (200 m above ground level) flight measPublished by NRC Research Press
866
urements difficult, e.g., proximity to steep terrain, proximity to antennas and (or) towers, or restricted airspace. Six of the seven airports and their associated plots are located in commercial forest. Only the 31 Radisson plots (site 7, Fig. 1), at 53.88N, 77.68W, are located in noncommercial forest. The forests north of Radisson are relatively inaccessible given the absence of roads north of the La Grande Rivie`re hydroelectric project. PALS data Ranging data from the profiling lidar PALS (Portable Airborne Laser System, Nelson et al. 2003) were used to tie ground-plot observations to the GLAS space lidar observations. PALS was flown over the ground plots and also over a subset of GLAS pulses. Using the ground-plot crossings that passed the quality assessment filters, PALS height and crown-closure measurements served as the independent variables in linear models to predict ground-measured total aboveground dry biomass. Once the ground–PALS model(s) was(were) developed, the model(s) could be used to predict biomass on those GLAS pulses intercepted by the PALS profiling lidar. PALS data were nominally collected at 200 m above ground level at 400 Hz. The system employed a single return system (i.e., one ranging measurement per pulse) that sequentially toggles between first and last returns. At a nominal flight speed of 97 knots, or 50 m/s, the system alternately acquired either a first or last return every 12.5 cm along the profiling track. In the postprocessing phase, adjacent first and last returns were combined to produce a pseudo first and last return at 200 Hz, with a postspacing of 25 cm. With a divergence of 1.7 mrad and a 10 cm transmit lens, the system illuminates a spot on the ground 44 cm in diameter at 200 m above ground level, so there is approximately 50% overlap between adjacent first and last returns. ICESat/GLAS data The GLAS L2a data collection over Quebec took place between 27 September and 18 November 2003. The L2a orbits are illustrated in Boudreau et al. (2008), see fig. 3.b. Because of the limited life expectancy of the three lasers aboard the ICESat/GLAS instrument, one of the three lasers is turned on within time windows that for the most part, reflect the interests of the community rather than the interests of those measuring vegetation. The dates were certainly not ideal from a vegetation standpoint given that some portion of the 138 010 pulses measured forest canopies under conditions that would affect the waveform returns and height measurements. Deciduous forests at different stages of senescence and (or) leaf drop were measured along the north– south GLAS orbits, resulting in increased pulse penetration and a downward bias to the canopy height measurements. Also, some unknown portion of the northerly GLAS height measurements may have been contaminated by recent snowfall, decreasing GLAS canopy height measurements. No attempt was made to mitigate these unknown, downward height biases. Their presence may tend to increase noise in the GLAS-based regressions used to predict regional biomass. We do not, however, expect these factors to impact regional means, since any height measurement changes
Can. J. For. Res. Vol. 39, 2009 Fig. 2. SRTM 90 m data used to characterize topography across Quebec.
would be incorporated into the GLAS-based biomass equations. SRTM data SRTM data were collected during a shuttle mission in February 2000 (Fig. 2). The C-band InSAR (Interferometric Synthetic Aperture Radar) digital topographic data over Quebec were used to infer local topography around individual GLAS pulses. The SRTM coverage extended to 608N; this coverage approximately coincides with the Quebec tree line, which follows a convoluted east–west path between *588 and 598N (the northernmost black line in Fig. 1). The SRTM 90 m data were used to calculate the vertical range (metres) and slope (degrees) in a 3 pixel 3 pixel window centered on each GLAS pulse used in the analysis. This coarse topographic correction was made based on results presented by Lefsky et al. (2005), who found that digital elevation model slope information could be used to mitigate the effects of pulse spreading caused by topography.
Methods The six data sets were integrated using a model-assisted, two-phase sampling design. The sampling process requires that quantitative relationships be established between the ground-based biomass estimates and the GLAS waveform measurements. The ground plot – GLAS pulse relationship is established by (1) flying the airborne laser profiler over selected ground plots; (2) developing an equation to predict plot biomass as a function of profiling heights; (3) flying the profiler over selected, individual GLAS pulses; and then (4) developing an equation relating airborne laser estimates of biomass to GLAS measurements. Once this predictive equation was developed, the GLAS/SRTM data sets were used, Published by NRC Research Press
Nelson et al.
in a model-assisted fashion, sensu Sa¨rndal et al. (1992), to estimate biomass throughout Quebec. Results are reported by cover type within vegetation zone. Using an airborne profiler to predict ground-measured dry biomass Airborne profiling passes (n = 413) were made over 295 MRNFPQ and CFS 400 m2 ground plots. Of the 413 overpasses, only those where the aircraft location was within 10 m of the plot center were considered, and each of these was checked further using coincident videography to ensure that significant land cover change (e.g., fire, harvesting, construction) had not taken place between the time of ground measurement (2000–2004) and the time of the PALS mission (August 2005). Those plots that made this cut were further assessed to make sure that the video record agreed with ground mensuration data and the airborne profiling observations. From 200 m above ground level, the video images an area approximately 60 m 80 m on a side. If the forest canopy within the video frames at the time the plane crosses the ground plot was homogeneous in this 60 m 80 m area, and if the laser profile mirrored the forest canopy characteristics in the video, then the point was used in the regression analysis. This additional data quality check was necessary because the PALS pulse locations contained an XY error — up to 5–20 m — as do the ground-plot center-point locations — up to 15 m. Of the 413 potentially useful overpasses, 207 were used to develop the linear models. Using GLAS and SRTM data to predict airborne profiling estimates of biomass PALS data were acquired along four GLAS orbits totaling 4600 km in length during the summer of 2005. Of the 27 000+ pulses that might have been imaged by the PALS profiler, 8841 GLAS pulses were intercepted. The PALS measurements on each of the 8841 GLAS pulses were extracted and used to predict dry biomass on each GLAS shot. The 8841 pulses were then put through five data quality filters (F1–F5) to delete spurious returns: F1 — GLAS pulses over water. F2 — GLAS pulses where the waveform was saturated (i.e., saturation flag enabled). F3 — GLAS pulses where the version 26 standard product height, h14, was >40 m. F4 — GLAS pulses where the waveform extent, wGLAS, was >50 m. F5 — GLAS pulses acquired over mixed cover types according to a 3 pixel 3 pixel window centered on the pulse location in the Landsat land cover classification. These 1325 pulses were used to develop two explicitly linear models: one based on nonstratified ground–PALS biomass estimates, and one based on stratified ground–PALS biomass estimates. Stratification at the ground–PALS level was considered to determine if stratification at this level improved the accuracy and precision of the associated GLASbased biomass equation. Using GLAS to inventory Quebec biomass The GLAS L2a data acquisition consisted of 138 010
867
GLAS pulses scattered across 97 orbits, both ascending and descending. The 138 010 GLAS pulses were pushed through data quality filters F2, F3, and F4 (listed directly above); 104 044 pulses remained after removing saturated pulses and those pulses with unreasonable, unrealistic height measures. Although land cover purity, i.e., data quality filter F5, was a prerequisite for a given pulse used to formulate GLAS predictive regressions, it was not a prerequisite for biomass inventory, as a purity test would have severely reduced the number of GLAS pulses available for the provincial measurement (see Boudreau et al. 2008, table 3). The land cover affiliation of a particular GLAS pulse was determined by a plurality rule. A 64 m diameter circle centered on the pulse centroid was queried in the geographic information system (GIS) to measure the percentages of the Landsat cover types comprising the pulse. The pulse was assigned to the land cover type that accounted for the greatest percentage of the 64 m circle. Similarly, a majority rule was used to assign pulses to a particular vegetation zone if a given pulse fell on a zone boundary. GLAS-pulse-level estimates of biomass were aggregated to the cover type, vegetation zone, and provincial levels. The sampling framework was designed such that the primary sampling unit was the orbit rather than the individual GLAS pulse. Considering orbital estimates of biomass per hectare for a given land cover type precluded the need to account for the effects of spatial autocorrelation between adjacent and near-adjacent GLAS pulses within orbit. In this context, the number of biomass observations for a given cover type within a vegetation zone equalled the number of orbits where one or more GLAS pulses actually intercepted that cover type. The biomass estimate for a given orbit was weighted by the number of GLAS pulses intercepting that cover type in that orbit. This brings up two important points. Firstly, an orbit that missed a given cover type was treated as ‘‘no information’’; it was not handled as a zero estimate of biomass. Secondly, given the first point and the fact that we do not know beforehand where exactly on the Earth’s surface a particular GLAS pulse lands, since repeated orbital ground tracks vary up to ±1 km, it becomes apparent that a priori knowledge concerning sample sizes does not exist. This particular GLAS sample then and, in fact, all GLAS samples, required poststratification because we cannot know ahead of time (1) the number of pulses collected (some may be saturated or contaminated by cloud), and (2) the land cover types that will be measured. We avoided poststratification issues (Scheaffer et al. 1990, p. 129) by conditioning on the particular sample obtained. This conditioning tied our findings most closely to GLAS samples of similar size (Rao 1985). From a practical standpoint, the random nature of a given cover type’s sample size affected our ability to monitor that cover type over time. Potential problems associated with a varying sample size include unstable means and variances, especially with respect to rare cover types. Variance estimation Our regional estimation procedure employed a modelassisted, two-phase sample design. The model-assisted portion of our statistical framework involved the predictions or estimates of ground-measured, total aboveground biomass based on airborne lidar measurements. The two-phase Published by NRC Research Press
868
sample design included (i) a primary phase, the systematically acquired, province-wide GLAS and SRTM satellite measurement, and (ii) a secondary phase, that subset of GLAS/SRTM locations measured by the airborne profiler. The primary phase was tied to the secondary phase via a second equation that estimated airborne lidar-based biomass based on GLAS and SRTM satellite measurements. We melded the model-based and design-based sampling procedures by adding GLAS/SRTM model prediction error to the design-based observations. The addition of model noise may not be an optimal approach given that the design-based inference was predicated based on the selection of a sample rather than an estimate from a fixed population. It is possible that a purely model-based, mixedmodel approach would work better, but that is a topic for future research. Including model-based variability — prediction error The GLAS/SRTM equation generated a mean biomass estimate for a given GLAS/SRTM observation. We introduced noise, or variability, to that prediction at the individual GLAS pulse level by adding prediction error (e.g., Ott 1988, p. 360; Myers 1990, p. 112). This prediction noise, assumed to be homoskedastic and normally distributed, was added to each GLAS pulse. We report biomass means unperturbed by prediction error throughout the report unless stated otherwise. Prediction error was incorporated in estimates of biomass variance. We employed two models, i.e., ground–air and air–satellite, but only incorporated prediction error from the air–satellite estimation phase. We handled only the primary-phase regression error (i.e., air–satellite) because incorporating both errors would require further processing steps such as the development of a Monte Carlo simulator, an undertaking beyond the scope of this study. Design-based variability Given a nominal 91 day orbit, no clouds, and continuous data acquisition on both ascending and descending orbits, GLAS will collect a very regular, near-global, systematic sample four times per year. The nominal distance between ascending or descending GLAS near-polar, 948 inclination orbits in a 91 day data collection will be cosine(latitude) 29.7 km, where 29.7 km is the distance between adjacent ascending or descending GLAS orbital paths at the equator. The average interorbit distance, considering both ascending and descending orbits, will be half this amount. The actual GLAS data collected over a 53 day period can be quite haphazard, and it is an open question as to whether such a collection of GLAS pulses should be treated as a systematic sample or a random sample. To investigate this question, we calculated the variance of our GLAS/SRTM-based estimates using three different variance estimators. Comparisons were made to see if the SRS estimator was inflated relative to the SD or NM estimators. Such a finding would indicate a generalized east–west trend in forest biomass across Quebec, a situation that would argue for use of the SD or NM estimators to mitigate the effects of this spatial autocorrelation among neighboring orbits.
Can. J. For. Res. Vol. 39, 2009
Covariance One might look at the GLAS L2a acquisition as a stratified random sample, where the 104 044 pulses are selected independently among the 23 land cover classes within vegetation zones. If this were the case, then we could estimate the variances of the vegetation-zone-level estimates (i) of biomass per hectare by summing the weighted variances of the component land cover classes (j), i.e., ni X vb a rðb biÞ ¼ ðwij 2 Þ ½vb a rðb b i j Þ . A similar calculation could j¼1
be made to generate provincial-level estimates. The assumption implicit in such a calculation is that GLAS pulses within a given land cover are selected independently of those in other land cover classes; hence, covariances equal zero. We hold that such an assumption is not merited. GLAS shots are not selected independently. There exists an orbital structure, i.e., a spatial dependence, associated with these GLAS observations. Given that big trees tend to grow near big trees and small near small, especially in studies involving ecotones, we argue that the covariance between land cover estimates cannot be ignored. In general, the argument extends to all regional surveys where profiling or scanning lidars are used as sampling tools. We account for and quantitatively describe the effects of adding these within- and across-orbit biomass covariances using equations A.6,7 and A.10,11 in Boudreau et al. (2008). A more general discussion concerning the inclusion of covariances in the calculation of the variance of a sum may be found in Mendenhall et al. (1990, p. 242). Stratification of the ground–PALS equations Two factors need to be considered with respect to stratification. The first factor involves how the PALS estimates of biomass are calculated. The second factor addresses how the nonstratified and stratified GLAS equations are applied to the 104 044 GLAS pulses that comprise the provincial sample. Regarding the first factor, PALS estimates of biomass for each of the 1325 GLAS pulses overflown by the airborne profiler can be calculated based on nonstratified or stratified ground–PALS equations. One nonstratified GLAS equation is developed by relating the 1325 GLAS/SRTM observations to a biomass value generated using a single, nonstratified biomassground–PALS equation. A second GLAS model, which we call a stratified GLAS equation, is developed by relating the 1325 GLAS/SRTM observations to biomass estimates calculated using a set of stratified ground–PALS models. The second factor deals with how the two biomass PALS = f(GLAS/SRTM) equations are applied to the 104 044 GLAS pulses spread across Quebec. Each of these GLAS pulses has assigned to it a cover type identity based on the Earth Observation for Sustainable Development (EOSD) Landsat ETM+ land cover classification. In the nonstratified approach, we used the nonstratified GLAS equation to estimate biomass on every GLAS shot in each of the 23 land cover types identified in the Landsat imagery. In this scenario, we dictated that any land cover type, including rock–rubble, moss–lichen, exposed – barren land, and water, can support measurable forest biomass. In essence, we placed more confidence in the GLAS measurement, not the Landsat identity. Published by NRC Research Press
Nelson et al.
Conversely, in the stratified approach, we placed our confidence in the Landsat identity rather than in the GLAS measurement. So in this approach, if a particular GLAS shot was identified as a nonforest cover type (water, snow–ice, rock– rubble, barren land, moss–lichen (i.e., bryoids), wetland– herb, herb), we assigned a biomass value of 0 t/ha to that pulse regardless of the GLAS/SRTM measurements. Attempts were made to develop actual stratified GLAS equations by cover type within vegetation zone, but these efforts were unsuccessful. In general, R2 values fell from *0.8 to 0.1 as one moved from the taller, closed canopy forests in the south to the stunted, open forests in the taiga and treed tundra. Given the lack of robustness of the northern GLAS equations, we could only test the effects of stratification in the ground–air phase.
869 Fig. 3. The nonstratified equation relating ground-measured total aboveground dry biomass to PALS height, used to predict total aboveground dry biomass for Quebec. (A) Total aboveground dry biomass, b b , as a function of PALS quadratic mean height, h qa . Black line indicates the predicted biomass value, b b ground ¼ 10:154ðh qa Þ 2:535 . (B) Biomass residuals versus predicted dry biomass in four vegetation zones: northern hardwood and mixedwood, south and north boreal, taiga.
Accuracy assessment Between 1998 and 2003, the MRNFPQ measured 17 852 temporary sample plots throughout the commercial forests of Quebec, i.e., south of the commercial forest (Fig. 1, below red line). Of these, 16 814 plots were located in the EOSD ETM+ land cover classification and were identified as belonging to the conifer, hardwood, or mixedwood cover type in each of the three southern vegetation zones: northern hardwood, mixedwood, and boreal south. Each of these 400 m2 circular ground plots were reanalyzed to estimate the total aboveground dry biomass, using the equations of Lambert et al. (2005). To calculate cover type estimates of biomass within vegetation zone and vegetation-zone-level estimates of biomass, the 16 814 TSPs were analyzed as if they are a stratified random sample. We assumed that the plots were randomly distributed within and among the forested cover types, i.e., cover type covariances = 0. As with the GLAS estimates, GIS weights were used to calculate the vegetation zone means and variances. GLAS-based biomass estimates for these same cover types and vegetation zones were compared with the averaged ground-plot values to determine if the GLAS 95% confidence limits incorporate the ground means.
Results Using an airborne profiler to predict ground-measured dry biomass Of the 413 overpasses over 295 plots, 207 crossings passed the data quality checks. These 207 ground plot – PALS observations were used to generate one nonstratified regression model (Fig. 3) and one set of stratified models (Boudreau et al. 2008, table 2). No patterned response was noted in the residuals of the models. Depending on the cover type and vegetation zone, the following airborne profiling laser variables proved most useful for predicting total aboveground dry biomass measured on the ground plots: h70 = 7th decile height, i.e., that height, in metres, below which 70% of the pulses fall; h a = mean canopy height, all pulses, in metres; h qa = quadratic mean height, all pulses, in metres; h qc = quadratic mean height, canopy hits > 3 m only, in metres; and g = crown closure, in percent (0–100), i.e.,
½ðnumber of canopy hits > 3 mÞ=ðtotal number of pulsesÞ 100:0 Two analyses were conducted to assess the quality of the nonstratified ground–air model. A cross-validation exercise indicated that the difference between biomass predicted using an equation developed using all 207 ground– PALS observations were not significantly different from cross-validated predictions, with the mean difference Published by NRC Research Press
870
Can. J. For. Res. Vol. 39, 2009
Table 1. The GLAS/SRTM equations used to estimate total aboveground dry biomass throughout Quebec. Type of GLAS equation Nonstratified Stratified
Equationa b b ¼ 3:846ðwGLAS Þ 6:587ðfGLAS Þ 0:753ðrSRTM Þ 4:515 b b ¼ 3:635ðwGLAS Þ 5:915ðfGLAS Þ 0:733ðrSRTM Þ þ 2:374
R2 0.60 0.58
n 1325 1325
RMSE 31.99 31.73
pcoeff
