L-Band SAR Backscatter Related to Forest Cover, Height and Aboveground Biomass at Multiple Spatial Scales across Denmark

: Mapping forest aboveground biomass (AGB) using satellite data is an important task, particularly for reporting of carbon stocks and changes under climate change legislation. It is known that AGB can be mapped using synthetic aperture radar (SAR), but relationships between AGB and radar backscatter may be confounded by variations in biophysical forest structure (density, height or cover fraction) and differences in the resolution of satellite and ground data. Here, we attempt to quantify the effect of these factors by relating L-band ALOS PALSAR HV backscatter and unique country-wide LiDAR-derived maps of vegetation penetrability, height and AGB over strong linear relation ( R 2 = 0.79 at 250-m scale). In areas of high fractional forest cover, there is a slight decline in backscatter as AGB increases, indicating signal attenuation. The two results demonstrate that accounting for spatial scale and variations in forest structure, such as cover fraction, will greatly beneﬁt establishing adequate plot-sizes for SAR calibration and the accuracy of derived AGB maps.


Introduction
Spatially explicit measures of aboveground forest biomass (AGB) are important for understanding the terrestrial carbon cycle and have been largely addressed in the context of global environmental change [1,2].Maps of AGB may be used as tools to monitor forests by capturing regrowth, deforestation and degradation processes and can support efforts towards conservation and enhancement of forest carbon stock, for example under the initiative of Reducing Emissions from Deforestation and Forest Degradation (REDD+) [3].Some studies have suggested that AGB must be mapped at 1-ha to 100-ha resolution with an accuracy of 20% or ±20 Mg ha −1 , the greater of the two, up to a maximum of ±50 Mg ha −1 for policy implementation [4,5].However, existing national, continental and global maps produced using Earth observation data still face challenges in presenting the spatial distribution of biomass with high accuracy [6][7][8].This may constrain evaluating the effectiveness of policies targeted towards forest carbon conservation.
Large-scale mapping of aboveground carbon stocks has benefited from the wide-coverage remote sensing satellite data, calibrated with ground measurements of AGB based on field inventories, e.g., [9][10][11].One of the sources of uncertainty in mapping AGB, however, concerns the incompatibility of the resolution (pixel size) of the satellite imagery and the size of forest plots measured on the ground [4,12,13].Forest inventories are expensive, time consuming and labor intensive and are thus geographically limited.Comparatively small field plots may be associated with larger sampling errors due to local spatial variability of AGB, as well as geo-location mismatches with satellite imagery [14,15].
Of the available remote sensing technologies, airborne light detection and ranging (LiDAR) is a powerful tool for measuring forest stand structure and extrapolating measurements spatially [16][17][18][19].Vegetation LiDAR instruments typically generate laser pulses at wavelengths between 900 and 1064 nm and record the time interval between emission and receipt to generate the sensor-object range (distance).Small-footprint (0.25 to 2 m in diameter) discrete-return systems digitize a combination of first, last and intermediate pulse returns, which can directly describe important forest characteristics, such as height distributions and gap fractions [19].These are subsequently related to AGB via correlative models calibrated with ground measurements [20].LiDAR techniques have hence been discussed in the context of use for REDD+ as a monitoring tool [1,16,21,22].However, in complex and heterogeneous forests, the use of LiDAR may be hindered by cloud cover, edge effects or overhanging neighboring tree canopy [13,18].Further, airborne LiDAR surveys are expensive, making them a relatively inefficient tool to map forest biomass periodically over large areas (although, see Mascaro et al. [22]).Spaceborne LiDAR offers an alternative, though it cannot achieve quasi-full coverage, like airborne systems, due to beam dispersal and power constraints, restricting it to large footprints separated by hundreds of meters.Although there are no current ongoing missions, the future ICESat-2, GEDIand MOLIproducts will be used for vegetation mapping [23].
There has also been growing interest in the use of microwave synthetic aperture radar (SAR) to estimate AGB.SAR offers some advantages over LiDAR, allowing continuous coverage over vast areas and consistent and frequent acquisitions achievable from spaceborne platforms [15].Microwave pulses are transmitted to the Earth's surface at shallow incidence angles (compared to airborne LiDAR scanning), and the backscattered signal received in return can be used to interpret information of land surface structure [24].Backscatter is characterized by polarized modes of the transmitted and received signals (e.g., horizontal send and horizontal receive, HH, and horizontal send and vertical receive, HV), sensitive to biophysical land parameters.Scattering mechanisms that contribute to backscattered energy include single bounce from rough surfaces and water, double bounce from edges (ground and trunks in forests) and volumetric scattering, such as that from inhomogeneous forest canopy and stems.Unlike LiDAR, long microwave wavelengths are of the same order as the sizes of elements of vegetation that contain the vast majority of aboveground biomass (stems and branches).Fundamentally, the larger the diameter of these elements, the greater the backscatter.Long wavelengths can penetrate further into the canopy, but deep and dense canopies, typical in high biomass forests, may cause the volumetric scattered signal to weaken (signal attenuation) [24][25][26].Differences in environmental variables, such as moisture, flooding, topography, etc., also cause variations in backscatter unrelated to the structure of vegetation elements.
Backscatter from long-wavelength SAR imaging (e.g., ∼23-cm L-band) has been used to retrieve AGB density in numerous previous studies [26][27][28][29][30][31].The research has also shown that SAR signals lose sensitivity to AGB, or "saturate", beyond a certain AGB value for a given wavelength, incident angle and forest type [32][33][34][35].AGB retrieval models have commonly used a range of ground-measured plot sizes (e.g., 0.01 ha [36] to >1 ha [31]), a variety of physical relations (e.g., simplified water cloud model [37,38], log relation [39,40] or power-law [11]) and have concluded a wide range of AGB saturation windows (e.g., 40 Mg ha −1 to 150 Mg ha −1 for L-band SAR [41]).Results of these studies are of immense importance to our understanding of quantifying and mapping AGB using SAR.However, it is well known that backscatter is not a direct measure of AGB density, and variability in vegetation structure may have a large impact on backscatter intensity [42].Studies modeling the interaction of microwaves with forests have demonstrated that accounting for stem distribution, size and number provides a stronger description of the AGB-backscatter relation [33,[43][44][45].The absence of extensive and large-scale independent field datasets that estimate vegetation structural properties, however, challenge characterizing the relation.Further development towards relating backscatter to forest structural properties, such as those readily measurable using airborne LiDAR, is hence still required.Such relations will aid producing coarse-resolution (compared to field-survey plot sizes) and regional-to national-scale AGB maps, with higher confidence in their associated uncertainties.
The scope of this paper is to investigate the relations between L-band SAR backscatter and measures of vertical (height) and horizontal (cover fraction) forest structure and AGB at multiple spatial scales across Denmark.The latter variables are derived using a unique dataset of small-footprint, country-wide, airborne LiDAR scanning and ground-measured AGB available for the country.Our study aims to answer the following questions: • How does L-band SAR backscatter relate comparatively to forest height, cover and AGB at fine (50 m × 50 m) to coarse resolutions (500 m × 500 m)?
• To what degree is the variation in the AGB-backscatter relation explained by the fraction of forest cover or canopy density?• What is the optimal scale to retrieve and map AGB using SAR in Denmark?

Study Area
Denmark lies in the central European temperate forest zone and contains broadleaf and conifer species types.Forested, wooded and protected forest areas covered 608,078 ha (20,355 ha unstocked) or 14.1% of the total land area in 2012 [46], while intensive agriculture dominates the rest.The amount of crop biomass on fields during late fall/spring, corresponding to the period of this study, is very low or absent [47].The topography of the country is flat, with <1% area with slopes >10 • at 100 m × 100 m map resolution.
Even-aged (planted and natural regeneration) forests cover 75% of the forested area; uneven-aged natural forests cover 7%; and the remaining forested areas are of mixed types, including operational uneven-aged forests, protected forests or those under other management forms [46].The forests are highly fragmented spatially, as a result of active government measures to increase the country's proportion of forested areas [47].Most forest patches are <10 ha in size, and only a few larger forests cover areas of >100 ha; these are also commonly split into smaller patches based on their history and establishment.Further, hedgerows are an important component of vegetated areas, covering 59,659 ha of land in 2012 and being accounted for as carbon sinks in the Danish emissions inventory [48].Given the complexity of vegetation distribution in the landscape, we have not distinguished between large forests, small plantations and hedgerows for this study.
Danish forests are characterized by a high degree of active management, with assisted regeneration, thinning and pre-commercial thinning.This results in forests with a large variation in structural characteristics (stem number, mean diameter and heights) and, at the same time, constrains the natural variation of biomass ranges.On the 2130 forested plots measured in the Danish National Forest Inventory (NFI) in 2006/2007, average AGB was 100 Mg ha −1 , and only 17% of plots had AGB >300 Mg ha −1 [49].

Species Trial Plots
Forest tree species trial plots (STPs) were established in Denmark in the 1960s, with each containing even-aged single-species trees distributed in clusters across 13 sites (Figure 1a).Although not set up for this purpose, the design of the STPs made them well suited for investigating relations between structural forest parameters and high-resolution airborne LiDAR data.Sections of STPs were selected to reduce the influence of neighboring crowns, by placing plots within the centers of STP stands, with edges bisecting the same species and relatively similar height trees (Figure 2).The plots were inventoried in 2007/2008, with all trees cross-calipered at breast height and a subsample of 30 randomly selected trees measured for height.Biomass was estimated using wood density and expansion factors described in Skovsgaard et al. [50] and Skovsgaard and Nord-Larsen [51] in the study by Nord-Larsen and Schumacher [49].Plot-level AGB estimates (AGB STP ), including aboveground stump, stem and branches (without foliage), were provided for this study.Plots that were lost during stand establishment or due to catastrophic wind throws were removed for this study, resulting in a dataset of 113 plots of an area of 0.07 ha to 0.23 ha (Table 1).The species represented in them are the most commonly grown deciduous and conifers in Denmark [46], covering ∼68% of the forest area in Denmark (conifers cover 86% of coniferous areas and broadleaves cover 49% of broadleaf areas).AGB STP ranged 47-629 Mg ha −1 , similar to the range measured in the NFI across the country in 2007/2008.

Airborne LiDAR Scans
Wall-to-wall discrete-return airborne LiDAR data with an average point density of 0.5 pulses/m 2 (Table 2) were collected over Denmark during leaf-off conditions in spring 2006 and fall 2006/spring 2007 by COWI.First and last return pulses were recorded and extrapolated to generate a digital elevation model (DEM) and to extract the elevation of structures above ground by the data provider.The elevation dataset has previously been shown to predict dominant tree height with a precision of ±2.5 m [52].For this study, uncalibrated flight strips were identified upon examining the dataset and removed (Figure 1a).

ALOS PALSAR Scenes
Eleven ALOS PALSAR scenes in fine-beam dual-polarization mode (horizontal send and receive (HH) and horizontal send and vertical receive (HV)), were acquired at Level 1.5 (see Section 2.4.1) in 2007 (Table 3).The scenes each have an approximate dimension of 59 km in azimuth and 70 km in range, together covering ∼78% of the extent of Denmark (Figure 1a).In selecting scenes, regional weather conditions for every day in the week prior to image acquisition were obtained from the Danish Meteorological Institute.Since Denmark has a typically wet climate, scenes with least precipitation around approximately the time of the airborne LiDAR scanning were sought.Pre-processing visual comparisons in areas of scene overlap, particularly over homogeneous land cover types, showed no distinguishable differences in backscatter.2.4.1.Pre-Processing Level 1.5 data are pre-processed for radiometric and geometric corrections after range and multi-look azimuth compression (i.e., 4 looks in azimuth and 1 in range, to a 12.5 m × 12.5 m pixel size).Terrain correcting and geo-referencing the images using a 10 m × 10 m resolution DEM (derived from a different LiDAR dataset over Denmark and available from Geodatastyrelsen) resulted in a reduction in terrain-induced distortion and geolocation errors.Georeferencing was verified by comparison to wall-to-wall, 0.16 m × 0.16 m resolution, aerial photography (available at COWI).Features, such as hedges of farms, small forest patches and islands, were found to be accurately geolocated to the sub-pixel level on the corrected images.The processing was done using the software, ASF MapReady 3.1.22(Alaska Satellite Facility, USA).
Normalized backscatter for both polarizations (σ 0 HH and σ 0 HV ) in the power domain (m 2 /m 2 ) was mapped at a 12.5-m pixel size.Preliminary analysis showed that σ 0 HV had a stronger relation to AGB than σ 0 HH (Supplementary Material, Figure S1), similar to the results of previous studies covering a variety of forest ecosystems, e.g., [31,[38][39][40]53,54].As vegetation tends to depolarize the radar signal [24], it was relevant to focus this study on AGB, vegetation cover and σ 0 HV relations.

Speckle Filtering
All SAR images have random speckle, which may mask the relationship between received backscatter signals and relevant surface interactions.Speckle must therefore be reduced before quantitative analysis, either by speckle filtering or extensive multi-looking, the latter dramatically reducing image resolution.To avoid coarsening resolution (since relations of SAR and biophysical vegetation parameters at multiple resolutions are investigated in this study), different spatial speckle reduction methods (e.g., median, Lee-sigma, local-region and Frost) and window sizes (e.g., 3 × 3, 5 × 5, 7 × 7 pixels) were tested.The aim was forest edge sharpness preservation, texture preservation and computational efficiency over the large study area.A comparison of filtered images over plantations showed that the enhanced Lee filter [55] at a 5 × 5 pixel window was most suited, as it retained stand edges, yet significantly reduced the standard deviation of σ 0 HV over individual plantations (Supplementary Material, Figure S2).Upon coarsening map resolutions (Section 2.5.2 and Section 2.5.3),there was a further reduction in speckle, particularly at hectare scales and larger.

LiDAR-to-AGB STP Model
Point cloud elevation data from each STP were used to build a LiDAR-to-AGB STP model.We derived variables, including the mean, maximum and minimum pulse height above ground; relative mean first pulse height (relative to maximum pulse height); variance, standard deviation, skewness, kurtosis and coefficient of variation of heights; interception ratio by canopy; and height distribution percentiles (ranging from the 25th, 50th, 75th, 95th to 99th percentile) for each STP.Most LiDAR metrics describing canopy elevation (height percentiles; maximum height; mean height; and coefficient variation of height) are highly correlated and have a similar correlation structure as vegetation biomass [56][57][58].To avoid model multicollinearity, canopy height metrics showing the highest correlation with AGB STP were selected based on the correlation matrix and cluster analysis conducted by Nord-Larsen and Schumacher [49].However, the uncalibrated return pulse intensity used in their study showed flight-path stripes when mapped.The use of the variable would have introduced spatial bias in mapping AGB in our study, where pixel-level accuracy is required for comparing maps, and was thus avoided.
Previous studies have demonstrated that structural parameters that combine height and gap fraction improve AGB estimates from LiDAR [57,59].A variable measuring the interception of LiDAR signals by vegetation, the vegetation interception ratio (VIR), was derived in our study.Several definitions of VIR were tested to find one that most captured differences within plantations of differing tree density and height, including the fraction of returns from various heights within plantations and the fraction of returns from the ground (Supplementary Material, Figure S3).It was found that VIR was best defined as the ratio of the number of first returns from above a 1-m height (FR >1m ) to the total number of first returns (FR) (Equation ( 1)).
The best LiDAR-to-AGB STP model was selected based on the lowest Akaike information criterion (AIC) and normality of residuals.Due to the small sample size of AGB STP , robust accuracy estimates were assessed using jackknife resampling, leaving out a random sample of 30% of data for testing in 5000 iterations.The jackknife mean relative absolute error (mRAE), mean bias error (mBias), mean root mean squared error (mRMSE) and the model coefficient of determination (R 2 ) were used to quantify uncertainty.Dependency on plot sizes was tested by examining fit residuals after model development.
For comparison, we also attempted relating AGB STP to σ 0 HV using 55 STPs contained within the PALSAR extent.However, the plot areas were small compared to SAR image resolution, such that the effects of speckle were not visibly minimized on each plot, and a large number of AGB STP values lay well above the range of SAR signal saturation.This made it impossible to map AGB with reasonable accuracy using only SAR data (Supplementary Material, Figure S4), requiring us to use the LiDAR-to-AGB STP model to derive a country-wide AGB map.2.5.2.SAR-to-VIR L and SAR-to-MH L Models Variable VIR was first derived over the fully-forested STPs to capture vegetation penetrability (Equation ( 1)).On pixels with sparse vegetation or low canopy cover, the variable summarizes the heterogeneity of vegetation.Some caution must be taken in the interpretation of VIR.VIR may be positively related to, but is not, however, a direct measure of the vertical component of canopy density, since signals may be intercepted at any height above the ground >1 m.As a result, a (hypothetical) lower height, but uniformly-vegetated pixel may have the same VIR as a higher height and uniformly-vegetated pixel.Further, a lower height, but uniformly-vegetated pixel may have the same AGB density, but higher average VIR than a pixel that has only 50% tree cover, but tall and large trees (Figure 3).This variation is further complicated by differences in species types, which may have similar AGB density at a pixel scale, but varying tree distribution, stem and branch number and canopy structure.For simplicity, this study interprets VIR only as a measure of the proportion of vegetation cover (as seen by LiDAR) in each mapped pixel, independent of mean canopy height or canopy density.This portion may also be low for stands containing very young trees or low stem density.The implications of height and canopy density being related to VIR are discussed in Section 4.3.VIR was first mapped at a 50 m × 50 m pixel size (50-m "scale") (Figure 1b), and pixels were then averaged to coarser scales, ranging from 100 m to 500 m.SAR backscatter, σ 0 HV , at the 12.5-m scale was averaged in the power domain (m 2 /m 2 ) to produce coarser scale images and regressed against average VIR both in the power and dB domain (SAR-to-VIR L model).For comparison, σ 0 HV was also regressed against mean canopy height (MH; SAR-to-MH L model).At each scale, large built-up and water areas were removed entirely using a high-resolution (10 m × 10 m pixel size) land use map (available at Geodatastyrelsen).Pixels with infrastructure surrounded by agricultural lands and forested areas were removed if the infrastructure occupied ≥25% of the pixel area (e.g., Figure 1d).We verified that this threshold removed most problematic areas from the analysis and did not impact our results.
For each model, the large number of pixels for calibration/validation ensures robustness without the need for jackknifing.Therefore, the dataset was randomly split into 50% training and 50% testing subsamples, which ensures confidence in the results, as the training data are independent of the testing data.All regressions were run using the CURVEFIT procedure in the software EXELIS IDL (Exelis Visual Information Solutions, Boulder, CO, USA).

SAR-to-AGB L Model
The LiDAR-to-AGB STP model was used to map AGB density across Denmark (AGB L in units of Mg ha −1 ) at 50-m scale first, comparable to the largest STP plot (0.23 ha), and AGB L pixels were subsequently aggregated by averaging to coarser scales.Built-up and water areas were removed using the same procedure as in Section 2.5.2, before testing the σ 0 HV to AGB L relation (SAR-to-AGB L model).There is a large volume of studies, e.g., [11,31,37,39,53,62,63], that have developed AGB-backscatter relations through either radiative transfer models or empirical regressions [64].The familiar loss of sensitivity of σ 0 HV to high AGB values (typically ∼100 Mg ha −1 for L-band SAR) is apparent in most of these studies.Further, backscatter is often interpreted as a function of biomass, and σ 0 HV is thus usually used as a response variable.While it is interesting to analyze the effect of using different models on our data, it is not relevant to the aims of this study.Instead, we focus on finding a single, statistically robust fit that can describe the backscatter-biomass relation best for the datasets.The fit with a balance between the least bias, least root mean squared error (RMSE) and well-distributed residuals was chosen as the best-fit model, using 50% of the dataset for training and 50% for testing.A sensitivity analysis to errors in AGB L (mRMSE and mBias obtained through the jackknife resampling, Section 2.5.1) was conducted using simulations discussed in Section 3.2.4.We assumed throughout that the model error is additive, i.e., the random error on the model fit does not increase as variable values increase.
Although we focus on σ 0 HV relations to vegetation cover in this study, we tested whether models, including both σ 0 HV and σ 0 HH , e.g., [11,65], provided lower RMSE during AGB L prediction than those using σ 0 HV alone.We further tested whether σ 0 HH had any relation to residuals of our selected model, to assess its ability to improve AGB L prediction.

Predicting AGB ST P with LiDAR
There was a strong relation between most LiDAR metrics and AGB STP , modeled using exponential or power functions.The regression with the lowest AIC was an exponential function, including two predictor variables: the 95th percentile of canopy height from returns above 1 m (P 95 ), multiplied by VIR, explaining 70% of the variation in AGB STP (Figure 4a, Table 4).

Multi-Scale SAR Relations to Vegetation Cover, Height and AGB L
The following sections detail relations between SAR backscatter and LiDAR-derived vegetation descriptors for the whole study area.Results from only a few spatial scales are reported to demonstrate the observed trends.

Up-Scaling and Forest Edges
In aggregating pixels from fine to coarse scales, consideration must be given to forest edges, which are lost by averaging over them.The inclusion or exclusion of edge pixels was found to influence model performance and retrieval accuracy.Thus, a number of methods were tested to fairly represent only forested areas and pixels along forest edges, as scales are coarsened from 50 m on, before any quantitative analysis.First, areas with AGB ≤ 5 Mg ha −1 were flagged and removed from the maps at 50-m scale.Further, some agricultural lands were identified with high VIR (>0.7) and low mean and maximum height (<1.5 m) and removed.These thresholds allowed a large portion of bare-ground or agricultural fields to be removed, i.e., ∼89% of pixels classified as 'potentially agriculture' in the Danish land use map (Geodatastyrelsen) discarded as a result.The remaining pixels identified as 'potentially agriculture' were found to be: (1) vegetated hedgerows; (2) forest-edges; and (3) areas of young forest plantations, divisions between individual plantations and, in some cases, harvested forests.These were retained for analysis at 50-m scale.For each coarser scale map, pixels were only included if they had AGB >5 Mg ha −1 and if less than a certain threshold of discarded pixels at 50 m were contained with them (Supplementary Material, Figure S5).This threshold was selected such that the average AGB represented across the landscape at each spatial scale remains nearly constant (Figure 5).Hence, primarily vegetated pixels are being examined at each scale and nearly the same forest stands being compared across scales, with the loss of comparatively small forest patches as scales coarsen.Once this was done, there was no significant influence of remaining forest-edge pixels on the general trends in SAR relations to LiDAR variables when comparing results across scales.

SAR-to-VIR L Model
Backscatter in the power domain (m 2 /m 2 ) was found to have a positive linear relation to average VIR (Figure 6, left).
The scatter in the data is relatively weakly explained at 50-m scale (R 2 = 0.63) and improves at coarser scales (up to R 2 = 0.79 at ≥ 250 m) (Table 6).
As there is no apparent loss of sensitivity of σ 0 HV to high VIR values, it was reasonable to build the model in the power domain, making it possible to distinguish high σ 0 HV at intervals comparable to low σ 0 HV values.However, since this trend is not expected in the SAR-to-MH L and SAR-to-AGB L models (due to signal saturation), we also tested and report SAR-to-VIR L performance in the dB domain (Figure 6, right) using the best-fitted exponential model (Y = α 0 e α 1 X + α 2 ).Trends across spatial scales were found to be the same, with prediction RMSE improving with coarsening scale (Table 6, model coefficients not shown).

SAR-to-AGB L and SAR-to-MH L Model
The σ 0 HV images showed good qualitative agreement with the AGB L map at 50-m scale, with no observable geolocation errors.The variation in σ 0 HV across the landscape was best explained (highest R 2 ) by an SAR-to-AGB L model defined as: where a, b and c are constants.Although the fit is empirical, it resembles the simplified form of the water cloud model [66] used in previous studies [37,38,53,67] to relate AGB and backscatter.Here, it was chosen due to the high R 2 values and low RMSE values achieved (Table 7), rather than for theoretical reasons.It must be noted that the fitted model accommodates for a decrease in backscatter with increasing AGB, which may physically be caused by increasing forest opacity (as explained by the water cloud model), as well as macroecological trends in stem sizes and number density as AGB increases [68].The strongest relations across all scales resulted from fitting Equation ( 4) both in dB (Figure 7) and the power domain (Supplementary Material, Figure S6).The scatter in the AGB L -backscatter relation was best explained at 100-150-m scales (R 2 = 0.62) and weakened in comparison at scales ≥250 m (R 2 ≤ 0.58), while the RMSE improves with coarsening resolution.This difference arises since the RMSE is more sensitive to larger and extreme values due to the squaring of differences, which are averaged out with coarsening resolution, as compared to R 2 , which represents how well the modeled curve explains the scatter in the data.The SAR-to-MH L model was also best developed using Equation (4).A loss of sensitivity of σ 0 HV to canopy heights of ∼7 m at 50 m (Figure 8) and a relatively weaker relation of backscatter to MH than backscatter to VIR (Supplementary Material, Figure S7, shows the SAR-to-MH L in the power domain) were observed.The scatter in the data is comparatively better explained at fine scales (R 2 = 0.62 at 100 m), with some decline in the robustness of the fit at scales ≥250 m (Table 8).Similar to the SAR-to-AGB L model, the prediction RMSE improves with coarsening resolution.4) to be inverted, resulting in a logarithmic function.
Since the inverted model is undefined for values where σ 0 HV ≥ c + a and predicts AGB < 0 Mg ha −1 where σ 0 HV ≤ c, the inversion requires the model to be constrained to a certain range of σ 0 HV values and a best-guess to be made where σ 0 HV is below/above the predictable range [38].Inverting the model to retrieve AGB results in a loss of 23%, 20%, 18% and 17% of the data where σ 0 HV fell above the predictable range and 8%, 6%, 4% and 3% where it fell below, at scales of 50 m, 100 m, 250 m and 500 m, respectively (the number of observations, N, are reported in Figure 7).For where it fell below the predictable range, 0 Mg ha −1 was considered the most reasonable best guess and assigned to predicted values.
In assessing overall retrieval accuracy, two approaches were taken.First, AGB prediction RMSE and bias were derived only for values where σ 0 HV fell below and within the predictable range.Second, a maximum retrieval value (MRV) was assigned to values where σ 0 HV fell above the predictable range.In this method, an approach similar to that described in Santoro et al. [63] was taken and a buffer zone set to remove outliers that were most likely resultant from local environmental factors.We then tested assigning the maximum successfully predicted AGB value and increasing/decreasing the value each time by 5 Mg ha −1 in multiple iterations until the minimum overall model RMSE was achieved.Furthermore, the percentile distribution of successfully predicted values was derived, and each time, a different percentile was set as the MRV until the minimum model RMSE was achieved.It should be noted that no consideration was given to bias in this method, since the overall model bias is strongly sensitive to the selected MRV, and low bias does not necessarily indicate that the best method has been chosen.Multiple test iterations showed that assigning an MRV of one of the percentile distributions of the successfully predicted data reduced the overall model RMSE most.The effect of assigning an MRV to the residual plots of predicted AGB against AGB L is displayed in Figure 7, right.
The differences in the two approaches for multi-scale comparisons showed that the method of assigning an MRV changes absolute and relative model errors, but the trends of improving model accuracy with spatial scale remain largely similar (Figure 9).AGB retrieval accuracy and residual bias were found to improve overall with coarsening map scale.The largest improvement is observed on mapping at 100-m instead of 50-m scale, and after 250-m scale, accuracy improvement is low and even decreases with increasing scale by some metrics thereafter.HV as a function of spatial scale before and after applying a maximum retrieval value (MRV) to predicted AGB.
It was also found that the chosen model (Equation ( 4)) provided lower RMSE values in AGB prediction overall, when compared to models including both σ 0 HV and σ 0 HH .Further, σ 0 HH showed a negligibly weak relation to AGB residuals, suggesting that it adds no further information to AGB L prediction than that which σ 0 HV has already provided (Supplementary Material, Figure S8).In accounting for forest edges in the multi-scale analysis (Section 3.2.1), the mean AGB of the pixels considered was nearly constant, but the distribution of measures changes as spatial variability decreases upon averaging (Figure 5).To fairly assess prediction trends across scales, we derived relative prediction bias and RMSE at each range of the independent AGB variable, which shows the degree to which the inverted model differs in accuracy for areas of different biomass ranges (Figure 10).In this analysis, we use the predicted AGB dataset before setting an MRV for the results data, to remove the stronger influence of the chosen MRV at higher AGB L ranges compared to lower AGB L ranges.A continuous moving bin of 50 Mg ha −1 was examined individually across the range of AGB L .It was found that high biomass values are under-predicted by the model, particularly at fine scales, indicating a decrease of sensitivity of σ 0 HV to AGB at an AGB value that varied with the mapping scale (Figure 10a,b).The model begins to under-predict AGB in areas with biomass >45 Mg ha −1 and >80 Mg ha −1 at 50-m and 500-m scale, respectively.RMSE in the bins indicated a stronger performance at 100-m scale for areas with AGB below 110 Mg ha −1 and some improvement of performance when mapping at scales ≥150 m for areas with higher AGB (Figure 10c,d).
Since AGB L was derived using VIR and P 95 at 50-m scale, an analysis of residual plots of the SAR-to-AGB L model was conducted to assess whether these variables offered an explanation for the scatter in the biomass-backscatter plot.VIR showed a weak (R 2 = 0.11), but significant (p < 0.0001) linear relation to residual σ 0 HV at coarse scales (≥250 m), while MH and height distribution percentiles were found to provide no additional statistical information (Figure 11a,b).It appears, therefore, that the AGB-backscatter curve is largely explained by the MH-backscatter relation, and VIR provides some additional information (i.e., 11% at 500 m) on the data scatter beyond the modeled AGB-backscatter curve.It is important to note that the analysis cannot assess whether this information can improve AGB retrieval accuracy, since AGB L is not independent of VIR in this study.As a means of demonstration of how the fraction of forest cover can explain variations in the AGB-backscatter curve, we stratified our study area into four VIR categories (0 to 0.25, 0.25 to 0.50, 0.50 to 0.75 and 0.75 to 1.0) and re-built the SAR-to-AGB L model on each.We allowed for the selection of both linear and asymptotic functions (e.g., Equation ( 4)) for the model based on R 2 values.We found significant differences in the shape of the biomass-backscatter curve for different VIR ranges (Figure 12).Areas of low VIR (primarily low stem density plantations and areas of low fraction of vegetation cover at VIR ≤ 0.5) were well fitted using Equation (4).In contrast, areas of high VIR (primarily fully-vegetated core plantation pixels at VIR > 0.75) showed a slight decline in σ 0 HV as AGB L increases.All model coefficients were significant (p < 0.0001), and σ 0 HV is shown in the power domain (m 2 /m 2 ) in Figure 12 to demonstrate the decline at high σ 0 HV values.

