Geostatistical estimation of forest biomass in interior Alaska combining Landsat-derived tree cover, sampled airborne lidar and field observations
Chad Babcock, Andrew O. Finley, Hans-Erik Andersen, Robert Pattison, Bruce D. Cook, Douglas C. Morton, Michael Alonzo, Ross Nelson, Timothy Gregoire, Liviu Ene, Terje Gobakken, Erik Næsset
GGeostatistical estimation of forest biomass in interior Alaska combiningLandsat-derived tree cover, sampled airborne lidar and field observations
Chad Babcock a, ∗ , Andrew O. Finley b , Hans-Erik Andersen c , Robert Pattison c , Bruce D. Cook d , Douglas C.Morton d , Michael Alonzo e , Ross Nelson d , Timothy Gregoire f , Liviu Ene g , Terje Gobakken g , Erik Næsset g a School of Environmental and Forest Sciences, University of Washington, Seattle, WA b Forestry Department, Michigan State University, East Lansing, MI c USDA Forest Service, Pacific Northwest Research Station, Seattle, WA d Code 618, Biospheric Sciences Branch, NASA/Goddard Space Flight Center, Greenbelt, MD e Department of Environmental Science, American University, Washington, DC f Yale School of Forestry and Environmental Studies, Yale University, New Haven, CT g Faculty of Environmental Sciences and Natural Resource Management, Norwegian University of Life Sciences, ˚As, Norway
Abstract
Lidar provides critical information on the three-dimensional structure of forests. However, collecting wall-to-wall laser altimetry data at regional and global scales is cost prohibitive. As a result, studies employinglidar for large area estimation typically collect data via strip sampling, leaving large swaths of the forestunmeasured by the instrument. The goal of this research was to develop and examine the performance of acoregionalization modeling approach for combining field measurements, strip samples of airborne lidar andLandsat-based remote sensing products to predict aboveground biomass (AGB) in interior Alaska’s TananaValley. The proposed modeling strategy facilitates mapping of AGB density across the domain. Additionally,the coregionalization framework allows for estimation of total AGB for arbitrary areal units within the studyarea—a key advance to support diverse management objectives in interior Alaska. This research focuses oncharacterization of prediction uncertainty in the form of posterior predictive coverage intervals and standarddeviations. Using the framework detailed here, it is possible to quantify estimation uncertainty for anyspatial extent, ranging from point-level predictions of AGB density to estimates of AGB stocks for thefull domain. The lidar-informed coregionalization models consistently outperformed their counterpart lidar-free models in terms of point-level predictive performance and total (mean) AGB precision. Additionally,including a Landsat-derived forest cover covariate further improved precision in regions with lower lidarsampling intensity. Findings also demonstrate that model-based approaches not explicitly accounting forresidual spatial dependence can grossly underestimate uncertainty, resulting in falsely precise estimates ofAGB. The inferential capabilities of AGB posterior predictive distribution (PPD) products extend beyondsimply mapping AGB density. We show how PPD products can provide insight regarding drivers of AGBheterogeneity in boreal forests, including permafrost and fire, highlighting the range of potential applicationsfor Bayesian geostatistical methods to integrate field, airborne and satellite data.
Keywords: lidar, Landsat, forest biomass, Bayesian hierarchical models, Gaussian process, cokriging,block kriging, carbon monitoring, linear model for coregionalization, small area estimation
1. Introduction
Coupling remote sensing data with field-based forest measurements via regression frameworks offers thepotential to increase the precision of inventory estimates and provides a mechanism for mapping the spatialdistribution of forest biophysical properties. A plethora of studies show strong relationships between lidarmetrics and forest variables (Asner et al., 2009; Babcock et al., 2013; Finley et al., 2014b, 2017; Lim et al.,2003; Næsset, 2004, 2011). These findings have spurred investment in collecting lidar data for large areas ∗ Corresponding author
Preprint submitted to Remote Sensing of Environment April 26, 2017 a r X i v : . [ s t a t . A P ] D ec rom aircraft and satellites alike. Of particular interest is the use of lidar to assist in the estimation offorest inventory parameters in high-latitude terrestrial ecosystems. From a carbon monitoring perspective,forests in boreal systems may contain large stores of aboveground biomass (AGB) and carbon, but theuncertainty associated with current estimates is extremely high (Bradshaw & Warkentin, 2015; Pan et al.,2011). Understanding that the taiga-tundra ecotone is one of the most vulnerable environmental systems toclimate change and that its boreal forests can contribute substantially to the global carbon cycle, methodsare needed to begin monitoring forest carbon stocks and fluxes for these systems (Gauthier et al., 2015;Magnani et al., 2007; Neigh et al., 2013).Current approaches used by the United States Forest Service’s (USFS) Forest Inventory and Analysis(FIA) program to quantify AGB and carbon stocks in temperate regions rely on extensive, spatially-balancedfield plot probability samples to generate forest inventory estimates with acceptable levels of precision(Bechtold & Patterson, 2005; Woodall et al., 2015). In vast remote landscapes, implementing the estimationtechniques used by the FIA in the contiguous United States becomes prohibitively expensive due to thehigh cost of collecting field inventory data in difficult-to-access boreal regions, e.g., interior Alaska (Barrett& Gray, 2011). A potential solution commonly put forward to reduce the expense of monitoring AGB inboreal forest systems is to augment sparse collections of field samples with remote sensing auxiliary data(Wulder et al., 2012). Lidar-derived measures of forest structure tend to be highly correlated with AGB fieldobservations and, thus, are prime candidates to supplement boreal field campaigns. Additionally, passivesensors such as Landsat can be used to derive remote sensing data products correlated with forest AGB(Kumar et al., 2015). Methodologies leveraging relationships between field and lidar can potentially befurther improved by incorporating Landsat-based products (Margolis et al., 2015; Pflugmacher et al., 2014;Powell et al., 2010; Zheng et al., 2004).Here, we address two challenges encountered when attempting to estimate forest AGB for large areasusing lidar coupled with other remote sensing information: 1) incomplete spatial coverage of remote sensingdata; and 2) prediction uncertainty quantification. Incomplete spatial coverage is a common problem forstudies using airborne or spaceborne lidar over sizable study domains (Andersen et al., 2011; Bolton et al.,2013; Nelson, 2010; Nelson et al., 2004). Model-based methodologies used to link field and lidar data toestimate and map AGB typically require laser altimetry information for the entire spatial domain of interest(Babcock et al., 2015, 2016; McRoberts et al., 2013). The expansive nature of boreal systems, make wall-to-wall collections of airborne lidar data unrealistic. Further, future spaceborne lidar systems are not designedto procure complete coverage information. Rather, these campaigns will collect data for relatively narrowbands along the orbital tracts of the sensors’ host satellite (GEDI, 2014; ICESat-2, 2015). In order to gleanany additional information provided by sampled remote sensing data in a statistically rigorous manner,estimation frameworks that can accommodate incomplete coverage auxiliary information are necessary.The second issue examined here is the problem of obtaining useful estimates of uncertainty about forestAGB stocks using model-based statistical procedures—necessary for decision making with imperfect predic-tions of forest AGB. In design-based estimation frameworks, error is assumed to arise from the samplingdesign, which can be appropriately characterized when plots are selected probabilistically (Cochran, 1977;Thompson, 2002). In model-based inference, error is attributed to the underlying process by which theresponse, e.g., AGB, is generated (Gregoire, 1998; Ver Hoef, 2002). Studies attempting to estimate meansand totals for areal units using ancillary data within a model-based paradigm need to specify frameworksthat reliably accommodate the structure of the data to be modeled. It can be the case that modelers whoattempt to use model-based forest inventory estimation approaches posit potentially unrealistic assumptionsabout the distributional characteristics of model errors, such as independent and identically distributed ( iid )errors. In a spatial context, it is likely the field observations of AGB will be spatially autocorrelated. If theauxiliary information used in the model fails to fully account for the spatial dependence among field ob-servations, model-based approximations of AGB uncertainty can be grossly underestimated (Cressie, 1993;Griffith, 2005).Coregionalization models constructed within a Bayesian hierarchical framework offer a solution to bothabove-mentioned challenges (Gelfand et al., 2004). This class of multivariate spatial regression modelsis designed to predict multiple response variables simultaneously while leveraging spatial cross-correlationstructures between error components of the responses. Further, the model can accommodate spatial mis-2lignment, i.e., missing response variable measurements at some locations. If the lidar data is treated asan explanatory variable (used on the right-hand side of the model as in most lidar studies) predictions areonly possible where lidar data is available. Within a coregionalization model, the lidar-derived metrics canbe treated as additional response variables (moved to the left-hand side) and jointly predicted with the re-sponse of interest, e.g., AGB, across the entire landscape while explicitly modeling the spatially co-varyingrelationship among the predictions within and across locations (Finley et al., 2014a). A coregionalizationframework also allows for the inclusion of wall-to-wall covariates derived from satellite data to assist in thejoint prediction of forest AGB and lidar information.When multivariate coregionalization models are estimated using a Bayesian hierarchical approach, uncer-tainty occurring at all levels of the model can be propagated through to prediction and subsequent estimationof means and totals for areal units (Berliner, 1996; Cressie & Wikle, 2011; Gelfand & Smith, 1990; Hobbs &Hooten, 2015). Forms of multivariate spatial prediction models have been in existence since the 1960s, e.g.,cokriging (Matheron, 1963). These non-hierarchical implementations, however, struggle to effectively dealwith uncertainty associated with spatial covariance parameters, e.g., spatial variances and decays (Diggle& Ribeiro, 2007, Section 7.1.1). Due to increased computational efficiencies gained by ignoring uncertaintyin spatial variability, ‘plug-in’ spatial covariance parameters are used in cokriging interpolation routinesavailable in popular GIS software packages. This limits their use for fully model-based predictive inference(Schelin & Sj¨ostedt-De Luna, 2010).The development of inferential approaches for complex spatial prediction within a statistical frameworkis an active area of research. In a hierarchical modeling context, coregionalization frameworks can be con-structed using random effects that arise from spatially correlated Gaussian processes and partition variabilityinto spatial and non-spatial components (Banerjee et al., 2014; Cressie et al., 2009). When formulated assuch, estimation approaches including Restricted Maximum Likelihood (REML) or Markov chain MonteCarlo (MCMC) become possible in frequentist and Bayesian paradigms of statistical model-based inference,respectively (Ver Hoef et al., 2004). There are advantages to choosing a Bayesian hierarchical approach toinference over counterpart frequentist methods. Access to the full posterior predictive distribution (PPD), aby-product of Bayesian inference, allows for easy posterior summarization of means or totals with associateduncertainty for the full spatial domain in addition to any sub-domains that may be of interest–even underback-transformation (Stow et al., 2006). Access to PPDs facilitate subsequent, i.e., post model-fitting, anal-ysis to inform ecological or management objectives while accounting for prediction uncertainty. However,these increases in flexibility come with substantial increases in computational demand.The aim of this study is to develop and examine the performance of a statistical modeling framework thatcan 1) incorporate partial coverage lidar data and wall-to-wall Landsat products to improve AGB densityprediction; and 2) accommodate spatially structured variability unaccounted for by covariates, therebyallowing for more reliable model-based characterizations of uncertainty (e.g., uncertainty intervals withintended coverage) and improved prediction accuracy. We look to the Tanana Inventory Unit (TIU) ininterior Alaska to explore the potential for the proposed coregionalization model to estimate forest AGBstocking by coupling spatially sparse field inventory, partial coverage lidar and Landsat-derived tree coverdata products in boreal landscapes. The USFS and National Aeronautics and Space Administration (NASA)collaborated on the collection of field inventory data using an augmented FIA sampling design and flight-linesamples of lidar data in 2014. Within this study region, four areal domains containing systematic samplesof field and lidar data serve as study sites in this analysis. We use model comparison to identify strengthsand limitations of candidate modeling approaches and information sources. These comparisons highlighttradeoffs associated with different estimation methods and benefits of coregionalization modeling for large-area estimation of AGB using sampled remote sensing data. We also demonstrate the potential of Bayesianspatial hierarchical models for generating small-area forest parameter estimates and analysis of variabilityin ecosystem structure via PPD summarization. At one site, generating AGB PPDs at a watershed scaleoffered insight regarding the drivers of biomass variability in interior Alaska, including permafrost and fire.We also compare uncertainty intervals generated using the Bayesian coregionalization model with classicalcokriging to highlight advantages to fitting spatial models using Bayesian inference.3 . Background There is increased interest in combining active and passive remote sensing data sources, such as lidarand Landsat, for forest variable mapping over large spatial extents. Lidar is attractive for forest height,volume and AGB prediction because it is able to gather information relevant to the vertical structure offorests. Using spatially explicit height information from lidar can lead to more accurate forest variable mapsthan those based solely on passive remote sensing information provided by Landsat or aerial photography.Ediriweera et al. (2014) explored the use of wall-to-wall lidar and Landsat data to map AGB and found thatmodels incorporating both remote sensing sources performed better than using either alone. However, it canbe expensive to collect wall-to-wall lidar useful for forest structure characterization over large areas, e.g.,state- or country-level. Given this limitation, there has been increased research examining the integrationof sampled lidar data, often along flight-lines or orbital tracts, with wall-to-wall readily available passiveimagery to improve inventory mapping efforts. Deo et al. (2017) explored multiple regression and RandomForest approaches fusing strip-sampled lidar and Landsat metrics to predict AGB and showed improvedmodel fits when both remote sensing sources were used as covariates. Their models, however, could notyield a wall-to-wall map due to incomplete spatial coverage of the lidar-derived covariates. Yava¸slı (2016)examined a two-stage processing approach; first predicting field measured AGB at lidar sample locations,then regressing those lidar based AGB predictions on wall-to-wall Landsat imagery to map AGB outside thelidar sampled areas. This type of two-stage approach has the potential to improve point predictions basedon field observations, lidar samples and wall-to-wall imagery. However, without a way to carry uncertaintyacross successive modeling stages, such approaches fail to provide realistic prediction uncertainty estimatesin final map products. Without a well designed statistical framework able to effectively propagate errorstemming from all processing steps through to prediction, multi-step interpolation routines can only rely onholdout validation metrics to assess uncertainty. Estimates of overall average error provided by such metricscan be inadequate when decision-makers attempt to implement policy or management initiatives using spa-tially explicit map products. Most proposed forest measurement, reporting and verification (MRV) systemsaim to provide interpretable measures of spatially explicit uncertainty that account for all error sources usingstatistical uncertainty intervals. Such MRV systems include the United Nations Programme on ReducingEmissions from Deforestation and Forest Degradation (UN-REDD) and NASA’s Carbon Monitoring System(CMS) (CMS, 2010; UN-REDD, 2009).
Forest inventory estimation is typically conducted using design-based statistical principles, relying onprobability samples of field plots to obtain estimates of inventory parameters, such as mean tree density,total growing stock volume or AGB. Most national forest inventory programs, including the USFS FIA, havebeen developed with a design-based approach to uncertainty quantification in mind. The advent of readilyavailable remote sensing information has led to new developments in design-based estimation approachesincorporating models to improve estimation. Concerning the use of incomplete spatial coverage lidar, studiesemploying model-assisted estimators within a design-based inferential paradigm have been proposed thatassume multi-phase data collection schemes (Gregoire et al., 2011; Saarela et al., 2015; Gregoire et al., 2016).In these studies, a model is employed to relate field and remote sensing data and differences are analyzed usingdesign-based variance estimators to approximate error associated with mean or total AGB estimates. Whenattempting to use remote sensing auxiliary data to improve AGB estimates in a design-based paradigm, oneonly requires remote sensing data to cover the plot locations, which can be accomplished by collecting flight-line strips of auxiliary information (Nelson et al., 2012). These model-assisted approaches show promise butthere are some significant limitations to their practical use in remote sensing-aided forest inventory. Design-based model-assisted approaches can be insufficient due to the restrictive sampling designs they require.We identify four shortcomings of design-based inference in reference to using wall-to-wall and/or sampledremote sensing data for forest inventory below. 4. Any model-assisted estimator depends on plots and sampled remote sensing data being establishedas probability samples, which in practice, is not easy to implement (S¨arndal et al., 1992, Chap. 8-9). Although technically a probabilistic sampling approach, the systematic sampling designs whichare typically conducted (or assumed) offer no means to derive a variance estimator. Rather varianceestimators for other sampling designs, such as simple random sampling, are used to approximatesystematic sampling variance with the understanding that uncertainty will likely be overestimated(Cochran, 1946). The simulation study comparing different model-assisted estimation approachespresented in Ene et al. (2012) exemplifies this phenomenon, showing that in the case of two-phasesystematic sampling, uncertainty was over estimated by a factor of four.2. Multi-phase designs require the primary sampling units (field observations) to be a subsample ofthe secondary sampling units (lidar strips/grid cells). This renders any AGB field samples missed bylidar flight-lines or orbital tracks useless for estimation purposes. This restriction becomes particularlyproblematic when combining national-scale forest inventory plots and space-based lidar data to improveAGB estimation. This can result in estimates relying on field observations alone (e.g., FIA data) beingmore precise than those incorporating sampled lidar information simply because the direct estimatoruses all available field samples, i.e., larger sample size.3. Model-assisted approaches can provide unstable small-area or post-stratified estimates of forest vari-ables when excessively low field sample sizes are encountered in sub-domains (Breidenbach & Astrup,2012). Further, design-based confidence intervals require assumptions about asymptotic normality.In small sample size situations, this assumption can be unrealistic. Also, design-based approachesoffer no sound methodology for estimating the variable of interest with uncertainty in sub-domainswhere no field samples were collected. Pfeffermann (2013) provides further general discussion on theshortcomings of design-based estimation strategies concerning small-area estimation.4. Model-assisted approaches leveraging sampled lidar data provide no mechanism for wall-to-wall map-ping of the variable of interest with uncertainty. The ability to generate mean or total estimates ofAGB alongside maps of AGB density, both with associated measures of error, can be advantageousfor forest management and decision-making.
Geostatistical techniques, such as kriging and cokriging (in all their variants), were originally developed inthe early 1960s as statistical geology tools to estimate underground ore reserves based on nearby field samples(Matheron, 1963; Cressie, 1990). Created specifically to obtain optimal interpolation at unknown locations(points or regions) based on ground-truth data, e.g., best linear unbiased predictions (BLUP), geostatisticalapproaches are uniquely relevant to predictive mapping problems throughout environmental remote sensing,including forest inventory applications. Hudak et al. (2002) analyzed five lidar and Landsat-based modelingapproaches to predict forest canopy height and showed that regressions incorporating spatial informationusing geostatistical techniques were better predictors than their non-spatial counterparts. Tsui et al. (2013)examined several spatial modeling approaches, including regression-kriging and cokriging, to interpolatelidar-based predictions of AGB between simulated flight-line strips leveraging wall-to-wall space-based radar.They found that geostatistical modeling incorporating wall-to-wall covariates, e.g., regression-kriging, washelpful for extrapolating between lidar strip samples.The geostatistical studies mentioned above, along with many others found in the remote sensing and forestinventory literature, focus only on optimal interpolation for accurate mapping, disregarding the potentialfor spatial modeling to characterize uncertainty associated with predictions (Meng et al., 2009; Mutanga &Rugege, 2006). For good reason, kriging and cokriging standard errors are often considered underestimatesof uncertainty and therefore ignored (Diggle & Ribeiro, 2007, Section 7). This is due to the use of ‘plug-in’semi-variogram (or autocovariance) model parameters within the kriging/cokriging prediction routine. Byinterpreting classical kriging standard errors and subsequent confidence intervals as statistical measures ofuncertainty, the researcher assumes the estimated spatial covariance parameters, e.g., nugget, partial sill andrange, are known without error. In typical applications this assumption is far from true. Spatial covarianceparameters based on semi- and cross-variogram models depend on many factors, including how observation5airs are binned and variogram model selected to fit to the empirical variogram points, e.g., exponential,spherical or Mat´ern. In addition, the method used to estimate semi- and cross-variogram model parameters,whether by maximum likelihood, eye-fitting or others, can lead to very different parameter estimates andall ignore the inherent variability around the resulting point estimates. Hierarchical Bayesian approachesto spatial model fitting allow researchers to relax these restrictive assumptions, leading to uncertaintycharacterizations that can be more safely interpreted as statistical measures of error. By allowing modelersto specify appropriately vague priors for spatial covariance parameters within a hierarchical framework,Bayesian estimation naturally propagates uncertainty into the posterior (predictive) distribution. Resultinguncertainty characterizations based on the posterior distribution of predictions therefore account for errorabout the spatial covariance parameters, especially when vague prior distributions are used during fitting(Gelman et al., 2013). Bayesian hierarchical modeling approaches are also easily extendable to multivariateresponse settings with incomplete observations, e.g., field AGB and airborne lidar samples (Gelfand et al.,2004).
3. Methods
The four study sites explored in this analysis all fall within the TIU in interior Alaska (Fig. 1). TheTanana Valley State Forest (TVSF) is a 730 000 hectare (ha) tract of predominately boreal forest stretchingalong the Tanana river basin. Nearly 90 percent of the TVSF is considered forested and close to 50 percentof all productive forestland in the TIU is contained within the boundary of the TVSF (Alaska Departmentof Natural Resources, 2016). Tree species typical of taiga forest can be found throughout, including whitespruce (
Picea glauca ), black spruce (
Picea mariana ), tamarack (
Larix laricina ), quaking aspen (
Populustremuloides ) and balsam poplar (
Populus balsamifera ).Tetlin National Wildlife Refuge (TNWR) is nearly 300 000 ha in size with lowland areas characterized byextensive wetlands and poorly drained soils. Wet upland sites are home to black spruce forests whereas drierlandscapes favor white spruce. Deciduous species including quaking aspen, paper birch (
Betula papyrifera )and balsam poplar persist on well-drained south-facing slopes. Shrub vegetation consisting of willow (
Salixspp. ), alder (
Alnus spp. ) and dwarf birch (
Betula spp. ) can be found in lowland areas around water bodies(U.S. Fish and Wildlife Service, 2016).Bonanza Creek Experimental Forest (BCEF) is a Long-Term Ecological Research (LTER) site withinthe TVSF, consisting of vegetation and landforms typical of interior Alaska. The BCEF domain delineatedfor this study is 21 000 ha and includes a section of the Tanana River floodplain along the southeasternborder. The BCEF is a mixture of forest and non-forest vegetation compositions featuring white spruce,black spruce, tamarack, quaking aspen and balsam poplar trees mixed with willow and alder shrublandspecies (Bonanza Creek LTER, 2016a).Caribou-Poker Creeks Research Watershed (CPCRW) is an intensively studied basin reserved for hydro-logical and ecological research (Bonanza Creek LTER, 2016b). CPCRW is approximately 10 600 ha in sizeand divided into 11 watershed units. Many research initiatives at CPCRW involve vegetation and hydrologycomparisons among watershed units (Amatya et al., 2016; Rinehart et al., 2015; Tanaka-Oda et al., 2016).Upland areas are dominated by paper birch and aspen on south-facing slopes. North-facing slopes are largelyoccupied by black spruce. Patches of alder exist in the understory. Lowland sites are composed of mossand dwarf shrubs. CPCRW has a rich source of associated mapped data products including permafrostpolygon layers and fire maps (Chapin & Hollingsworth, 2006; Rieger et al., 1972; Verbyla, 2011). In 2004,the Boundary Fire in interior Alaska burned a significant portion of CPCRW, predominantly along thesoutheastern border of the research area (Hollingsworth et al., 2013, Fig. 1).At all sites, AGB density was not predicted over
Water and
Barren Lands grid cells classified by theNational Land Cover Database (Homer et al., 2015). These areas were not considered to be part of thestudy domains. 6 .2. Field measurements
At TVSF and TNWR, a spatially-balanced systematic design was implemented using a tessellation ofhexagons covering the study areas. The hexagons were approximately 12 141 ha (30 000 acres) in size. Thesepolygons are five times larger than the hexagons used to establish plots in the continental United States. Inthe summer of 2014, FIA field crews established standard FIA plots at the center of each hexagon completewith four subplots; one at the midpoint and three additional subplots extended radially approximately 36meters (m) at 0 ◦ , 120 ◦ and 240 ◦ . Each subplot has an approximate 7.3 m radius and an area of 168.11 m .Field measurements were taken using an augmented FIA inventory design (Bechtold & Patterson, 2005). Anotable change to the typical protocol included the addition of a second micro-plot to be inventoried in eachsubplot. Adding a second micro-plot helped to ensure that a sufficient number of small diameter trees weretallied (between 2.5 and 12.7 centimeter (cm) diameter at breast height) (Andersen et al., 2011). Pattisonet al. (In prep) provide detailed field protocols for the TIU inventory pilot project. The same samplingprotocol was used at BCEF and CPCRW, although at a greatly increased sampling intensity. Field plots atBCEF and CPCRW were each inventoried once in the summer months of 2011, 2012 or 2014. We recognizethe mismatch between field plot measurement and lidar acquisition years for some plots at BCEF andCPCRW is not ideal, but given that boreal forests in this region are slow growing, we argue that this levelof temporal misalignment will have a negligible affect on subsequent results. The total number of subplotsmeasured at each site was 263, 123, 292 and 149 at TVSF, TNWR, BCEF and CPCRW, respectively (notethat some subplots at each plot were not measured for various logistical reasons or fell outside the study area Fig. 1:
Map showing locations of the four study sites within the Tanana Inventory Unit (TIU). plot will be used in subsequent sectionsto refer to the individual FIA subplots described above.
Lidar data were collected using a flight-line strip sampling approach with NASA Goddard’s LiDAR,Hyperspectral and Thermal (G-LiHT) airborne imager in the summer of 2014 (Cook et al., 2013). G-LiHTis a portable multi-sensor system that is used to collect fine spatial-scale image data of canopy structure,optical properties, and surface temperatures. G-LiHT’s on-board laser altimeter (VQ-480, Riegl LaserMeasurement Systems, Horn, Austria) provides an effective measurement rate of up to 150 kilohertz alonga 60 ◦ swath perpendicular to the flight direction using a 1 550 nanometer laser. At a nominal flying altitudeof 335 m, laser pulse footprints were ≈
10 cm diameter and sampling density was ≈ − .Pulse densities > − are able to more accurately characterize complex terrain and vertical distribution ofcanopy elements in dense stands, and tree height metrics (e.g., 90 th percentile height used here) are largelyunaffected by pulse densities above 1 pulse m − (Jakubowski et al., 2013; White et al., 2013). G-LiHT’slidar is capable of producing up to eight returns per pulse. Point cloud information was summarized to a13 x 13 m grid cell size (grid cell area equal to 169 m ) to approximate field plot areas. Over each grid cell,percentile heights were calculated at 10 percent intervals ranging from 10 percent to 100 percent. Maximumheight (100 th percentile height) relative densities were also calculated at 10 equal width intervals. Additionalmetrics including point cloud skewness and kurtosis, among others, were calculated for each grid cell as well.Identical lidar metrics were obtained using point clouds extracted over each field plot. Exploratory regressionmodel fits indicated that 90 th percentile height alone accounted for substantial amounts of variability insquare-root transformed AGB at all sites (0 . ≤ R ≤ . th percentile height and AGB, it was decided that only this metric will be considered for subsequent modelingefforts. G-LiHT data for the study area are available online at https://gliht.gsfc.nasa.gov .To reduce computational demand associated with fitting the proposed models, the gridded lidar covariatesets were subsampled. At the four study sites, approximately 0.5 percent of the original lidar grid cells wererandomly selected. The lidar subset grid cells were combined with the 90 th percentile height values calculatedusing lidar point clouds over the field plots to construct the lidar metric sets for subsequent model-fitting.Specifically, 1 216 of 199 045, 813 of 132 635, 2 011 of 491 621 and 7 069 of 1 374 118 lidar observations wereused for model-fitting at BCEF, CPCRW, TNWR and TVSF, respectively. Increases in computational powerin the future will reduce the need for substantial thinning of the lidar metric set, allowing for the modelstested here to be fitted using more observations. Also, see Section 6 for discussion about future modelingresearch avenues that may lead to increased computational efficiencies, thereby alleviating the need tosubstantially thin sampled remote sensing datasets to implement the modeling frameworks proposed here.We examined several Landsat-derived products to evaluate the potential for additional predictive gainsover the use of lidar data alone. Candidate data layers with wall-to-wall coverage over the study area includedpercent tree cover for 2010 (Hansen et al., 2013), individual reflectance bands from Landsat 8 OperationalLand Imager (OLI) composite surface reflectance products for 2014, and vegetation indices (e.g., TasseledCap values and Normalized Difference Vegetation Index, NDVI) derived from the 2014 Landsat composites.The 2010 percent tree cover data product exhibited the strongest relationship with field observations ofAGB. Further, including 2014 Landsat composite bands and indices in addition to 2010 percent tree coverin regressions did not result in appreciable gains in fit performance. For these reasons, only the 2010 percenttree cover metric was considered for this analysis. R values for preliminary models relating square-roottransformed AGB density and percent tree cover ranged between 0 .
25 and 0 .
55 for the study areas.8 ig. 2:
Overview schematic of the full hierarchical model, labeled as
Coregionalization + Tree Cover in the results tables. Solidblack arrows indicate stochastic linkages for model parameterization, and dashed black arrows show deterministic links basedon input data. To fully specify this framework as a Bayesian hierarchical model, all parameters stochastically linked to otherparameters or observations need to be assigned probability distributions. The blue rectangle contains all model elements (pinkrectangles) associated with the aboveground biomass component and the green rectangle contains all model elements associatedwith the lidar component of the full model. The Landsat tree cover data and cross-covariance matrix belong to both modelcomponents, providing the link between the primary and secondary responses. The levels of the hierarchy are delineated withgray bars highlighting the multi-level structure of the model. .4. Model Overview3.4.1. Coregionalization model explanation Fig. 2 is a graphical depiction of the full model tested in this analysis, labeled as
Coregionalization + TreeCover in the results tables. The other models tested can be viewed as special cases of the full model, formedby setting different model terms to zero. The goal of this section is to provide a high-level description ofthe model framework used in this study. Further details concerning each of the individual tested modelsare provided in Section 3.5 and references therein. We point interested readers to Hobbs & Hooten (2015)for an accessible overview of Bayesian modeling principles in general and Banerjee et al. (2014) for specificdiscussion about the models tested in this analysis—including the coregionalization geostatistical framework.Fig. 2 shows distinct, but connected, components for each response variable—AGB (shown in blue) andlidar derived 90 th percentile height (shown in green)—highlighting that this is a joint modeling framework.Level 1 of the hierarchical framework contains three elements for each response. The trend model elementfor both response components is a simple linear regression with an intercept and slope parameter associatedwith the Landsat derived tree cover data x ( s ). We define s to be a vector of geographic coordinates, e.g.,easting and northing, in the spatial domain D . In addition to the trend model, each component containsa spatial random effect ( v y ( s ) and v z ( s )) designed to model any spatially structured variability left afteraccounting for the trend model. The variance parameters ( τ y and τ z ) account for any uncorrelated erroronce the trend models and spatial random effects are considered.Level 2 of the hierarchical model contains the makings of the Level 1 spatial random effects, i.e., the u y ( s ) and u z ( s ) spatial correlation effects and a shared cross-covariance matrix K . The cross-covariancematrix parameters estimate the spatial variability for each spatial random effect ( k and k ) and also theircovariance k . The covariance parameter, k , is the link between the AGB and 90 th percentile heightmodel components and is informed by observations of spatially coinciding responses (i.e., locations whereboth AGB and lidar derived 90 th are observed). For our setting, estimates of k provide two benefits. First,for locations where neither response is observed, the joint prediction of AGB and lidar 90 th percentile heightwill maintain their observed covariance. Second, in spatial misalignment settings the value of the observedresponse will inform the prediction of the missing response, e.g., prediction of AGB will be improved forlocations where 90 th percentile height is observed along the lidar transects.In addition to informing prediction where one or both of the responses are missing (via the K ) theproposed model informs prediction by pooling information from spatially nearby observations through aspatial correlation function which is informed by a set of spatial correlation parameters that reside inLevel 3. As formalized in the subsequent section, these parameters control the geographic range of spatialdependence exhibited by the random effects. Recall, our goal is to produce spatially explicit predictions of AGB density and estimates of total (mean)AGB, both with quantifiable measures of uncertainty that reflect the fact that none of the estimated pa-rameters are known without error. Such error quantification requires a probabilistic modeling frameworkand some parameter distributional assumptions. Fig. 2 identifies the parameters to be estimated and thehierarchal dependence of the model components, but like most graphical depictions, the distributional as-sumptions are omitted (Lunn et al., 2012). We elect to fit the model described in Fig. 2 using Bayesianmodel-based inference to handle the process of propagating uncertainty stemming from all stages of thehierarchy through to prediction and subsequent estimation of total (mean) AGB.Beginning with the Data Level, we assume that square-root transformed AGB, y ( s ), arises from a normaldistribution with a mean equal to the trend model plus a spatial random effect, i.e., β y + β y x ( s )+ v y ( s ), anda variance parameter, τ y . Similarly, we assume lidar derived 90 th percentile height, z ( s ), to be distributednormally with a mean equal to the trend model plus a spatial random effect, i.e., β z + β z x ( s ) + v z ( s ), anda variance parameter, τ z .In Level 1 of the hierarchy, we see the trend model intercept and slope in addition to the spatial randomeffect are unknown for both components. Moreover, the variance parameters, τ y and τ z , also need tobe estimated. To estimate these parameters within a Bayesian statistical paradigm, we need to set prior10istributions for them. Further, these prior distributions should reflect the researchers belief about whatthese parameters are before observing the data used to fit the model. In our case, we have no prior knowledgeabout the value of the intercept and slope parameters so we establish prior distributions to be normal witha mean equal to 0 and a variance of 10 000, i.e., N (0 ,
10 000). The large variance for these priors makesthem sufficiently vague, i.e., this choice of prior should not influence posterior inference beyond the observeddata. We also have no prior understanding about the true values for the τ y and τ z variance parametersaside from knowing they are non-negative. For the τ y and τ z parameters we assume an Inverse-Gamma( IG ) prior with a mean equal to a reasonable starting value and an infinite variance. The infinite variancesetting establishes these as vague priors. The spatial random effects in Level 1, v ( s ) = ( v y ( s ) , v z ( s )) (cid:48) , areset equal to Au ( s ), where u ( s ) = ( u y ( s ) , u z ( s )) (cid:48) and A is the square-root of the cross-covariance matrix K . This is a deterministic link, meaning that we do not need to specify a distribution for v ( s ) directly.Rather we specify distributions for the unknown parts of v ( s ), namely u y ( s ), u z ( s ) and K .In Level 2, both the u y ( s ) and u z ( s ) random effects are assumed to be distributed according to a Gaussianprocess ( GP ) with a mean equal to zero and a variance defined by a spatial correlation function. Thecorrelation functions chosen for this analysis were exponential decay functions which depend on separationdistances between points in the domain and spatial decay parameters. This class of correlation function isvery useful for modeling spatially autocorrelated variability. When points in the domain are close together,the correlation is near one. The further the two points are from each other, the closer the correlation is tozero. The rate at which the functions approach zero depends on the spatial decay parameters (the Level3, φ y and φ z parameters). The K matrix controls the variance of each of the random effects and also thecovariance between them. Again, K is unknown so we set a prior for it to be an Inverse-Wishart ( IW )distribution. The IW distribution is a multivariate extension of an IG distribution that is useful for settingpriors for variance-covariance matrices. The IW distribution has two parameters that define its shape. Thefirst is a degrees of freedom parameter which we set to two because we are modeling a two-variable (bivariate)variance-covariance matrix. The second parameter for the IW prior is a variance-covariance matrix that weset to . I , where I is a 2x2 identity matrix. This is a vague prior specification. Without standardization,the k covariance parameter (an element of the variance-covariance matrix K ) can be difficult to interpret.Because of this, we present a correlation estimate between the spatial random effects in the results tableslabeled cor ( u y ( s ) , u z ( s )) which is a simple transformation of the estimated cross-covariance matrix.Since the Level 3 spatial decay parameters, φ y and φ z , need to be estimated, sensible priors need to beset for them. We set uniform prior distributions for φ y and φ z with a lower bound near zero and large upperbound in effort to make these priors as vague as possible. Spatial decay parameters are difficult to interpretdirectly, so we present effective range ’s in the results tables ( er y and er z ) which are simply transformationsof φ y and φ z to distance units (km).Now that we have defined the likelihood, random effect and prior distributions for all unknown parameterswe can construct a sampling algorithm to draw from the posterior distribution of our model parameters.For this analysis we used the spBayes package written for the R statistical computing environment whichcontains a class of functions designed to fit models as defined in Fig. 2, namely, the spMisalign function. To assess the utility of coupling plot-level AGB measurements with sampled airborne lidar and the 2010Landsat-based tree cover index developed by Hansen et al. (2013), six candidate models were compared.Candidate models were evaluated on fit to observed data, prediction performance and total (mean) AGBestimation precision. The following sections provide the statistical modeling details required to implementthe six candidate models. The mathematical notation used in the following sections closely mirrors that ofBanerjee et al. (2014).
The
Null model is designed to be a baseline regression where no information beyond plot measurementsis considered. The
Null model is y ( s ) = β y + (cid:15) y ( s ) , (1)11here y ( s ) is a square-root transformed field-based measurement of AGB at location s in spatial domain D . β y is an intercept parameter to be estimated. Since the Null model contains no additional regressionparameters, the estimated value for β y should approximate the overall mean of square-root transformedAGB from the field plots, i.e., (cid:80) ni =1 y ( s i ) /n , where n is the total number of field plots measured. The errorterm (cid:15) y ( s ) captures any departure from the overall mean at site s . Imposing a distributional assumption on (cid:15) y iid ∼ N ( , τ y I )—where (cid:15) y = ( (cid:15) y ( s ) , (cid:15) y ( s ) , ..., (cid:15) y ( s n )) (cid:48) , τ y is a variance parameter to be estimated and I is an n x n identity matrix—allows for model-based statistical inference concerning model parameters andpredictions. Any statistical inference, e.g., credible intervals, posterior standard deviations, etc., evaluatedfor model parameters or predictions can be considered reliable if the distributional assumptions about (cid:15) y are not seriously violated. The
Tree Cover model is y ( s ) = β y + β y x ( s ) + (cid:15) y ( s ) , (2)where y ( s ) and (cid:15) y ( s ) are defined as in model (1). The intercept, β y , and regression slope parameter, β y ,together describe the linear relationship between the Landsat-based tree cover product, x ( s ), and y ( s ). Incorporating a spatial random effect into model (1) establishes the
Spatial model and is y ( s ) = β y + w y ( s ) + (cid:15) y ( s ) , (3)where y ( s ), β y and (cid:15) y ( s ) are defined as in model (1). Here, w y ( s ) is modeled as a Gaussian process ( GP )with a zero mean and a spatial covariance function that captures the covariance between any pair of locations s and s ∗ , i.e., w y ( s ) ∼ GP (0 , C ( s , s ∗ ; θ y )). We specify C ( s , s ∗ ; θ y ) = σ y ρ ( s , s ∗ ; φ y ) where ρ ( s , s ∗ ; φ y ) isa valid spatial correlation function and θ y = { σ y , φ y } , where φ y is a correlation decay parameter and σ y = V ar ( w y ). For this analysis, we assume an exponential correlation function, i.e., ρ ( || s − s ∗ || ; φ y ) =exp( − φ y || s − s ∗ || ), where || s − s ∗ || is the Euclidean distance between locations s and s ∗ in kilometers (km).This specification for the spatial correlation function requires the modeler to assume the spatially structurederror variability to be stationary and isotropic, meaning that spatial variability left after accounting forthe trend does not depend on location or direction. This is true for the spatial random effects definedin subsequent sections as well. To ease interpretation of the φ y estimates, corresponding effective rangeestimates, labeled er y , are presented in the results tables. We define er y as the distance (km) where thespatial correlation between locations drops to 0.05. See Gelfand et al. (2004) for specifics on how to calculateeffective spatial ranges using coregionalization models. + Tree Cover
Introducing a spatial random effect in model (2) defines the
Spatial + Tree Cover candidate written as y ( s ) = β y + β y x ( s ) + w y ( s ) + (cid:15) y ( s ) , (4)with all terms defined previously. A bivariate coregionalization model for the joint prediction of square-root transformed AGB and lidar-derived 90 th percentile height is (cid:20) y ( s ) z ( s ) (cid:21) = (cid:20) β y β z (cid:21) + A (cid:20) u y ( s ) u z ( s ) (cid:21) + (cid:20) (cid:15) y ( s ) (cid:15) z ( s ) (cid:21) , (5)where y ( s ) , β y and (cid:15) y ( s ) are defined previously. z ( s ) is 90 th percentile height recorded at location s and,because no other covariates appear in the sub-model for z ( s ), β z approximates the overall mean of the12ampled 90 th percentile heights. Here, u y ( s ) and u z ( s ) are modeled as GP s with zero means and spatialcorrelation functions (rather than covariance functions as with w y ( s ) in Spatial and
Spatial + Tree Cover models), i.e., u y ( s ) ∼ GP (0 , ρ ( s , s ∗ ; φ y )) and u z ( s ) ∼ GP (0 , ρ ( s , s ∗ ; φ z )), where φ y and φ z are spatialdecay parameters. We again assume an exponential correlation function for both random effects, i.e., ρ ( || s − s ∗ || ; φ y ) = exp( − φ y || s − s ∗ || ) and ρ ( || s − s ∗ || ; φ z ) = exp( − φ z || s − s ∗ || ). Effective range estimates( er y and er z ) are again presented in the results tables to facilitate interpretation. The matrix K = AA (cid:48) models the cross-covariance between the spatial random effects u y ( s ) and u z ( s ). The 2 x 2 parameter matrix A is the lower triangular square-root of K . The diagonal and off-diagonal elements of K capture the spatialprocesses variances and covariances, respectively. We assume the errors associated with the lidar sub-modelare independent and normally distributed, i.e., (cid:15) z iid ∼ N ( , τ z I ), where (cid:15) z = ( (cid:15) z ( s ) , (cid:15) z ( s ) , ..., (cid:15) z ( s m )) (cid:48) and τ z is a variance parameter to be estimated. + Tree Cover
The
Coregionalization model described in model (5) can be extended to include wall-to-wall Landsat-derived tree cover information as such, (cid:20) y ( s ) z ( s ) (cid:21) = (cid:20) β y + β y x ( s ) β z + β z x ( s ) (cid:21) + A (cid:20) u y ( s ) u z ( s ) (cid:21) + (cid:20) (cid:15) y ( s ) (cid:15) z ( s ) (cid:21) , (6)where all terms with the exception of β z have been previously defined. The lidar sub-model intercept, β y ,and regression slope parameter, β y , together describe the linear relationship between the Landsat-basedtree cover product, x ( s ), and z ( s ). For all six candidate models, a Bayesian paradigm of statistical inference was pursued, which requiredus to specify prior distributions for all model parameters. Then inference proceeded by sampling from theposterior distribution of the parameters. We set prior distributions to be as vague as possible to minimizetheir influence on posterior inference. For the regression intercept and slope parameters, i.e., β y , β y , β z and β z , we assumed a N (0 ,
10 000). Any spatial and non-spatial variance components, i.e., σ y , τ y and τ z ,were assigned an inverse Gamma prior IG ( a, b ). The a hyperparameter was set equal to 2, which resultsin a prior distribution mean equal to b and infinite variance. The b hyperpriors were determined usingpreliminary semi-variograms fit to the residuals of the non-spatial models (1) and (2). The spatial decayparameters, φ y and φ z were assigned uniform priors with support over the geographic range of the studyareas. The matrix K was assigned an inverse Wishart prior with hyperparameter degrees of freedom setto 2 and diagonal scale matrix equaling 0 . I , where I is a 2 x 2 identity matrix. Algorithms for efficientestimation of parameters for all six candidate models are detailed in Banerjee et al. (2014) and Finleyet al. (2014a). All six models were fitted using the spBayes package written for the R statistical computingenvironment (Finley et al., 2015). After collecting a sufficient number of samples from the posterior distribution of a model’s parametersvia MCMC (following typical sampling and convergence diagnostics in Gelman et al. (2013)), composi-tion sampling was used to obtain samples from the posterior predictive distribution (PPD) of square-rootAGB density at all grid cell locations within the study area (Banerjee et al., 2014). Then each sam-ple from the PPD was back-transformed (squared) to approximate spatially explicit AGB density PPDs(˜ y = (˜ y ( s ) , ˜ y ( s ) , . . . , ˜ y ( s n )) (cid:48) , where n is the number of prediction units, e.g., grid cells). Useful sum-maries of ˜ y ( s ) for each grid cell, such as, median, standard deviation, or credible interval width, can bemapped to examine the spatial distribution of AGB density and associated uncertainty. ˜ y can also be sum-marized to estimate average AGB densities for arbitrary areal units. For example, to estimate mean AGBdensity for the full spatial domain we need to integrate the PPD of AGB density over the entire study region,i.e., ˜ y ( D ) = |D| − (cid:82) D ˜ y ( s )d s where |D| is the area of the domain D . We can approximate this integral via13CMC integration using a fine grain systematic grid of ˜ y ( s ) samples, i.e., ˜ y ( D ) ≈ |D| − (cid:80) L l =1 ˜ y ( s l ), where L is the number of equally-spaced prediction grid points covering domain D . To convert ˜ y ( D ) to a PPD fortotal AGB, we multiply each sample by |D| . This process also can be used to generate mean AGB densityPPDs for sub-domains by integrating grid cell-level PPDs over redefined areal blocks, i.e., |B| − (cid:82) B ˜ y ( s )d s ,where B is a sub-region of D . See Banerjee et al. (2014, Chap. 7) for a more thorough discussion on sum-marizing PPDs over areal units. Total (mean) AGB estimates presented in the results tables were obtainedby calculating the median of the total (mean) AGB PPD for the study area (labeled Est ). The standarddeviation of the total (mean) AGB PPD is presented in the results tables as well and serves as a model-baseduncertainty estimate for total (mean) AGB (labeled SD ). Relative standard deviations (labeled RSD in theresults tables) were calculated as
SD/Est ∗ Cokriging interpolation is commonly used in situations where multiple response variables need to bepredicted across space and those responses are posited to co-varying across the domain. A cokriging approachis able to spatially predict the responses of interest accounting for spatial autocorrelation. Additionally,cokriging variance estimates can be calculated at prediction locations. Further, block-cokriging can generatesmall- and large-area estimates of the mean of a response. Given the similar capabilities of cokriging tothe coregionalization framework examined here and the documented use of cokriging in forestry and remotesensing, it is useful to conduct a comparison of the two approaches to potentially identify the merits andshortcomings of each (Meng et al., 2009; Mutanga & Rugege, 2006; Tsui et al., 2013; Wang et al., 2009).Our interest is in understanding if either method produces more accurate point predictions and whetheruncertainty intervals produced by each approach exhibit intended coverage.Cokriging interpolation of the BCEF, CPCRW, TNWR and TVSF data sets was implemented usingthe gstat package for the R statistical computing environment, similarly to the cokriging estimation ap-proach detailed in Tsui et al. (2013). Specifically, the variogram function was used to develop semi- andcross-variograms for the square-root AGB and 90 th percentile height responses. Subsequent exponentialsemi- and cross-variogram models were fitted using the gstat package’s fit.lmc function. The fit.lmc function forces restrictions on semi- and cross-variogram parameter estimates to ensure resulting cokrigingvariance-covariance matrices are positive definite. Positive definiteness is a requirement for resulting cok-riging prediction variances to be valid, i.e., variances be non-negative. One particularly strong restrictionimposed by fit.lmc is that ranges for all semi- and cross-variogram models need to be equal. This is not arestriction for cokriging in general but approaches for building valid cross-variograms relaxing these typesof restrictions can be complex and are not typically used in applications of cokriging in practice (Hoef &Barry, 1998). Once the models were fitted, the gstat package’s predict.krige function was used to garnerpredictions and associated prediction variances for unobserved locations. To examine grid cell-level predictive performance of the six candidate models, for each study site, a 10-fold holdout design was constructed by randomly assigning AGB plot observations to 10 approximately equalsize groups. Square-root transformed AGB for each holdout group was sequentially predicted given modelparameters estimated using data in the remaining nine groups. 10-fold holdout cross-validation root meansquared prediction error (
CV-RMSE ) was calculated using back-transformed holdout posterior predictedmedians and observed AGB for each model at all four sites. The model with the lowest
CV-RMSE wasconsidered the best grid cell-level predictor. The holdout predictions used for
CV-RMSE calculation werealso used to generate the holdout residual semi-variograms presented in Fig. 3.A similar 10-fold cross-validation procedure was used to obtain out-of-sample predictions using cokrigingfor the four sites. To avoid potential bias introduced due to back-transformation of cokriging predictions andvariances, we elected to compare prediction accuracy and uncertainty interval coverage with the
Coregion-alization model on the transformed ( (cid:112)
M g/ha ) scale. Back-transformation in a frequentist setting is muchmore complicated than with the Bayesian models (Stow et al., 2006). Empirical 95% coverage probabilitywas assessed by dividing the number of 95% prediction uncertainty intervals containing y ( s ) by the total14umber of observations at the site. For intervals exhibiting proper coverage, the empirical 95% coverageprobability should be near 0.95.
4. Results and Discussion
Tables 1, 2, 3 and 4 present parameter posterior distribution summaries and prediction accuracy estimatesfor the six candidate models applied to the BCEF, CPCRW, TNWR and TVSF datasets. At all four sites,results show the
Null models to be more precise estimators of total (mean) AGB (labeled
Est in Tables 1, 2,3 and 4) than their
Spatial counterpart models, evidenced by lower total (mean) AGB standard deviations(labeled SD in Tables 1, 2, 3 and 4). However, the Null models also show dramatically higher
CV-RMSE accuracy assessments compared to the
Spatial models at all sites, suggesting that the models includingspatial random effects produced more accurate predictions.The increased accuracy coupled with apparent decreased precision of the
Spatial model predictionshighlight the role of distributional assumptions concerning the SD and CV-RMSE estimates. The
CV-RMSE metric employed here does not rely on distributional assumptions concerning stochastic model componentsbecause it only assesses the ability of the PPD median point estimate to approximate observed AGB at plotlocations—see Efron & Gong (1983) for discussion on robustness of cross-validation strategies to modelingassumptions. On the other hand, model-based summaries of PPDs, such as, the SD precision metricpresented in Tables 1, 2, 3 and 4, heavily rely on modeling assumptions. Strong violations of distributionalassumptions made during model-fitting can lead to unrealistic uncertainty estimates. Looking at the Null model holdout residual semi-variograms in Fig. 3, we see that, at all four sites, there is strong residualspatial autocorrelation. This means that the
Null models likely violate the iid error assumption imposedto conduct model-based inference about any parameters or predictions using standard deviations or otherposterior distribution summaries. The SD metrics for total (mean) AGB resulting from the Null models areunderestimated because the errors are falsely assumed to be independent of one another (Griffith, 2005).The
Spatial model semi-variograms, however, show no signs of residual spatial structure, indicating that w ( s ) has effectively absorbed any extraneous spatial variability (Fig. 3). Relaxing the iid error assumptionby incorporating spatial random effects allows for modelers to better interpret model-based uncertaintymeasures. Total (mean) AGB SD estimates garnered via the Spatial models are more reliable than the
Null model SD s.A similar phenomenon is observed when comparing the Tree Cover and
Spatial + Tree Cover modelpredictive performance measures, although differences were smaller. Comparing the
Null and
Tree Cover holdout residual semi-variograms in Fig. 3 shows that introducing the Landsat-derived tree cover variableabsorbs a portion of the spatially structured variability in AGB at all four sites. However, spatial structurestill exists in the holdout residuals of the
Tree Cover models. We argue that, due to positive residualspatial autocorrelation, SD estimates are underestimated for the Tree Cover models as well. Comparisonsbetween the univariate non-spatial (
Null and
Tree Cover ) and spatial (
Spatial and
Spatial + Tree Cover )models underscore the need to adhere to posited distributional assumptions concerning stochastic modelcomponents when attempting to interpret model-dependent measures of predictive performance, includinggrid cell-level and areally integrated PPDs.
At all four sites, the lidar-informed
Coregionalization models outperformed their counterpart lidar-free
Spatial models by producing lower
CV-RMSE accuracy assessments and total (mean) AGB SD metrics. Sim-ilarly, the Coregionalization + Tree Cover models produced better accuracy measures than the
Spatial + TreeCover models. Prediction gains resulting from the incorporation of lidar information are not surprisingconsidering that lidar-derived measures of forest height are strongly related to AGB. It is encouraging, how-ever, that this relationship produces appreciable gains in total (mean) AGB estimation precision at suchlow lidar sampling intensities—see Section 3.3 for description of the lidar dataset subsampling. Tables 1, 2,3 and 4 show the correlation between the spatial random effects for the AGB and lidar sub-models (labeled15 able 1:
Bonanza Creek Experimental Forest (BCEF) candidate model parameter estimates, prediction accuracy metricsand total (mean) aboveground biomass (AGB) estimates.
Est = estimated total (mean) AGB with associated 95% credibleinterval in parentheses. SD = posterior predictive standard deviation for AGB total (mean) estimate. RSD = relative standarddeviation for total (mean) AGB estimate (
SD/Est ∗ CV-RMSE = 10-fold holdout cross-validation root mean squaredprediction error accuracy assessment.
Null Tree Cover Spatial Spatial+Tree Cover Coregionalization Coregionalization+Tree Cover P a r a m e t e r p o s t e r i o r q u a n t il e s u mm a r i e s % ( . % , . % ) β y β y — 0.10 (0.08, 0.12) — 0.08 (0.05, 0.10) — 0.08 (0.06, 0.10) β z — — — — 9.35 (7.81, 11.13) 3.16 (1.81, 4.78) β z — — — — — 0.09 (0.08, 0.11) τ y τ z — — — — 2.14 (1.54, 2.89) 2.03 (1.36, 2.82) σ y — — 17.83 (12.42, 23.46) 10.19 (6.34, 14.83) — — K [1 ,
1] — — — — 20.71 (16.50, 25.57) 13.60 (10.40, 17.98) K [1 ,
2] — — — — 25.16 (21.25, 29.50) 14.77 (11.56, 18.66) K [2 ,
2] — — — — 41.39 (36.12, 48.35) 29.51 (25.76, 34.92)cor( u y , u z ) — — — — 0.86 (0.77, 0.92) 0.74 (0.61, 0.86) er y (km) — — 0.40 (0.23, 2.88) 0.32 (0.17, 3.03) 0.78 (0.59, 1.01) 0.36 (0.23, 0.54) er z (km) — — — — 2.37 (1.37, 5.83) 2.15 (1.49, 3.61)TotalAGB (Tg) Est SD RSD
Est SD RSD
CV-RMSE (Mg/ha) 68.58 60.53 47.62 47.44 35.75 36.40
Table 2:
Caribou-Poker Creeks Research Watershed (CPCRW) candidate model parameter estimates, prediction accuracymetrics and total (mean) aboveground biomass (AGB) estimates.
Est = estimated total (mean) AGB with associated 95%credible interval in parentheses. SD = posterior predictive standard deviation for AGB total (mean) estimate. RSD = relativestandard deviation for total (mean) AGB estimate (
SD/Est ∗ CV-RMSE = 10-fold holdout cross-validation root meansquared prediction error accuracy assessment.
Null Tree Cover Spatial Spatial+Tree Cover Coregionalization Coregionalization+Tree Cover P a r a m e t e r p o s t e r i o r q u a n t il e s u mm a r i e s % ( . % , . % ) β y β y — 0.09 (0.08, 0.11) — 0.08 (0.05, 0.09) — 0.05 (0.03, 0.07) β z — — — — 5.25 (3.14, 7.15) 0.66 (-1.28, 2.24) β z — — — — — 0.08 (0.07, 0.10) τ y τ z — — — — 1.82 (1.28, 2.52) 2.15 (1.51, 2.91) σ y — — 8.96 (6.10, 13.58) 2.98 (1.20, 5.41) — — K [1 ,
1] — — — — 9.63 (6.83, 12.95) 4.86 (2.59, 8.59) K [1 ,
2] — — — — 14.03 (10.56, 18.38) 7.10 (4.74, 10.03) K [2 ,
2] — — — — 26.32 (20.46, 34.39) 14.60 (11.61, 19.55)cor( u y , u z ) — — — — 0.89 (0.81, 0.94) 0.85 (0.71, 0.92) er y (km) — — 1.28 (0.55, 2.91) 0.33 (0.10, 2.19) 1.96 (1.40, 2.95) 1.20 (0.76, 1.92) er z (km) — — — — 2.49 (1.66, 6.35) 2.31 (1.32, 6.01)TotalAGB (Tg) Est SD RSD
Est SD RSD
CV-RMSE (Mg/ha) 41.55 27.41 23.88 24.06 20.78 21.19 able 3: Tetlin National Wildlife Refuge (TNWR) candidate model parameter estimates, prediction accuracy metrics andtotal (mean) aboveground biomass (AGB) estimates.
Est = estimated total (mean) AGB with associated 95% credible intervalin parentheses. SD = posterior predictive standard deviation for AGB total (mean) estimate. RSD = relative standarddeviation for total (mean) AGB estimate (
SD/Est ∗ CV-RMSE = 10-fold holdout cross-validation root mean squaredprediction error accuracy assessment.
Null Tree Cover Spatial Spatial+Tree Cover Coregionalization Coregionalization+Tree Cover P a r a m e t e r p o s t e r i o r q u a n t il e s u mm a r i e s % ( . % , . % ) β y β y — 0.07 (0.05, 0.10) — 0.06 (0.04, 0.09) — 0.07 (0.05, 0.09) β z — — — — 3.09 (2.38, 3.68) -0.39 (-1.00, 0.14) β z — — — — — 0.08 (0.07, 0.09) τ y τ z — — — — 1.96 (1.39, 2.70) 2.24 (1.63, 2.95) σ y — — 12.90 (8.57, 21.39) 8.59 (5.46, 12.98) — — K [1 ,
1] — — — — 17.70 (12.69, 24.90) 11.50 (8.12, 15.24) K [1 ,
2] — — — — 13.58 (11.27, 16.92) 8.48 (7.05, 10.16) K [2 ,
2] — — — — 15.31 (13.57, 17.77) 9.05 (7.80, 10.73)cor( u y , u z ) — — — — 0.83 (0.75, 0.89) 0.83 (0.76, 0.89) er y (km) — — 0.36 (0.19, 0.72) 0.28 (0.15, 0.83) 0.57 (0.44, 0.74) 0.39 (0.30, 0.50) er z (km) — — — — 5.23 (2.96, 9.71) 5.19 (2.80, 19.17)TotalAGB (Tg) Est SD RSD
Est SD RSD
CV-RMSE (Mg/ha) 51.89 44.70 37.37 38.01 34.03 30.89
Table 4:
Tanana Valley State Forest (TVSF) candidate model parameter estimates, prediction accuracy metrics and total(mean) aboveground biomass (AGB) estimates.
Est = estimated total (mean) AGB with associated 95% credible interval inparentheses. SD = posterior predictive standard deviation for AGB total (mean) estimate. RSD = relative standard deviationfor total (mean) AGB estimate (
SD/Est ∗ CV-RMSE = 10-fold holdout cross-validation root mean squared predictionerror accuracy assessment.
Null Tree Cover Spatial Spatial+Tree Cover Coregionalization Coregionalization+Tree Cover P a r a m e t e r p o s t e r i o r q u a n t il e s u mm a r i e s % ( . % , . % ) β y β y — 0.09 (0.07, 0.11) — 0.07 (0.04, 0.10) — 0.08 (0.06, 0.09) β z — — — — 8.35 (7.68, 9.06) 3.03 (2.50, 3.68) β z — — — — — 0.08 (0.07, 0.09) τ y τ z — — — — 2.61 (2.03, 3.20) 2.48 (2.05, 2.96) σ y — — 18.60 (14.04, 24.08) 13.15 (9.67, 18.23) — — K [1 ,
1] — — — — 23.80 (18.93, 29.98) 17.38 (13.50, 21.67) K [1 ,
2] — — — — 22.61 (19.30, 26.37) 16.48 (14.50, 18.90) K [2 ,
2] — — — — 40.00 (36.73, 43.89) 27.86 (25.89, 29.63)cor( u y , u z ) — — — — 0.73 (0.67, 0.80) 0.75 (0.70, 0.81) er y (km) — — 0.56 (0.35, 1.28) 0.42 (0.26, 0.87) 0.73 (0.60, 0.91) 0.47 (0.40, 0.58) er z (km) — — — — 7.10 (5.14, 11.17) 4.49 (3.39, 6.72)TotalAGB (Tg) Est SD RSD
Est SD RSD
CV-RMSE (Mg/ha) 72.37 63.46 42.68 42.15 31.00 28.95 ig. 3: Holdout residual semi-variograms for the six candidate models at all four study sites. or ( u y , u z ) in Tables 1, 2, 3 and 4). Median point estimates of correlation range between 0.73 and 0.89indicating strong correspondence between square-root transformed AGB density and 90 th percentile height.Since the relationship between the AGB and lidar responses is so substantial, considerable improvementsin predictive performance and total (mean) AGB estimation precision are realized. At all sites, relativestandard deviations (labeled RSD in Tables 1, 2, 3 and 4) for the
Coregionalization models are nearly halfthat of the
Spatial models.At the sites with sparser field and lidar sampling, i.e., TVSF and TNWR, we see further gains in gridcell-level prediction accuracy and total (mean) AGB certainty when the Landsat-based tree cover product isincluded—evidenced by
Coregionalization + Tree Cover models having lower
RSD and
CV-RMSE estimatesthan counterpart
Coregionalization frameworks. Examining effective range estimates (labeled er y and er z inTables 1, 2, 3 and 4) offers insight concerning this phenomenon. Effective range estimates help us understandhow far-reaching the underlying spatial dependence structure is in the responses, i.e., y ( s ) and z ( s ), afteraccounting for covariates. The longer the effective range, the larger the proximate neighborhood each gridcell draws information from. The er z median point estimates for the TNWR and TVSF Coregionalization and
Coregionalization + Tree Cover models extend less than the 10 km separation distance between lidarflight-lines at those sites (4.49 km < er z < .
10 km). However, the effective range estimates for thecorresponding models exceed the average flight-line separation distances of < . < er z < .
41 km). On average, grid cell-level predictions at BCEF and CPCRW are borrowingmore proximate lidar data to inform prediction because there is more available within the effective range ofspatial dependence at these sites. Since field plots are also closer together at CPCRW and BCEF comparedto TVSF and TNWR, grid cell predictions are potentially tapping into more information provided by nearbyAGB observations as well. Any AGB explanatory power offered by the Landsat tree cover product is alreadyexplained by simply borrowing strength from nearby field and lidar data to inform prediction at BCEF andCPCRW. However, since field and lidar data are much farther apart at TVSF and TNWR, we do see benefitfrom the wall-to-wall information provided by the Landsat-based tree cover product. + Tree Cover mapped predictions
Figs. 4, 5, 6 and 7 show maps of PPD median point estimates and standard deviations for AGB densityusing the
Coregionalization and
Coregionalization + Tree Cover models. The flight-lines where lidar wascollected are easily discernible in the
Coregionalization model prediction maps (Figs. 4 and 5). We seedifferentiation of AGB within lidar strips and near field plot locations. As one moves away from regionswhere information was collected, predicted AGB densities retreat to the global mean. Additionally, gridcell-level uncertainty maps show a trend of higher predictive precision in close proximity to field and lidarobservations compared to areas farther away. Even though the
Coregionalization models provide reliableestimates of total (mean) AGB and show favorable
CV-RMSE accuracy assessments, the resulting predictionmaps leave something to be desired. Having no information outside the flight-lines, the coregionalizationmodels only rely on borrowing strength from nearby observations to inform prediction between strips. Theaddition of the wall-to-wall Landsat tree cover product allows for better differentiation of median pointpredictions outside lidar strips, leading to maps that are potentially more useful for assessing the spatialdistribution of AGB across the landscape. At TVSF and TNWR,
CV-RMSE estimates suggest that theaddition of the Landsat tree cover product leads to better map accuracy. However, at BCEF and CPCRW,the inclusion of Landsat tree cover appears to only offer cosmetic adjustments to resulting maps, not actualincreases in prediction accuracy. It seems that the intensity of the field and lidar sample at BCEF andCPCRW mixed with constructing models to borrow information from neighboring locations offered morepredictive advantage than simply including Landsat tree cover as a covariate.
Fig. 8 shows estimated mean AGB densities for the 11 watershed units at CPCRW predicted using the
Coregionalization + Tree Cover model with a permafrost polygon layer overlaid. Average AGB density wasnegatively correlated with permafrost presence at the watershed unit level. By appropriately summarizingthe model’s PPD surface, it is possible to obtain an estimate of the correlation between AGB density and19ermafrost proportion at the watershed unit scale while accounting for uncertainty in the AGB densitypredictions. By calculating the correlation between the proportion of the watershed unit covered in per-mafrost with each MCMC sample of average AGB density, we can obtain a PPD summary characterizingtheir degree of correlation. The median correlation point estimate representing the relationship betweenpermafrost and watershed-unit AGB density is − . − .
587 and − . (a) BCEF: Coregionalization predictions (b)
BCEF: Coregionalization prediction standard deviations (c)
CPCRW: Coregionalization predictions (d)
CPCRW: Coregionalization prediction standard deviations
Fig. 4:
Mapped grid cell-level predictions of aboveground biomass density (a and c) and associated standard deviations (band d) (Mg/ha) using the
Coregionalization model for Bonanza Creek Experimental Forest (BCEF) and Caribou-Poker CreeksResearch Watershed (CPCRW). The red polygon boundary shown in c and d defines the area burned during the 2004 BoundaryFire at CPCRW. . Comparing Coregionalization model and cokriging predictions
Table 5 presents holdout prediction accuracies and 95% coverage probabilities for the cokriging inter-polations and
Coregionalization models at BCEF, CPCRW, TNWR and TVSF. The results in Table 5indicate that cokriging and
Coregionalization holdout prediction accuracies were similar at all four studysites, although prediction accuracy estimates were slightly better for the
Coregionalization models at BCEF,CPCRW and TNWR. This is likely due to the restrictions imposed on the fitted variogram models neces-sary for valid cokriging prediction variances. In contrast, the
Coregionalization framework estimated usingBayesian inference ensures positive definiteness by construction, meaning that, as long as the covariancefunctions used for the spatial correlation elements and the K variance-covariance matrix are positive def-inite, the model is estimable. Positive definiteness for the spatial random effect variances and covariancesis satisfied through appropriate prior specification, e.g., the IW prior on K is a distribution of only pos-itive definite matrices and the spatial decay priors only consider positive values. The increased flexibilityin modeling spatial autocorrelation provided by the Coregionalization framework estimated using Bayesianinference may be leading to slightly improved prediction over the gstat implementation of cokriging atBCEF, CPCRW and TNWR. As shown in Table 5, the empirical 95% coverage probabilities for the
Core-gionalization model holdout uncertainty intervals better match the intended 95% coverage than the cokriging (a)
TNWR: Coregionalization predictions (b)
TNWR: Coregionalization prediction standard deviations (c)
TVSF: Coregionalization predictions (d)
TVSF: Coregionalization prediction standard deviations
Fig. 5:
Mapped grid cell-level predictions of aboveground biomass density (a and c) and associated standard deviations (b andd) (Mg/ha) using the
Coregionalization model for Tetlin National Wildlife Refuge (TNWR) and Tanana Valley State Forest(TVSF).
6. Conclusions and Next Steps
The goal of this analysis was to develop and test the performance of a statistical modeling framework thatcan 1) incorporate partial coverage lidar data and wall-to-wall Landsat products to improve AGB densityprediction; and 2) accommodate spatially structured error, thereby allowing for more reliable model-basedcharacterization of uncertainty and improved prediction. Through model comparison we were able to show (a)
BCEF: Coregionalization+Tree Cover predictions (b)
BCEF: Coregionalization+Tree Cover prediction standarddeviations (c)
CPCRW: Coregionalization+Tree Cover predictions (d)
CPCRW: Coregionalization+Tree Cover prediction stan-dard deviations
Fig. 6:
Mapped grid cell-level predictions of aboveground biomass density (a and c) and associated standard deviations (band d) (Mg/ha) using the
Coregionalization + Tree Cover model for Bonanza Creek Experimental Forest (BCEF) and Caribou-Poker Creeks Research Watershed (CPCRW). The red polygon boundary shown in c and d defines the area burned during the2004 Boundary Fire at CPCRW. able 5: Cokriging and
Coregionalization holdout prediction accuracies and 95% coverage probabilities. To avoid potentialbias introduced due to back-transformation of cokriging predictions and variances, prediction and uncertainty interval coverageis compared on the transformed ( (cid:112)
Mg/ha ) scale.
CV-RMSE ( (cid:112)
M g/ha ) 95% Coverage Probabilitycokriging coregionalization cokriging coregionalizationBCEF 2.23 2.14 0.739 0.947CPCRW 1.75 1.63 0.651 0.946TNWR 2.04 1.96 0.722 0.963TVSF 1.94 1.96 0.620 0.945that a coregionalization framework can effectively be used to couple sampled lidar and field data to improvegrid cell-level AGB density prediction accuracy and increase confidence in total AGB estimates. Rigorousmodel comparison also demonstrated the adverse effects of spatial autocorrelation and how including appro-priately specified random effects within a Bayesian hierarchical framework can absorb spatial dependence (a)
TNWR: Coregionalization+Tree Cover predictions (b)
TNWR: Coregionalization+Tree Cover prediction standarddeviations (c)
TVSF: Coregionalization+Tree Cover predictions (d)
TVSF: Coregionalization+Tree Cover prediction standarddeviations
Fig. 7:
Mapped grid cell-level predictions of aboveground biomass density (a and c) and associated standard deviations (b andd) (Mg/ha) using the
Coregionalization + Tree Cover model for Tetlin National Wildlife Refuge (TNWR) and Tanana ValleyState Forest (TVSF).
23n the response variable not accounted for by the covariates.An examination of
Coregionalization and
Coregionalization + Tree Cover model mapped predictions re-vealed benefits of including wall-to-wall information to better visualize spatial variability in AGB. We alsosaw that greater differentiation of AGB median point estimates provided by Landsat-based tree cover datadoes not necessarily lead to improved grid cell-level predictive performance. It can be the case that simplyconstructing models to leverage proximate observations may improve prediction more than supplementingauxiliary data, especially when sites are intensively sampled and ancillary information is only weakly relatedto the response variable.We demonstrated the ability to summarize PPDs from the
Coregionalization + Tree Cover model togenerate small-area estimates of mean AGB density at CPCRW. Small-area estimates are a challenge fortraditional estimation methods (Breidenbach & Astrup, 2012; Goerndt et al., 2011), and have also limitedthe ability to investigate fine-scale spatial heterogeneity within larger management units such as the TVSFor TNWR. In contrast, our results demonstrate how small-area integrated PPDs can be used to examinecorrelation with other information about sub-domains while accounting for prediction uncertainty. Relatingwatershed unit-level AGB density to permafrost and fire history polygon layers uncovered possible driversof AGB variability among watershed units at CPCRW while explicitly accounting for uncertainty of AGBdensity predictions. This exercise shows the inferential potential of AGB PPD products generated usingBayesian hierarchical spatial models beyond simply mapping AGB density.Upcoming satellite lidar missions, e.g., ICESat-2 and GEDI LiDAR, are designed to collect data alongorbital transects but not provide complete spatial coverage. Results from this study suggest that it may bepossible to implement a coregionalization approach to augment the FIA program’s network of permanentsample plots over the contiguous United States with sampled space-based lidar to improve forest inventoryestimates. Further, Landsat products can be incorporated to potentially improve point-level prediction
Fig. 8:
Left figure maps watershed unit-level mean aboveground biomass (AGB) (Mg/ha) with associated standard deviationsin parentheses using the
Coregionalization + Tree Cover model for Caribou-Poker Creeks Research Watershed. Solid polygonboundaries delineate watershed units. Translucent purple polygon identifies area covered in permafrost. The black dashedline delineates the border of the 2004 Boundary Fire. The red outlined watershed units are identified as the two units withthe highest proportion of area burned during the 2004 Boundary Fire. The right figure shows a scatter plot highlighting therelationship between the proportion of permafrost in a watershed unit and AGB density. The dashed trend-line shows thegeneral relationship between AGB density and permafrost proportion. The points shown in red correspond to the red outlinedwatershed units in the left figure.
Coregionalization model took twodays to complete a single MCMC chain of 50 000 iterations on a Linux workstation equipped with aneight-core processor leveraging threaded BLAS and LAPACK C++ libraries ( and ). When lidar and/or field observations exceed 10 000, fitting this class of multi-variate spatial models within a Bayesian paradigm becomes intractable on typical desktop computers dueto computational difficulties related to inverting massive matrices. Recent advances in GP modeling theoryand computation have unearthed new approaches for estimating spatial random effects that circumvent thelimitations imposed by the necessity to invert large matrices (Datta et al., 2016a,b; Finley et al., 2017).Implementing Nearest-neighbor Gaussian process priors ( N N GP ) in place of traditional GP s alleviates theneed to invert large covariance matrices, allowing modelers to effectively estimate complex Bayesian spatialmodels using many remote sensing and field observations on standard workstations. In the future, we willinvestigate the approximation of GP s with N N GP s to allow for the inclusion of more field and remotesensing samples within coregionalization frameworks to further improve grid cell-level prediction and arealestimation certainty. These techniques will allow broader applicability of robust modeling techniques toother disciplines, based on the large data volumes typically associated with spatially extensive airborne andsatellite remote sensing products. We also plan to extend this class of multivariate spatial models to allow forthe simultaneous prediction of multiple forest variables, e.g., tree density, basal area and AGB, while leverag-ing spatial cross-correlations between all responses. By jointly predicting many forest inventory parameterswe will be able to preserve the inherent relationships between them and use those intrinsic correlations toaid spatial prediction and subsequent areal estimation of all included forest inventory parameters.
7. Acknowledgments
The research presented in this study was partially supported by NASA’s Arctic-Boreal VulnerabilityExperiment (ABoVE) and Carbon Monitoring System (CMS) grants (proposal 13-CMS13-0006 funded viasolicitation NNH13ZDA001N-CMS and proposal 15-CMS13-0012 funded via solicitation NNH15ZDA001N-CMS). Additional support was provided by the United States Forest Service Pacific Northwest ResearchStation. Andrew Finley was supported by National Science Foundation (NSF) DMS-1513481, EF-1137309,EF-1241874, and EF-1253225 grants. 25 eferences
Alaska Department of Natural Resources (2016). Alaska’s state forests. http://forestry.alaska.gov/stateforests . Accessed:10-16-2016.Amatya, D., Campbell, J., Wohlgemuth, P., Elder, K., Sebestyen, S., Johnson, S., Keppeler, E., Adams, M., Caldwell, P., &Misra, D. (2016). Hydrological processes of reference watersheds in experimental forests, USA., in: Amatya, D., Williams,T., Bren, L., & de Jong, C. (Eds.), Forest Hydrology: Processes, Management, and Applications. CABI Publishers, UK, pp.219–239.Andersen, H.E., Strunk, J., & Temesgen, H. (2011). Using airborne light detection and ranging as a sampling tool for estimatingforest biomass resources in the upper Tanana Valley of interior Alaska.
Western Journal of Applied Forestry , 26, 157–164.Asner, G., Hughes, R., Varga, T., Knapp, D., & Kennedy-Bowdoin, T. (2009). Environmental and biotic controls overaboveground biomass throughout a tropical rain forest.
Ecosystems , 12, 261–278.Babcock, C., Finley, A.O., Bradford, J.B., Kolka, R., Birdsey, R., & Ryan, M.G. (2015). Lidar based prediction of forestbiomass using hierarchical models with spatially varying coefficients.
Remote Sensing of Environment , 169, 113–127.Babcock, C., Finley, A.O., Cook, B.D., Weiskittel, A., & Woodall, C.W. (2016). Modeling forest biomass and growth: Couplinglong-term inventory and lidar data.
Remote Sensing of Environment , 182, 1–12.Babcock, C., Matney, J., Finley, A., Weiskittel, A., & Cook, B. (2013). Multivariate spatial regression models for predictingindividual tree structure variables using lidar data.
IEEE Journal of Selected Topics in Applied Earth Observations andRemote Sensing , 6, 6–14.Banerjee, S., Carlin, C., & Gelfand, A. (2014). Hierarchical Modeling and Analysis for Spatial Data. Chapman & Hall/CRCMonographs on Statistics & Applied Probability. second ed., Chapman & Hall/CRC.Barrett, T.M., & Gray, A.N. (2011). Potential of a national monitoring program for forests to assess change in high-latitudeecosystems.
Biological Conservation , 144, 1285–1294.Bechtold, W.A., & Patterson, P.L. (2005). The Enhanced Forest Inventory and Analysis Program: National Sampling Designand Estimation Procedures. US Department of Agriculture Forest Service, Southern Research Station Asheville, NorthCarolina.Berliner, L.M. (1996). Hierarchical bayesian time series models, in: Maximum Entropy and Bayesian Methods. Springer, pp.15–22.Bolton, D.K., Coops, N.C., & Wulder, M.A. (2013). Measuring forest structure along productivity gradients in the Canadianboreal with small-footprint lidar.
Environmental monitoring and assessment , 185, 6617–6634.Bonanza Creek LTER (2016a). Bonanza Creek Experimental Forest. .Accessed: 10-16-2016.Bonanza Creek LTER (2016b). Caribou-Poker Creeks Research Watershed. . Accessed: 10-16-2016.Bradshaw, C.J., & Warkentin, I.G. (2015). Global estimates of boreal forest carbon stocks and flux.
Global and PlanetaryChange , 128, 24–30.Breidenbach, J., & Astrup, R. (2012). Small area estimation of forest attributes in the Norwegian national forest inventory.
European Journal of Forest Research , 131, 1255–1267.Chapin, F., & Hollingsworth, J. (2006). GIS file of permafrost distribution in Caribou-Poker Creeks Research Watershed,Bonanza Creek LTER - University of Alaska Fairbanks. .CMS (2010). NASA carbon monitoring system. http://carbon.nasa.gov . Accessed: 8-11-2016.Cochran, W. (1946). Relative accuracy of systematic and stratified random samples for a certain class of populations.
TheAnnals of Mathematical Statistics , 17, 164–177.Cochran, W. (1977). Sampling Techniques. Wiley Series in Probability and Mathematical Statistics: Applied Probability andStatistics, Wiley.Cook, B., Corp, L., Nelson, R., Middleton, E., Morton, D., McCorkel, J., Masek, J., Ranson, K., Ly, V., & Montesano, P.(2013). NASA Goddard’s lidar, hyperspectral and thermal (G-LiHT) airborne imager.
Remote Sensing , 5, 4045–4066.Cressie, N. (1990). The origins of kriging.
Mathematical Geology , 22, 239–252.Cressie, N. (1993). Statistics for Spatial Data. Wiley series in probability and mathematical statistics: Applied probabilityand statistics, John Wiley & Sons.Cressie, N., Calder, C.A., Clark, J.S., Ver Hoef, J.M.V., & Wikle, C.K. (2009). Accounting for uncertainty in ecologicalanalysis: The strengths and limitations of hierarchical statistical modeling.
Ecological Applications , 19, 553–570.Cressie, N.A.C., & Wikle, C.K. (2011). Statistics for Spatio-Temporal Data. Wiley series in probability and statistics, Hoboken,N.J. Wiley.Datta, A., Banerjee, S., Finley, A.O., & Gelfand, A.E. (2016a). Hierarchical nearest-neighbor gaussian process models for largegeostatistical datasets.
Journal of the American Statistical Association , 111, 800–812.Datta, A., Banerjee, S., Finley, A.O., & Gelfand, A.E. (2016b). On nearest-neighbor gaussian process models for massivespatial data.
Wiley Interdisciplinary Reviews: Computational Statistics , 8, 162–171.Deo, R.K., Russell, M.B., Domke, G.M., Andersen, H.E., Cohen, W.B., & Woodall, C.W. (2017). Evaluating site-specific andgeneric spatial models of aboveground forest biomass based on landsat time-series and lidar strip samples in the eastern usa.
Remote Sensing , 9.Diggle, P., & Ribeiro, P. (2007). Model-based Geostatistics. Springer Series in Statistics, Springer.Ediriweera, S., Pathirana, S., Danaher, T., & Nichols, D. (2014). Estimating above-ground biomass by fusion of lidar andmultispectral data in subtropical woody plant communities in topographically complex terrain in north-eastern australia.
Journal of Forestry Research , 25, 761–771. fron, B., & Gong, G. (1983). A leisurely look at the bootstrap, the jackknife, and cross-validation. The American Statistician ,37, 36–48.Ene, L.T., Næsset, E., Gobakken, T., Gregoire, T.G., St˚ahl, G., & Nelson, R. (2012). Assessing the accuracy of regionallidar-based biomass estimation using a simulation approach.
Remote Sensing of Environment , 123, 579 – 592.Finley, A.O., Banerjee, S., & Cook, B.D. (2014a). Bayesian hierarchical models for spatially misaligned data in R.
Methods inEcology and Evolution , 5, 514–523.Finley, A.O., Banerjee, S., & Gelfand, A.E. (2015). spBayes for large univariate and multivariate point-referenced spatio-temporal data models.
Journal of Statistical Software , 63, 1–28. URL: .Finley, A.O., Banerjee, S., Weiskittel, A.R., Babcock, C., & Cook, B.D. (2014b). Dynamic spatial regression models forspace-varying forest stand tables.
Environmetrics , 25, 596–609.Finley, A.O., Banerjee, S., Zhou, Y., Cook, B.D., & Babcock, C. (2017). Joint hierarchical models for sparsely sampledhigh-dimensional LiDAR and forest variables.
Remote Sensing of Environment , 190, 149–161.Finley, A.O., Datta, A., Cook, B.C., Morton, D.C., Andersen, H.E., & Banerjee, S. (2017). Applying nearest neighbor gaussianprocesses to massive spatial data sets: Forest canopy height prediction across Tanana Valley Alaska.
ArXiv e-prints , arXiv:1702.00434 .Gauthier, S., Bernier, P., Kuuluvainen, T., Shvidenko, A.Z., & Schepaschenko, D.G. (2015). Boreal forest health and globalchange. Science , 349, 819–822.GEDI (2014). Global ecosystem dynamics investigation lidar. http://science.nasa.gov/missions/gedi/ . Accessed: 8-11-2016.Gelfand, A.E., Schmidt, A.M., Banerjee, S., & Sirmans, C. (2004). Nonstationary multivariate process modeling throughspatially varying coregionalization.
Test , 13, 263–312.Gelfand, A.E., & Smith, A.F. (1990). Sampling-based approaches to calculating marginal densities.
Journal of the Americanstatistical association , 85, 398–409.Gelman, A., Carlin, J., Stern, H., & Rubin, D. (2013). Bayesian Data Analysis, 3rd Edition. Chapman & Hall/CRC Texts inStatistical Science, Chapman & Hall/CRC.Goerndt, M.E., Monleon, V.J., & Temesgen, H. (2011). A comparison of small-area estimation techniques to estimate selectedstand attributes using LiDAR-derived auxiliary variables.
Canadian Journal of Forest Research , 41, 1189–1201.Gregoire, T.G. (1998). Design-based and model-based inference in survey sampling: Appreciating the difference.
CanadianJournal of Forest Research , 28, 1429–1447.Gregoire, T.G., Næsset, E., McRoberts, R.E., St˚ahl, G., Andersen, H.E., Gobakken, T., Ene, L., & Nelson, R. (2016). Statisticalrigor in lidar-assisted estimation of aboveground forest biomass.
Remote Sensing of Environment , 173, 98–108.Gregoire, T.G., St˚ahl, G., Næsset, E., Gobakken, T., Nelson, R., & Holm, S. (2011). Model-assisted estimation of biomass ina lidar sample survey in Hedmark County, Norway.
Canadian Journal of Forest Research , 41, 83–95.Griffith, D.A. (2005). Effective geographic sample size in the presence of spatial autocorrelation.
Annals of the Association ofAmerican Geographers , 95, 740–760.Hansen, M., Potapov, P., Moore, R., Hancher, M., Turubanova, S., Tyukavina, A., Thau, D., Stehman, S., Goetz, S., Loveland,T., Kommareddy, A., Egorov, A., Chini, L., Justice, C., & Townshend, J. (2013). High-resolution global maps of 21st-centuryforest cover change.
Science , 342, 850–853.Hobbs, N., & Hooten, M. (2015). Bayesian Models: A Statistical Primer for Ecologists. Princeton University Press.Hoef, J.M.V., & Barry, R.P. (1998). Constructing and fitting models for cokriging and multivariable spatial prediction.
Journalof Statistical Planning and Inference , 69, 275–294.Hollingsworth, T.N., Johnstone, J.F., Bernhardt, E.L., & Chapin III, F.S. (2013). Fire severity filters regeneration traits toshape community assembly in Alaska’s boreal forest.
PloS one , 8, e56033.Homer, C.G., Dewitz, J.A., Yang, L., Jin, S., Danielson, P., Xian, G., Coulston, J., Herold, N.D., Wickham, J., & Megown, K.(2015). Completion of the 2011 national land cover database for the conterminous United States-representing a decade ofland cover change information.
Photogramm. Eng. Remote Sens , 81, 345–354.Hudak, A.T., Lefsky, M.A., Cohen, W.B., & Berterretche, M. (2002). Integration of lidar and landsat etm+ data for estimatingand mapping forest canopy height.
Remote Sensing of Environment , 82, 397–416.ICESat-2 (2015). Ice, cloud, and land elevation satellite-2. http://icesat.gsfc.nasa.gov/icesat2 . Accessed: 8-11-2016.Jakubowski, M.K., Guo, Q., & Kelly, M. (2013). Tradeoffs between lidar pulse density and forest measurement accuracy.
Remote Sensing of Environment , 130, 245 – 253.Kumar, L., Sinha, P., Taylor, S., & Alqurashi, A.F. (2015). Review of the use of remote sensing for biomass estimation tosupport renewable energy generation.
Journal of Applied Remote Sensing , 9, 097696.Lim, K., Treitz, P., Wulder, M., St-Onge, B., & Flood, M. (2003). Lidar remote sensing of forest structure.
Progress in physicalgeography , 27, 88–106.Lunn, D., Jackson, C., Best, N., Thomas, A., & Spiegelhalter, D. (2012). The BUGS Book: A Practical Introduction toBayesian Analysis. Chapman & Hall/CRC Texts in Statistical Science, Taylor & Francis. URL: https://books.google.com/books?id=Cthz3XMa_VQC .Magnani, F., Mencuccini, M., Borghetti, M., Berbigier, P., Berninger, F., Delzon, S., Grelle, A., Hari, P., Jarvis, P.G., Kolari,P., et al. (2007). The human footprint in the carbon cycle of temperate and boreal forests.
Nature , 447, 849–851.Margolis, H.A., Nelson, R.F., Montesano, P.M., Beaudoin, A., Sun, G., Andersen, H.E., & Wulder, M.A. (2015). Combiningsatellite lidar, airborne lidar, and ground plots to estimate the amount and distribution of aboveground biomass in the borealforest of North America.
Canadian Journal of Forest Research , 45, 838–855.Matheron, G. (1963). Principles of geostatistics.
Economic Geology , 58, 1246–1266.McRoberts, R.E., Næsset, E., & Gobakken, T. (2013). Inference for lidar-assisted estimation of forest growing stock volume.
Remote Sensing of Environment , 128, 268–275. eng, Q., Cieszewski, C., & Madden, M. (2009). Large area forest inventory using landsat etm+: A geostatistical approach. ISPRS Journal of Photogrammetry and Remote Sensing , 64, 27 – 36.Mutanga, O., & Rugege, D. (2006). Integrating remote sensing and spatial statistics to model herbaceous biomass distributionin a tropical savanna.
International Journal of Remote Sensing , 27, 3499–3514.Næsset, E. (2004). Practical large-scale forest stand inventory using a small-footprint airborne scanning laser.
ScandinavianJournal of Forest Research , 19, 164–179.Næsset, E. (2011). Estimating above-ground biomass in young forests with airborne laser scanning.
International Journal ofRemote Sensing , 32, 473–501.Neigh, C., Nelson, R., Ranson, K., Margolis, H., Montesano, P., Sun, G., Kharuk, V., Næsset, E., Wulder, M., & Andersen,H. (2013). Taking stock of circumboreal forest carbon with ground measurements, airborne and spaceborne lidar.
RemoteSensing of Environment , 137, 274–287.Nelson, R. (2010). Model effects on GLAS-based regional estimates of forest biomass and carbon.
International Journal ofRemote Sensing , 31, 1359–1372.Nelson, R., Gobakken, T., Næsset, E., Gregoire, T., St˚ahl, G., Holm, S., & Flewelling, J. (2012). Lidar sampling – using anairborne profiler to estimate forest biomass in Hedmark County, Norway.
Remote Sensing of Environment , 123, 563–578.Nelson, R., Short, A., & Valenti, M. (2004). Measuring biomass and carbon in Delaware using an airborne profiling lidar.
Scandinavian Journal of Forest Research , 19, 500–511.Pan, Y., Birdsey, R.A., Fang, J., Houghton, R., Kauppi, P.E., Kurz, W.A., Phillips, O.L., Shvidenko, A., Lewis, S.L., Canadell,J.G., Ciais, P., Jackson, R.B., Pacala, S.W., McGuire, A.D., Piao, S., Rautiainen, A., Sitch, S., & Hayes, D. (2011). A largeand persistent carbon sink in the world’s forests.
Science , 333, 988–993.Pattison, R., Andersen, H.E., Gray, A.N., Schulz, B.K., Smith, K., Jovan, S., & Manies, K.L. (In prep). Forests of the TananaValley State Forest and Tetlin National Wildlife Refuge Alaska - results of the 2014 pilot inventory. US Department ofAgriculture Forest Service, Pacific Northwest Research Station.Pfeffermann, D. (2013). New important developments in small area estimation.
Statist. Sci. , 28, 40–68.Pflugmacher, D., Cohen, W.B., Kennedy, R.E., & Yang, Z. (2014). Using Landsat-derived disturbance and recovery historyand lidar to map forest biomass dynamics.
Remote Sensing of Environment , 151, 124–137.Powell, S.L., Cohen, W.B., Healey, S.P., Kennedy, R.E., Moisen, G.G., Pierce, K.B., & Ohmann, J.L. (2010). Quantificationof live aboveground forest biomass dynamics with Landsat time-series and field inventory data: A comparison of empiricalmodeling approaches.
Remote Sensing of Environment , 114, 1053–1068.Rieger, S., Furbush, C.E., Schoephorster, D.B., Summerfield Jr, H., & Geiger, L.C. (1972). Soils of the Caribou-Poker CreeksResearch Watershed, interior Alaska. Technical Report. DTIC Document.Rinehart, A.J., Jones Jr, J.B., & Harms, T.K. (2015). Hydrologic and biogeochemical influences on carbon processing in theriparian zone of a subarctic stream.
Freshwater Science , 34, 222–232.Saarela, S., Grafstr¨om, A., St˚ahl, G., Kangas, A., Holopainen, M., Tuominen, S., Nordkvist, K., & Hyypp¨a, J. (2015). Model-assisted estimation of growing stock volume using different combinations of lidar and Landsat data as auxiliary information.
Remote Sensing of Environment , 158, 431–440.S¨arndal, C., Swensson, B., & Wretman, J. (1992). Model Assisted Survey Sampling. Springer series in statistics, Springer-Verlag.Schelin, L., & Sj¨ostedt-De Luna, S. (2010). Kriging prediction intervals based on semiparametric bootstrap.
MathematicalGeosciences , 42, 985–1000.Sj¨ostedt-De Luna, S., & Young, A. (2003). The bootstrap and kriging prediction intervals.
Scandinavian Journal of Statistics ,30, 175–192.Stow, C.A., Reckhow, K.H., & Qian, S.S. (2006). A bayesian approach to retransformation bias in transformed regression.
Ecology , 87, 1472–1477.Tanaka-Oda, A., Kenzo, T., Toriyama, J., & Matsuura, Y. (2016). Variability in the growth rates and foliage δ
15 N values ofblack spruce trees across a slope gradient in the Alaskan interior.
Canadian Journal of Forest Research , .Thompson, S. (2002). Sampling. Wiley Series in Probability and Statistics, Wiley.Tsui, O.W., Coops, N.C., Wulder, M.A., & Marshall, P.L. (2013). Integrating airborne lidar and space-borne radar viamultivariate kriging to estimate above-ground biomass.
Remote Sensing of Environment , 139, 340–352.UN-REDD (2009). The UN-REDD programme. . Accessed: 8-11-2016.U.S. Fish and Wildlife Service (2016). Tetlin National Wildlife Refuge, Alaska. . Accessed: 10-16-2016.Ver Hoef, J. (2002). Sampling and geostatistics for spatial data.
Ecoscience , 9, 152–161.Ver Hoef, J., Cressie, N., & Barry, R.P. (2004). Flexible spatial models for kriging and cokriging using moving averages andthe fast fourier transform (FFT).
Journal of Computational and Graphical Statistics , 13, 265–282.Verbyla, D. (2011). Statewide 2-km raster of year since last widlfire, Bonanza Creek LTER - University of Alaska Fairbanks. .Wang, G., Oyana, T., Zhang, M., Adu-Prah, S., Zeng, S., Lin, H., & Se, J. (2009). Mapping and spatial uncertainty analysis offorest vegetation carbon by combining national forest inventory data and satellite images.
Forest Ecology and Management ,258, 1275 – 1283.White, J.C., Wulder, M.A., Varhola, A., Vastaranta, M., Coops Nicholas, C., Cook, B.D., Pitt, D., & Woods, M. (2013). A bestpractices guide for generating forest inventory attributes from airborne laser scanning data using an area-based approach.
The Forestry Chronicle , 89, 722–723.Woodall, C.W., Coulston, J.W., Domke, G.M., Walters, B.F., Wear, D.N., Smith, J.E., Andersen, H.E., Clough, B.J., Cohen,W.B., Griffith, D.M., et al. (2015). The US forest carbon accounting framework: Stocks and stock change, 1990-2016, .Wulder, M.A., White, J.C., Nelson, R.F., Næsset, E., Ørka, H.O., Coops, N.C., Hilker, T., Bater, C.W., & Gobakken, T. Remote Sensing of Environment , 121, 196–209.Yava¸slı, D.D. (2016). Estimation of above ground forest biomass at mu˘gla using icesat/glas and landsat data.
Remote SensingApplications: Society and Environment , 4, 211 – 218.Zheng, D., Rademacher, J., Chen, J., Crow, T., Bresee, M., Moine, J.L., & Ryu, S.R. (2004). Estimating aboveground biomassusing Landsat 7 ETM+ data across a managed landscape in northern Wisconsin, USA.
Remote Sensing of Environment ,93, 402–411.,93, 402–411.