Uncertainty Analysis in SAR-to-AGB L
The methods of this study are primarily empirical, involving large datasets, each with inherent uncertainty.It must be specified here that the purpose of our study is to examine trends in relations and retrieval uncertainty across scales, rather than to achieve the highest possible prediction accuracy at any one scale; we do not attempt to estimate the total AGB stock of Denmark.As such, we do not integrate datasets, do not restrict our dataset to specific AGB ranges (e.g., AGB L < 100 Mg ha −1 , where σ 0 HV signal saturation may not be problematic) and assess the models' performances by using a subset of our data for testing.
We distinguish between uncertainty related to bias and measurement errors in the original datasets.The former are only likely to cause a systematic shift in AGB L values consistently across scales and, hence, are not relevant to multi-scale comparisons.This includes allometric-related bias in AGB STP and systematic errors in the LiDAR data calibration.In estimating AGB STP , all trees were cross-calipered, and hence, the sampling errors are negligible.Measurement errors in LiDAR returns used to produce AGB L are considered negligible given the high point density of returns (Table 2) pooled together at 50-m scale.
Errors in the AGB L map derived from the LiDAR-to-AGB STP model (i.e., mRMSE and mBias) result in higher uncertainty in AGB retrieved from the SAR-to-AGB L model.Two complications arise in propagating errors from one map to the other.First, the accuracy of AGB predicted by LiDAR has been previously shown to increase with coarsening resolution [19,69], but we lacked ground data to build multi-scale LiDAR-to-AGB STP models.Since the best possible map was produced at 50-m scale and then averaged to coarser scales, we take a conservative approach and assume no reduction in AGB prediction error (mBias and mRMSE) at coarse scales.This assumption is only appropriate if systematic errors (mBias) are low at very fine scales and the residuals of the LiDAR-to-AGB STP model are normally distributed (Figure 4b).
Second, computing uncertainty in nonlinear models with measurement errors in both the response (σ 0 HV ) and predictor (AGB L ) variables is nontrivial.There is no one method for a robust error estimation, and common substitute methods include simulation or Monte Carlo-type tests [64].In this study, we conduct a sensitivity analysis by using different combinations of the dataset.We ran SAR-to-AGB L models on the extremes of the values of AGB L 's average error range (i.e., on AGB L +mRAE% and AGB L −mRAE%) to verify that the chosen best-fit models remained suitable for use at these ranges.Further, we simulated AGB values with an added random noise (AGB L + N (0, mRMSE)) and ran the models on these values to derive the standard deviation of prediction bias, RMSE and R 2 .We found that the chosen best fit remained the strongest fit at the extreme ranges AGB L +mRAE% and AGB L −mRAE%.Models with simulated AGB values showed an average of <0.1-Mg ha −1 change in AGB retrieval accuracy, primarily due to the larger number of observations available for regression.
Finally, σ 0 HV may vary widely between different radar scenes due to differences in acquisition dates and variations in environmental conditions across the landscape.Our SAR-to-AGB L model was tested on each scene separately.Chosen best fits remained the same, and trends in AGB retrieval uncertainty when comparing models across different spatial scales were found to be largely similar to using the whole dataset.However, there was higher variation in scenes with a smaller proportion of land at coarser scales due to fewer data points (Supplementary Material, Figure S9).Hence, it was reasonable to use the full dataset for model training and validation, particularly for coarse scales.

Biomass Retrieval Using LiDAR
Biomass retrieval models using LiDAR calibrated with data from the species trial plots showed a satisfactory level of precision for the objectives of this study, explaining 70% of the variation in AGB with an RMSE of 53.8 ± 6.4 Mg ha −1 .Some of the performance may be attributed to the model being trained on plots with even-aged trees, the same species composition, buffer-areas along perimeters and with forests of the same biophysical structure and management procedures.These factors can significantly reduce co-registration errors [69] and plot-edge errors [13,70] related to the use of LiDAR on small-sized plots.
Our study demonstrates the suitability of a two-step approach, in which radar images are calibrated against a higher-resolution LiDAR product, itself calibrated with ground measurements of AGB [14].It is important to discuss the use of such an approach and the associated errors in predicting AGB across Denmark and other forest types.Since the majority of Danish forests are plantations represented by the STPs [46], we can be confident that the model performs well in these areas.The model may however not perform as well over mixed forests types, particularly by under-predicting AGB in areas with AGB >300 Mg ha −1 , as observed in previous studies [49].This may result, since low-density LiDAR pulses (here, 0.5 pulses/m 2 ) are unable to sufficiently penetrate dense canopies and extract distinguishable structural characteristics in mixed forests.Hence, replicating our study in structurally complex or high biomass forests would need higher point-density LiDAR data or models built on larger plot sizes and varied sampling designs to improve AGB retrieval accuracy [14,19,69].Further, in averaging the LiDAR-derived AGB map (here, AGB L ) to coarser scales, we could not account for changes in model RMSE and bias due to the lack of large-area ground-data.Most studies have, however, shown marked improvements in accuracies as plot sizes increase [9,69,71,72], and hence, assuming no improvement with coarsening scale was considered a reasonable conservative approach in our study.Models built for forests with high spatial variability in vegetation structure would need to be scale-invariant [73] or need further analysis on the averaging of systematic errors with scale.

SAR Relation to Aboveground Biomass at 50 m-500 m Scales
As is common for most national, continental and global maps, our study involves building AGB-backscatter relations that cut across a complex landscape of patchy forests and agricultural lands, providing insight into AGB retrieval accuracy trends observed in areas dominated by human use.On coarsening resolution, the numerous small forest patches cause an increasing proportion of forest-edge pixels, each with high AGB variability within, but low overall AGB density.Trends across mapping scales may only be fairly compared after adequately removing these areas from the analysis (Section 3.2.1).The SAR-to-AGB L relation examined then demonstrated that mapping scale is a critically important parameter in AGB retrieval.Retrieval accuracy improved as the pixel sizes of the SAR image and LiDAR-derived AGB map used to calibrate it (AGB L ) were increased from 50 m to 250 m, and negligibly thereafter (Figure 9).A large part of the initial improvement (50 m to 100 m) is expected from the averaging of errors inherent in radar images, such as those from speckle, thermal noise, geolocation, canopy layover and variations due to moisture or topography [11].Improvements up to larger scales are also expected due to the reduction of local spatial variability of σ 0 HV and AGB by averaging, reducing the bias in calibration equations used to relate these variables [14].The negligible improvement beyond the 250-m scale indicates that there is an upper limit up to which mapping resolution needs to be coarsened before a large part of the underlying forest heterogeneity is lost and retrieval accuracy is constant.This is a relevant finding for studies that need to establish adequately-sized ground-plots to calibrate backscatter, acquired for example from the recently launched L-band ALOS PALSAR-2 mission, to AGB.Previous studies using up to 1-ha plots have noted improving AGB-backscatter relations with scale [11,74].Our study shows that plot scales of 150-250 m (i.e., sizes of 2.25-6.25 ha) may be most suitable for mapping AGB, significantly larger than most forest inventory plots (e.g., the 707-m 2 area plots of the Danish NFI).
Observing results of the SAR-to-AGB L model for regions of different biomass reveals that while AGB retrieval is more accurate at fine scales (50 m or 100 m) for areas with AGB below ∼110 Mg ha −1 , retrieval at coarser scales is more accurate for areas with higher AGB (Figure 10d).The AGB value at which backscatter begins to lose sensitivity to AGB in the models also increases with coarsening spatial scale as the data distribution changes (Figure 10a).A comparison of model R 2 's (Table 7) suggests that suitable mapping resolutions would be 100 m for low biomass areas and 250 m for high biomass areas in Denmark.These results are of relevance for studies aimed at AGB change detection using SAR, where accurate AGB estimates pre-and post-disturbances are desired [4].Establishing a statistically optimal plot size for data calibration and an optimal biomass mapping resolution are dependent on the range of biomass considered, desired error bounds and the distribution of forests being studied across landscapes.

SAR Relation to Vegetation Height and Cover
The backscatter relation to mean vegetation height displayed similar trends of signal saturation and improvement in model accuracy with scale as the AGB-backscatter relation.This result is expected, since a height distribution percentile (P 95 ) was used to derive AGB L , and height was thus representative of the AGB across the landscape.Mean height and other height percentiles were unable to explain any of the scatter in AGB-backscatter plot.
In contrast, relating the LiDAR-derived canopy interception ratio (VIR) to σ 0 HV showed no saturation and a strong linear relation in the power domain.If VIR is interpreted simply as the proportion of vegetation cover in a pixel, the result shows that as more signals are scattered off of vegetation, more energy is returned to the satellite sensor.However, this interpretation alone makes no distinction between backscatter returned from forests with different canopy or stem density.Since VIR does not provide a measure of the vertical component of canopy density, some interpretation of the relation is required and discussed here.SAR signal saturation is observed even in areas of low VIR (≤ 0.5), such as low stem density plantations and areas of low fraction of vegetation cover (Figure 12).This indicates that canopy opacity may not be the only cause of signal saturation, but other factors, such as the size of scattering elements (trunks and stems) may influence saturation [43,68].Further, a slight decreasing trend in σ 0 HV is observed as AGB increased on forested pixels with VIR > 0.75.If increasing AGB is associated with increasing height and correspondingly larger and more numerous branches encountered in the vertical dimension, then signal attenuation may be associated with increased vertical canopy density.For long-wavelength SAR, changes in both the cross-sectional area and number density of branching elements with depth into the canopy have been theoretically shown to influence backscatter [43,45].Our LiDAR data alone were not suitable to study the depth of the attenuating layer due to its relatively low point density.The analysis and discussion regarding VIR only provide a preliminary and demonstrative insight into the degree to which interaction with canopies explains variability in the AGB backscatter plot.To investigate how information on forest structural variables may improve AGB retrieval, however, requires mapping AGB independent of VIR and adding species-specific structural information to the models.

Summary and Conclusion
In this study, we have provided a characterization of the relations between cross-polarized ALOS PALSAR backscatter and LiDAR-derived vegetation cover, height and aboveground biomass at multiple map resolutions over Denmark.The study uses unique stepping-stone datasets, including a high-resolution DEM, ground-measured plots over even-aged plantations, country-wide airborne LiDAR, aerial photography and land use maps to develop these relations through empirical models.
Our key results show that in complex landscapes dominated by patchy forests and human use, optimal scales to retrieve and map AGB from SAR are approximately 100 m for areas with low (< 110 Mg ha −1 ) and 250 m for areas with high AGB.Improvement in retrieval at scales coarser than 250 m is negligible.Furthermore, backscatter loses sensitivity to height and biomass where these variables are relatively high, but showed a linear relation to a metric measuring the fraction of vegetation cover at >1 m height and LiDAR signal penetrability (VIR).VIR explained some of the observed scatter in the biomass-backscatter relation at very coarse scales (≥ 250 m), such that in areas of high VIR (≥0.75), there is a slight decline in backscatter as biomass increased.
Broadly, our study reinforces that spatial scale and vegetation structural properties are essential to consider when using SAR backscatter as an indicator of biomass.Particularly, for large-scale national or continental maps that cut across heterogeneous or dense vegetation types, information on the distribution of forest cover and structure will support accurate AGB retrieval.Airborne or spaceborne LiDAR surveys can readily provide such information, and future products from ICESat-2, GEDI and MOLI will hence be an important data source for AGB mapping.Further examination of the relationships explored in the study is also certain to provide a stronger basis for designing protocols to map vegetation biomass using radar data from ALOS PALSAR-2 and the future BIOMASS mission.

Figure 1 .
Figure 1.(a) Forested areas in Denmark (map provided by Geodatastyrelsen), locations of the species trial plots (STPs) in clusters and the extent of the ALOS PALSAR scenes and LiDAR strips that were removed before analysis.(b,c,d) An example of a forest plantation in Denmark with extracted LiDAR variables, vegetation interception ratio (VIR), mean height (MH) and AGB L , at 50 m× 50 m pixel resolution.Built areas and water bodies were removed from all maps prior to analysis.

Figure 2 .
Figure 2. Example of a cluster of STPs, overlaid on a high-resolution (0.16 m × 0.16 m) airborne optical image obtained during leaf-on season in 2010 by COWI.Each STP contains a single tree species and is placed within the center of a forest stand.

Figure 3 .
Figure 3.A hypothetical demonstration of the variation of vegetation interception ratio (VIR) with vegetation distribution and aboveground biomass (AGB).(a) Pixels with different AGB density (due to varying tree height), but similar VIR.(b) Pixels with similar AGB density, but varying VIR due to the different distribution of stems.

Figure 4 .
Figure 4. (a) Results of the LiDAR-to-AGB STP model (Equation (2)) predicting AGB STP with the 95th height percentile multiplied by VIR.(b) Normal probability plot of residual distribution.Normality was assessed using the Shapiro-Wilk test (p = 0.17).(c) Predicted vs. measured biomass for STPs of different sizes demonstrate the lack of dependence of bias on plot size.

Figure 5 .
Figure 5. Distribution of AGB L for forested areas at different spatial scales (50 m to 500 m, indicated on each sub-figure).Frequency (y-axis) refers to the number of observations in each range of AGB (x-axis).

Figure 6 .
Figure 6.SAR-to-VIR L model in the power domain (Left panel) and dB (Right panel) for the training dataset at a few of the tested spatial scales in this study (50 m to 500 m, indicated on each sub-figure).Data are represented on a 2D histogram density plot, with values byte-scaled (0 to 255) to the color bars shown.

Figure 7 .
Figure 7. SAR-to-AGB L model (Equation (4)) (Left panel) and the performance of the inverted model (Equation (5)) (Right panel) used to retrieve AGB at multiple spatial scales (50 m to 500 m, indicated on each sub-figure).Observations (N) are represented on a 2D histogram density plot, with values byte-scaled (0 to 255).High density observations that run across the residual plots at ∼130 Mg ha −1 show where a maximum retrieval value (MRV) was assigned to predicted AGB when σ 0HV fell above the predictable range.

Figure 8 .
Figure 8. SAR-to-MH L model at multiple spatial scales (50 m to 500 m, indicated on each sub-figure).Observations (N) are represented on a 2D histogram density plot, with values byte-scaled (0 to 255).

Figure 9 .
Figure 9. AGB retrieval accuracy ((Left panel), RMSE; (Centre panel), relative RMSE; (Right panel), bias) from σ 0HV as a function of spatial scale before and after applying a maximum retrieval value (MRV) to predicted AGB.

Figure 10 .
Figure 10.AGB retrieval bias (a,b) and RMSE (c,d) at different ranges of AGB using the SAR-to-AGB L model.AGB bins of 50 Mg ha −1 were analyzed individually through the dataset and the prediction bias and RMSE computed for each bin.Relative prediction bias (b) is calculated as the ratio of the average value of the difference between predicted and measured AGB, to average measured AGB, in each bin.Relative prediction RMSE (d) is calculated as the ratio of RMSE to average measured AGB in each bin.

Figure 11 .
Figure 11.Residual σ 0 HV from the SAR-to-AGB L model plotted against VIR (a) and MH (b).Observations (N) are represented on a 2D histogram density plot, with values byte-scaled (zero to 255).

Figure 12 .
Figure 12.AGB L and σ 0 HV relation for areas of different VIR across the country at 500-m scale (VIR ranges are indicated on each sub-figure).The power (m 2 /m 2 ) scale is used for σ 0 HV to demonstrate the linear and decreasing relation at high σ 0 HV ranges.All model coefficients were significant (p-value < 0.001).N indicates the number of observations.

Table 2 .
Summary of the airborne LiDAR and flight data collected over Denmark in spring

Table 3 .
List of ALOS PALSAR scenes and acquisition dates used in the study.All scenes were taken in fine beam dual polarization mode (HH + HV) in ascending orbit.

Table 4 .
Summary of the results of the LiDAR-to-AGB STP model (Equation (2)).

Table 5 .
Summary of the LiDAR-to-AGB STP model validation using jackknife resampling.R i and P i are the ground-measured and model-predicted values for sample i, respectively, and n is the number of samples.

Table 6 .
Results of SAR-to-VIR L model for the validation dataset.dandg are regression coefficients used in Equation (3), reported with the standard error.The last two columns show regression results when building the model with σ 0 HV in dB.Significance is provided as * p-value < 0.01.

Table 7 .
Results of the SAR-to-AGB L model.a,b and c are regression coefficients used in Equation (4), reported with the standard error.Significance is provided as * p-value < 0.01.