Estimation of Forest Growing Stock Volume with Synthetic Aperture Radar: A Comparison of Model-Fitting Methods

Abstract

Satellite-based estimation of forest variables including forest biomass relies on model-based approaches since forest biomass cannot be directly measured from space. Such models require ground reference data to adapt to the local forest structure and acquired satellite data. For wide-area mapping, such reference data are too sparse to train the biomass retrieval model and approaches for calibrating that are independent from training data are sought. In this study, we compare the performance of one such calibration approach with the traditional regression modelling using reference measurements. The performance was evaluated at four sites representative of the major forest biomes in Europe focusing on growing stock volume (GSV) prediction from time series of C-band Sentinel-1 and Advanced Land Observing Satellite Phased Array L-band Synthetic Aperture Radar (ALOS-2 PALSAR-2) backscatter measurements. The retrieval model was based on a Water Cloud Model (WCM) and integrated two forest structural functions. The WCM trained with plot inventory GSV values or calibrated with the aid of auxiliary data products correctly reproduced the trend between SAR backscatter and GSV measurements across all sites. The WCM-predicted backscatter was within the range of measurements for a given GSV level with average model residuals being smaller than the range of the observations. The accuracy of the GSV estimated with the calibrated WCM was close to the accuracy obtained with the trained WCM. The difference in terms of root mean square error (RMSE) was less than 5% units. This study demonstrates that it is possible to predict biomass without providing reference measurements for model training provided that the modelling scheme is physically based and the calibration is well set and understood.

1. Introduction

Knowledge of the spatial distribution and magnitude of forest resources is of paramount importance to understand pathways of climate, to support the economy of communities, and to provide societal benefits [1]. Forest resources are best quantified with a variable that expresses the biomass stored by vegetation, this being a direct proxy for its carbon content. Estimation of forest biomass defined either as the structural variable growing stock volume (GSV, unit m³ ha⁻¹) or as organic mass of dry matter, i.e., above ground biomass (AGB, unit Mg ha⁻¹), requires measurements of both the horizontal (i.e., basal area) and vertical (i.e., tree height) properties of a forest, as well as the number of trees, to obtain the biomass density. To estimate the volume, in addition, the tree form factor is needed whereas the organic mass can be estimated from the wood volume by accounting for the wood density and the total-to-stem biomass ratio [2]. Field measurements provide accurate information locally. Such measurements, however, are not easily scalable and require statistical tools to obtain spatially explicit estimates [3].

Satellite-based active microwave remote sensing enables frequent and systematic monitoring of the Earth’s surface. Over vegetated land, the observations contain information on both the structure and water content of plants and soil [4]. Identifying each component in the microwave signal is therefore essential for quantifying the quantification of the biomass stored in the vegetation. Typically, this is achieved by preferring observations where the contribution of the forest structure to the recorded signal is predominant. Research demonstrated that data acquired by Synthetic Aperture Radar (SAR) operating at longer wavelengths such as in the Very High Frequency (VHF) band [5] and P-band [6] are more strongly related to forest biomass than data acquired by SAR sensors operating at shorter wavelengths ([7,8]). In space, however, SAR sensors have operated only up to the L-band (23 cm) so far, so that the sensitivity to biomass of the SAR backscattered intensity in satellite images is mostly moderate or weak due to the limited capability of the microwave signal to penetrate the canopy and the influence of environmental conditions and terrain slope on the backscattered signal [9,10,11,12,13,14]. Stronger agreement with biomass was obtained with polarimetric features [15,16,17,18]. To overcome the limitations intrinsic in the data, a wide range of parametric and non-parametric model-based approaches that aimed at estimating biomass from SAR images have been developed and tested [9,10,19]. With the availability of SAR image datasets acquired at different frequencies and polarizations, these approaches have been advanced to improve the accuracy of the estimation by exploiting complementarity of the observations, for example, through multiple linear regression techniques [20,21,22], physical-based models [23,24,25], or machine learning [20,26,27,28].

While the choice of the estimation approach, i.e., whether parametric or non-parametric, has mostly shown only a marginal impact on the estimation accuracy [20,29], retrieval approaches based on machine learning techniques have received more attention in the recent literature. This is due to their capability to handle diverse datasets and provide a reliable estimate of the biomass even at levels when the SAR predictors are reported to saturate [27,30,31,32,33,34,35,36,37,38]. Although the variance of the estimates was mostly limited, they were often erroneous at the tails of the biomass distribution, with over- and underestimation affecting the low and high biomass range, respectively. These issues were overcome in a few studies in which either the predicting features were first optimized with a pre-selection based on seasonal conditions [33] or forest types [38] or supported by a genetic algorithm [39]. Increasing the feature space of predictors with additional SAR backscatter values simulated with a Water Cloud Model was instead proven to stabilize the retrieval based on Neural Networks across sites [26].

One common observation in the cited literature on the effective use of non-parametric estimation models is that they typically require a densely populated dataset of biomass values that act as a reference to establish the relationship between the predictors and the biomass variable of interest. This aspect is crucial, especially when satellite measurements have limited sensitivity to forest biomass (as currently available satellite sensors are). In addition, these studies were undertaken at local sites. While a model trained on field reference data can predict biomass within forests for which the reference data are representative, the extension of the same model to dissimilar areas might be unreliable [26]. The scarcity of ground reference data on forest biomass for most regions of the world means that wide-area model-based biomass mapping from satellite imagery, e.g., continents, requires strong generalizations during the model calibration phase to prevent unrealistic biomass retrievals. Eventually, this results in strongly biassed estimates of biomass [40,41]. The high spatial sampling and spatial resolution of LiDAR-based estimates of AGB from the Global Ecosystem Dynamics Investigation (GEDI) mission have recently been exploited to train retrieval models [25,38] and to derive nation-wide maps of AGB [35]. Although the results were not as accurate as when the retrieval model was trained with field plot measurements due to the superior quality of the latter, spaceborne LiDAR-based datasets may have a prominent role in future studies as their accuracy will improve.

As an alternative to calibrating the retrieval model with a training dataset of GSV or AGB values, the parameters of parametric models can be derived from a statistical analysis of the SAR image used as a predictor and masked for given forest conditions [23,24]. Despite several approximations in how the retrieval model is calibrated, the backscatter predicted by the model was within the range of the observations. These results suggested that the calibration method, referred to as BIOMASAR, could support large-scale mapping of biomass using spaceborne SAR backscatter observations at C- and L-bands [42] and thus overcoming the issues of limited availability of samples to train the retrieval model. Although BIOMASAR is used to create wall-to-wall estimates of GSV [43] and AGB [24,44], the retrieval performance with a calibrated model has been seldom explored elsewhere [45] and an assessment of factors that impact the calibration of the retrieval model is still missing.

In this study, we therefore compared the retrieval of forest GSV from SAR backscatter data acquired at C- and L-bands using the Water Cloud Model (WCM) trained on one side with field measurements and calibrated on the other side with the aid of auxiliary data products. The study was undertaken at multiple study sites representative of boreal, temperate, and Mediterranean forest landscapes in Europe to disentangle the effect of the training dataset, SAR backscatter measurements, forest structure, and modelling framework.

The focus is on GSV because it is a key parameter in the context of forest resource management in Europe [1] and major proxy of AGB and carbon stored in vegetation [46,47,48,49,50] to address global climate issues [51]. GSV represents the volume of the tree stems for all living species per unit area, including bark but excluding branches and stumps. The study sites and the data supporting our study are presented in Section 2 and Section 3, respectively. For each site, plot-level GSV data from field inventories were used to evaluate the performance of the calibration and validate the GSV estimates (Section 3.1). The study in addition utilized sub-national averages of GSV values reported by National Forest Inventories (NFIs), which are easily accessible and provide information on the spatial distribution of GSV at a continental scale (Section 3.2). Sub-national averages were considered to simulate scenarios where GSV reference values, such as plot inventory data, are unavailable. The proposed framework for GSV retrieval foresaw the use of SAR backscatter acquired by Sentinel-1 (C-band) and Advanced Land Observing Satellite 2 (ALOS-2) Phased Array L-band SAR 2 (PALSAR-2) as predictors (Section 3.3) and canopy height metrics based on satellite LiDAR observations as an additional source to calibrate a set of forest structural functions that were integrated in the main retrieval model (Section 3.4). Our analysis focused on (i) understanding the relationship between SAR backscatter values and GSV, and (ii) discussing the training and the calibration of the WCM. By computing the accuracy of the GSV estimates from the trained and calibrated WCM with respect to the plot-level data, our aim was to assess the capability of a WCM to estimate GSV in a scenario when no reference data are available to train the model.

2. Study Sites

The study sites (Figure 1) included forests from major European ecological zones and different types of topographic conditions from gently undulating landscapes to mountains (

Table 1. Summary of study sites and corresponding field measurements.

Site (Geographic Location)	Biome	Topography	Field Measurements
			# Sample Plots	Year of Inventory	Ownership
Finland N	Boreal forests and taiga	Hilly	1004	2018	Public
Finland S	Boreal forests and taiga	Gently undulating	1064	2018	Public
Catalonia	Mediterranean woodlands, forests, and scrub	Hilly to mountainous	663	2016	Public
Romania	Temperate broadleaf and mixed forests	Hilly to mountainous	1306	2019	Private (Tornator Oyj)

). The Finland N and Finland S sites were located north of the Arctic Circle and by the Baltic Sea, respectively. At both sites, the forests consisted of semi-natural broadleaf and coniferous forest. The Finland N site was characterized by a hilly landscape with scattered open peatlands. Some of the hilltops had very sparse forests or no trees at all. The Finland S site was gently undulating and consisted of partly fragmented forest and agricultural landscapes. The Catalonian site was in the east of Catalonia, by the Mediterranean Sea. The forests consisted of semi-natural broadleaf and coniferous forest on hilly to mountainous landscape. The Romanian site was in the Carpathian Mountains. The forests consisted of semi-natural broadleaf (mainly beech and oak) and coniferous (mainly spruce and fir) species on hilly to mountainous landscape.

3. Material

3.1. Plot-Level GSV

Field measurements conducted by the Finnish Forest Centre (https://www.metsakeskus.fi/en, last accessed on 4 April 2024) at sample plots were used as a reference for the two Finnish sites. The plots were arranged in clusters, each containing six plots and chosen to cover the variation in forest types and site conditions. Three different plot radii were used in the sampling: 9 m in young and advanced managed forests with relatively large tree densities; 12.62 m in forest with small stem densities but usually a large volume due to the mature development stage, and 5.64 m in seedling stands. For the Catalonian site, sample plot data measured by the Spanish NFI were used and consisted of circular sample plots with a radius varying from 5 to 25 m. For the Romanian site, field plot measurements were provided by a private company, which had contracted a field inventory on their forest area to support forest management decisions. The sample plots were circular and covered areas of 500 m² (radius of 12.6 m).

Among several features, the database of measurements from each site reported values for the basal area, tree height, diameter at breast height, growing stock volume, site type, and species proportions (pine, spruce, and broadleaf). The statistical distribution of the GSV from the field measurements is summarized in

Table 2. GSV statistics for each site.

Site	GSV (m3 ha−1)
	Average	Quartiles (1, 2, 3)	Min/Max
Finland N	95	44/83/135	0/498
Finland S	157	47/129/233	0/751
Catalonia	108	55/92/144	1/480
Romania	418	204/379/582	0/1677

. The overall largest GSVs were measured in Romania. GSV at the two Finnish sites increased from north to south. The Catalonian site was mostly characterized by moderate-to-low GSV, with some patches of very high GSV.

3.2. Sub-National Average GSV Values

European NFIs regularly publish forest variable statistics at administrative unit levels, thus providing comprehensive information on the spatial distribution of biomass in Europe, which otherwise would not be attainable because the original measurements at the plot level are usually not accessible [52]. For this study, GSV average values were gathered from reports and table statistics published by the NFIs for 449 administrative or jurisdictional regions in 28 countries (Table S1). The reporting period is heterogeneous and spans the entire 2010 decade. An update with growth factors to the study years was not undertaken to unnecessarily bias the values since the effect of disturbances could not be accounted for.

3.3. SAR Dataset

The SAR dataset consisted of images acquired at the C-band (5.6 cm wavelength) by the European Space Agency (ESA) Sentinel-1 satellite and at the L-band (23 cm wavelength) by Phased Array L-band Synthetic Aperture Radar 2 (PALSAR-2) on Advanced Land Observing Satellite 2 (ALOS-2) operated by the Japan Aerospace Exploration Agency (JAXA, Tokyo, Japan). In correspondence with the area of each sample plot, we computed the area-weighted backscatter from the pre-processed Sentinel-1 and the ALOS-2 PALSAR-2 images.

3.3.1. Sentinel-1 Dataset and Pre-Processing

The Sentinel-1 dataset consisted of all images acquired during the same year as the corresponding forest inventory (Table 1). Images were acquired in the Interferometric Wide Swath mode (IWS) in VV and VH polarization. For the sites of Catalonia, Finland N, and Romania, 59 dual-polarized pairs of images were available. For the site of Finland S, the dataset included 45 dual-polarized image pairs. The Sentinel-1 imagery consisted of detected amplitudes of the radar backscattered intensity and was obtained in ground-range geometry, with a pixel size of 10 m in the range and azimuth direction. Pre-processing of the SAR images followed a standard approach and consisted of (i) multi-looking to 20 m to be closer to the native spatial resolution of the data, (ii) terrain geocoding using the 30 m Copernicus Digital Elevation Model (DEM) as a reference [53] with an automated matching procedure [54], and (iii) correction for the pixel area effect on the radar backscatter [55]. For each site, the Sentinel-1 images were clipped and geocoded to a common geometry and then filtered for speckle with a multi-channel approach [56].

3.3.2. ALOS-2 PALSAR-2 Dataset and Pre-Processing

The ALOS-2 PALSAR-2 dataset consisted of yearly mosaics of the radar backscatter at HH and HV polarization with a pixel size of 25 m × 25 m [57]. Mosaics have been generated by JAXA for each year since 2015.

The mosaics are based primarily on ALOS-2 Fine Beam Dual-polarization (FBD) images acquired between May and October of a given year. Each of the mosaics is provided in the form of 1° × 1° tiles and includes the (i) two backscatter images; (ii) an image of the local incidence angle with respect to the orientation of the pixel, derived from a DEM (3-arcsec Shuttle Radar Topography Mission (SRTM) or 1-arcsec ASTER DEM); (iii) a radar layover/shadow mask; (iv) an image of the date of acquisition; and (v) a mask detailing whether a pixel is labelled as land or water. Prior to mosaicking, the SAR images were geocoded, orthorectified, and calibrated by JAXA. The mosaics were compensated for variations in the pixel scattering area due to topography and for the dependence of backscatter on the local incidence angle [58]. Although the pre-processing of the ALOS-2 PALSAR-2 and the Sentinel-1 images followed the same sequence, the geocoding and terrain normalization procedures implemented different approaches. The visual inspection of the mosaics revealed partially uncompensated slope-induced effects on the backscatter. The SAR images were finally resampled to a pixel size of 1/4000th of a degree in both latitude and longitude, corresponding to roughly 25 m at the Equator.

For this work, we used version 2.1 of the mosaics and resampled them with bilinear interpolation to the geometry of the Sentinel-1 images. Besides the mosaic corresponding to the year of the field inventory data, we also considered the other mosaics to understand the impact of an L-band multi-temporal dataset on the estimation of GSV.

3.3.3. Uncertainty of the SAR Backscatter

The uncertainty of a SAR backscatter measurement was represented by a measure of its precision, i.e., its standard deviation. It was estimated empirically with the Equivalent Number of Looks (ENL) [59] where μ represents the expected value of the SAR backscatter and σ its standard deviation.

$$ENL={\displaystyle \frac{{\mu }^2}{{\sigma }^2}}$$

(1)

To compute the ENL, we superimposed a grid with a window size of 100 × 100 pixels on the stack of the SAR images and computed Equation (1) in each window and for each image. To avoid that heterogeneous landscapes or sloped terrain would impact the estimation, the ENL was defined as the 90th percentile of all values obtained for the grid cells. To account for the seasonal variability of C-band backscatter [42,60] and obtain an ENL representative of the Sentinel-1 dataset at a given site, we computed the median of the ENL estimates for the individual images.

Table 3. Estimates of the ENL from the Sentinel-1 and ALOS-2 PALSAR-2 datasets.

Site	ENL
	Sentinel-1		ALOS-2 PALSAR-2
	VV	VH	HH	HV
Finland N	9	14	13	9
Finland S	6	13	8	7
Catalonia	4	7	8	9
Romania	4	11	8	9

lists the estimates of ENL for each polarization and site for both sensors. The lower signal-to-noise ratio of the VH-polarized channel compared to the VV-polarized channel implied that the ENL of the Sentinel-1 VH-polarized data was higher than for VV polarization. The VV-polarized ENL of Sentinel-1 was rather constant between 4 and 6 except for the Finland N site where the estimate was slightly higher. The VH-polarized ENL for Sentinel-1 was also rather constant, between 11 and 14, except for Catalonia. For the ALOS-2 PALSAR-2 dataset, the ENL was mostly between 7 and 9, with the exception of Finland N where we estimated an ENL of 13 for the HH polarization. Compared to the Sentinel-1 dataset, the ENLs were lower. Although the ALOS-2 mosaics were obtained by multi-looking the original Single Look Complex images, the multi-channel filtering could not be applied.

Finally, a single global value for each sensor was computed by taking the median of the site-specific estimates. For Sentinel-1, the ENL was set as equal to 6 for the VV polarization and 11 for the VH polarization, corresponding to a standard deviation of 1.5 dB and 1.1 dB, respectively. For the ALOS-2 PALSAR-2 dataset, the ENL was set as equal to 8 for both the HH and HV polarizations, corresponding to a standard deviation of 1.3 dB.

3.4. ICESat-2 Dataset and Pre-Processing

The Advanced Topographic Laser Altimeter System (ATLAS) onboard the Ice Cloud and Elevation Satellite (ICESat-2) is a photon-counting instrument that splits the laser into six beams. Beams are arranged in pairs approximately 3.3 km apart [61], each pair consisting of a strong and weak energy beam (4:1 ratio). Geophysical parameters related to vegetation and terrain heights are available in the ATL08 product as data segments with a 100 m step size along the flight direction [62]. From the original photon data, the original files were reformatted into segments of a 100 m length and 25 m width. Data are available starting at October 2018; for this study, observations from 2019 were used. To reduce the number of segments with potentially erroneous height metrics, we discarded segments that were (i) acquired by the weak beams because of the very small number of photons recorded from the forest floor, in particular in dense canopies; (ii) characterized by less than three photons reflected by the canopy; (iii) flagged as not belonging to natural vegetation in the metadata; and (iv) exhibited an elevation differing by more than 25 m from the reference DEM used in all ATL products. In addition, values larger than 50 m were filtered out because they frequently occurred in sparsely vegetated or unvegetated regions. As suggested by the product developers, the canopy height was represented by the 98th percentile of all height values within a segment (RH98).

From this dataset, an estimate of the maximum canopy height was computed for each site to constrain the retrieval of GSV (

Table 4. Maximum canopy height for each site from the field inventory data and the ICESat-2 metric.

Site	Field Inventory		ICESat-2
	Height [m]	# Plots	Height [m]	# Segments
Finland N	23	1005	21	26,049
Finland S	35	1065	28	18,867
Catalonia	33	663	30	26,751
Romania	44	1306	42	10,900

). This parameter was set as equal to the 99th percentile of all valid canopy height values to avoid that individual unrealistic values would bias the estimate if set as equal to the largest of all canopy height values. The difference between the ICESat-2-based values and the maximum canopy height derived from the field inventory dataset (Table 4) was small at Finland N, Catalonia, and Romania. At Finland S, the discrepancy was due to the different sampling in either the ICESat-2 dataset or the plot dataset, respectively, in correspondence with the tallest forests. The uncertainty of the maximum height was assessed by computing the same value with ICESat-2 data acquired in 2020 and 2021. The difference between the estimates of maximum canopy height for the three years was less than 1 m so that the uncertainty associated with this variable was considered negligible in this study.

4. Methods

4.1. Forest Backscatter Model

The relationship between the SAR backscatter and GSV was expressed with the Water Cloud Model with gaps [63] integrated with forest structural functions [64]. The original WCM described the backscattered intensity from a forest as a function of the backscatter from the forest floor through gaps in the canopy, the backscatter from the forest floor attenuated by the canopy, and the backscatter from the canopy.

$${\sigma }_{for}^0=\left(1-\eta \right){\sigma }_{gr}^0+{\eta \sigma }_{gr}^0{T}_{tree}+{\eta \sigma }_{veg}^0\left(1-T_{tree}\right)$$

(2)

The forest backscatter was expressed as a function of two structural parameters describing horizontal and vertical properties of the vegetation. The variable η represented the canopy density from a microwave perspective; as a first approximation, we considered it equal to the canopy density in the optical domain even though the gaps should decrease for an increasing wavelength. The two model parameters σ⁰_gr and σ⁰_veg represent the backscattering coefficient of the ground and vegetation layer, respectively. T_tree represents the two-way tree transmissivity and was expressed in terms of a two-way attenuation per metre through a tree canopy, α, and the depth of the attenuating layer, which was assumed to correspond to the canopy height, h.

$$T_{tree}=e^{-\alpha h}$$

(3)

To express the backscatter from a forest as a function of GSV, we used two structural functions relating canopy density, canopy height, and GSV. Using spaceborne LiDAR metrics of canopy height, h, and density, CD, the exponential model in Equation (4) was demonstrated as a reliable fit to the data both at regional [64] and global [65] scales. The coefficient q was of an empirical nature.

$$CD=1-e^{-qh}$$

(4)

The power-law relationship between canopy height, h, and GSV, V, in Equation (5) with a and b being coefficients of an empirical nature was demonstrated when relating either field measurements or LiDAR-based metrics of canopy height to GSV [63,64].

$$V=a·h^b$$

(5)

With η = CD, the integrated GSV retrieval model was therefore expressed as

$$\begin{gathered} {\sigma }_{for}^0={\sigma }_{gr}^0\left(e^{-q{\left(aV\right)}^b}+e^{-\alpha {\left(aV\right)}^b}-e^{-\left(q+\alpha \right){\left(aV\right)}^b}\right) \\ +{\sigma }_{veg}^0\left(1-e^{-q{\left(aV\right)}^b}-e^{-\alpha {\left(aV\right)}^b}+e^{-\left(q+\alpha \right){\left(aV\right)}^b}\right) \end{gathered}$$

(6)

4.2. Estimation of Model Parameters

The forest backscatter model in Equation (6) contained six parameters. The structural parameters q, a, and b were estimated from a least squares regression between LiDAR-based metrics of canopy density, canopy height, and GSV. The WCM parameters σ⁰_gr, σ⁰_veg, and α described the interaction of the microwave signal with the woody vegetation and the soil underneath and were estimated from the SAR backscatter data following different strategies.

The estimation of the coefficient q in Equation (4) required values of canopy density and canopy height as, for example, obtained from the ICESat LiDAR measurements [66]. Here, we used the spatially explicit dataset of the coefficient presented in [66] and associated with each site and the most frequent estimate of q. A site-specific estimate based on footprint-level data was discarded because of the sparse sampling by ICESat. The use of more recent LiDAR observations was not considered because it either would not cover all sites (Global Ecosystem Dynamics Investigation, GEDI) or would not provide a suitable canopy density metric (ICESat-2).

The coefficients of the structural function in Equation (5) were estimated by regressing to the sub-national averages of GSV and the corresponding averages of LiDAR-based canopy height. To account for the spatial variability of the association between height and volume, the equation was fitted for each of the four largest terrestrial biomes in Europe according to [67], i.e., “temperate broadleaf and mixed forests”, “temperate conifer forests”, “boreal forests and taiga”, and “Mediterranean forests, woodlands, and scrub”. The smaller biomes “temperate grasslands, savannas, and scrublands” and “tundra” were merged with temperate broadleaf forests and boreal forests, respectively. Estimates of the two coefficients for the respective biome were assigned to each site accordingly (Table 1). Although measurements of forest height and GSV were available for each plot and at each site, we did not create site-specific estimates because this study aimed at assessing the performance of the retrieval model in a scenario where training data in the form of plot measurements are not available.

Having set the two structural functions’ set, the three unknown parameters of the WCM are estimated. Traditionally, the WCM is regressed to a training set of GSV measurements and corresponding backscatter values. In practice, two major issues are associated with the regression.

• The weak sensitivity of the C- and L-band backscatter to GSV implies that a model fit with three degrees of freedom may be characterized by large uncertainties. It is also likely that some of the physical constraints of the three model parameters are violated so that they may become pure regression parameters. This applies especially to the coefficient α that governs the rate of change of the backscatter with GSV.
• When the training data are not available or the area of interest is far away from the location of the training data, it is likely that the estimates of the three parameters are not representative of the area.

To overcome these issues, the estimation of the three parameters was performed in two steps. Firstly, the coefficient α was estimated. Values and estimates of tree transmissivity reported in the literature are scarce and uncertain [68]. The only alternative was to rely on the generic values reported in the literature and corroborate them with local estimates by means of a least squares’ regression of the WCM with a training set of volume and backscatter measurements to capture variations dependent on the frequency, polarization, look angle, and seasonal conditions. Given an estimate of α, the other two parameters, σ⁰_gr and σ⁰_veg, were estimated either with a regression to a training dataset or with calibration using auxiliary data.

In absence of a training dataset, one approach to estimate the model parameters σ⁰_gr and σ⁰_veg was based on image statistics in areas where backscatter stems from either unvegetated or densely vegetated terrain, respectively [23]. Since σ⁰_gr represented the backscatter from an unvegetated surface, it was reasonable to assume that an estimate might be derived from the backscatter measurements observed over surfaces that are unvegetated. A mean or median value from such measurements was then associated with σ⁰_gr. The parameter σ⁰_veg represented the backscatter of the fully opaque canopy. The backscatter for the densest forests within the area of interest was the type of measurement closest to the definition of the model parameter. However, the measurements had to be corrected for an offset that removed the residual ground contribution since the measured backscatter contained a contribution from the soil underneath the canopy that is visible through the canopy gaps or is attenuated by the canopy.

The compensation of the backscatter for such an offset used the model in Equation (1) and expressed σ⁰_veg as a function of the backscatter representative of dense forests, σ⁰_df; the canopy density of dense forests, η_df; and the corresponding canopy height of dense forests, h_df, which is contained in the tree transmissivity term T_tree,df, together with the attenuation coefficient α and the backscatter coefficient σ⁰_gr.

$${\sigma }_{veg}^0={\displaystyle \frac{{\sigma }_{df}^0-{\left(1-{\eta }_{df}+{\eta }_{df}{T}_{tree,df}\right)\sigma }_{gr}^0}{{\eta }_{df}\left(1-T_{tree,df}\right)}}$$

(7)

In previous developments, σ⁰_gr and σ⁰_df were estimated by selecting pixels in correspondence with sparse and dense forest cover according to a map of canopy density, e.g., [69], and calculating the mean observed backscatter in such areas, respectively. The estimation of σ⁰_veg from σ⁰_df with Equation (7) required an estimate of the canopy density and canopy height of dense forests. In [43], the values corresponded to the 90th percentile of the histogram of each variable derived from waveform data acquired by ICESat GLAS. This percentile minimized the difference in the estimate of σ⁰_veg from Equation (7) and the estimate of σ⁰_veg when fitting the backscatter model to field measurements of biomass. The empirically defined sensor- and region-specific thresholds to delineate sparse and dense forest cover and a definition of η_df and h_df based on LiDAR data, independent of the definition of σ⁰_df, implied a risk of obtaining biassed estimates of σ⁰_veg and σ⁰_gr locally.

To reduce the impact of uncertainties on the estimation of σ⁰_gr and σ⁰_veg, we investigated here an alternative calibration approach that regressed between backscatter measurements and corresponding canopy density values with Equation (1). To express the backscatter as a function of canopy density only, the structural function in Equation (4) was first inverted. By substituting canopy height with a function of canopy density in Equation (1), the WCM was rewritten as

$${\sigma }_{for}^0=\left(1-\eta \right){\sigma }_{gr}^0+{\eta \sigma }_{gr}^0{e}^{\alpha {\displaystyle \frac{\mathrm{l}\mathrm{o}\mathrm{g}\left(1-\eta \right)}{q}}}+{\mathsf{\eta }\sigma }_{veg}^0\left(1-e^{\alpha {\displaystyle \frac{\mathrm{l}\mathrm{o}\mathrm{g}\left(1-\eta \right)}{q}}}\right)$$

(8)

This simplification was motivated by the availability of canopy density map products [69,70], a variable that is closely related to observations by optical remote sensing. The advantage of this calibration approach over the previous one was the independence from empirical thresholds used to select “ground” and “dense forest” pixels. The caveat was that uncertainties affecting the canopy density observations would have propagated to the estimates of σ⁰_gr and σ⁰_veg.

According to its definition, the estimate of σ⁰_veg from the regression of Equation (8) to observations of the SAR backscatter and canopy density expressed the backscatter of a canopy density of 100%. However, a canopy density of 100% corresponded to a large range of canopy heights and GSV values. Assuming that the backscatter for a canopy density of 100% was somewhat sensitive to forest structure, the initial estimate of σ⁰_veg from the regression was offset with a term here represented by the standard deviation of the SAR backscatter for the largest canopy density.

$${\sigma }_{veg}^0={\widehat{\sigma }}_{veg}^0+{2\ast SD}_{CD=100\%,meas}$$

(9)

In Equation (9), ${\widehat{\sigma }}_{veg}^0$ represents the initial estimate of σ⁰_veg obtained from the regression and ${SD}_{CD=100\%,meas}$ represents the standard deviation of the backscatter for a canopy density of 100%. The factor 2 was used to simulate a confidence interval from the standard deviation; herewith, it was assumed that the backscatter measurements for a given canopy density range were normally distributed. Since the observed standard deviation consisted of both the measurement noise and signal related to biomass variability, the term ${SD}_{CD=100\%}$ was obtained by subtracting the measurement noise from the observed value ${SD}_{CD=100\%,meas}$:

$${SD}_{CD=100\%}=\sqrt{{\left({SD}_{CD=100\%,meas}\right)}^2-{\displaystyle \frac{{{\widehat{\sigma }}_{veg}^0}^2}{ENL}}}$$

(10)

The measurement noise required the original measurement representative of the backscatter for a canopy density of 100%, ${\widehat{\sigma }}_{veg}^0$, and the ENL.

4.3. Retrieval of GSV

GSV was estimated from the SAR backscatter with the WCM in Equation (6). Since this form of the WCM could not be inverted mathematically, an estimate of GSV was obtained by minimizing the difference between the measured backscatter and the modelled backscatter.

Because of the loss of sensitivity of the backscatter for increasing GSV, the retrieval was constrained within a range of values between 0 and an estimate of the maximum GSV so that unnatural values were avoided [24]. An estimate larger than the maximum GSV was replaced with the maximum GSV value. Likely, a negative estimate of GSV was replaced with 0. In this study, we constrained the retrieval to a modelled value, V_max, using Equation (11). The first term represented the value obtained from the maximum canopy height, h_max, and Equation (5). As this value corresponded to an average value of several potential maximum GSVs, we introduced the second term that accounted for the range of maximum GSV estimates. This term was expressed with the standard deviation of the maximum GSV estimate corresponding to the maximum canopy height, $\delta V\left(h_{max}\right)$.

$$V_{max}=a·{\left(h_{max}\right)}^b+2\ast \delta V\left(h_{max}\right)$$

(11)

With several backscatter measurements from one sensor, the corresponding estimates of GSV were then combined to form a new estimate that averaged out individual noises. The weighted average in Equation (12) combined the individual GSV estimates from SAR data of a given sensor and defined the weights as the difference between the backscatter parameters of the WCM, w_i = σ⁰_veg − σ⁰_gr [23].

$${GSV}_{mt}={\displaystyle \frac{\sum _{i=1}^N{w}_i{GSV}_i}{\sum _{i=1}^N{w}_i}}$$

(12)

To alleviate systematic deficiencies in GSV estimates from multiple sensors, the sensor-specific estimates of GSV were eventually combined with an additional weighting scheme [43] where the definition of the weights accounted for (i) differences in the sensitivity of C- and L-band data to GSV and (ii) residual topographic effects in the ALOS-2 PALSAR-2 mosaics.

The first weighting scheme accounted for the sensitivity of the backscatter to GSV at a given frequency. It was modelled with the transmissivity term in Equation (1).

$$T=1-\eta +{\eta \sigma }_{gr}^0$$

(13)

A simple weighting scheme that reflected the difference in sensitivity between C- and L-bands consisted of differencing in the derivatives of Equation (13):

$$w_s=\left|{\displaystyle \frac{\partial T_L}{\partial V}}\right|-\left|{\displaystyle \frac{\partial T_c}{\partial V}}\right|$$

(14)

where the derivatives of the transmissivity at L- and C-bands, T_L and T_C, were determined using the structural functions in Equations (4) and (5). Equation (14) was then constrained with the number of images available at L- and C-bands, N_L and N_C, respectively:

$$w_s=\left|{\displaystyle \frac{\partial T_L}{\partial V}}\sqrt{N_L}\right|-\left|{\displaystyle \frac{\partial T_c}{\partial V}}\sqrt{N_C}\right|$$

(15)

This weight was finally rescaled to the range 0 to 1 to obtain normalized weights:

$$w_{s,norm}={\displaystyle \frac{w_s-min\left(w_s\right)}{max\left(w_s-min\left(w_s\right)\right)}}$$

(16)

The second weighting scheme was introduced because the effect of the pixel scattering area on the SAR backscatter was compensated for differently in the C- and L-band datasets. The additional weight was defined by calculating the percent difference in pixel scattering area between flat terrain (i.e., for the nominal incidence angle of 38° of ALOS-2 PALSAR-2) and the pixel area estimated from the local incidence angle, θ_i.

$${\Delta }_{area}=100·\left|1-{\displaystyle \frac{\mathit{sin}\left(38°\right)}{\mathit{sin}{\theta }_i}}\right|$$

(17)

The weight, w_topo, was then obtained by linearly scaling Δ_area between 0 and 1 with 1 representing flat terrain (i.e., for a 38° incidence angle) and 0 reflecting Δ_area values beyond 30, i.e., for steep slopes.

The final weight for the L-band was then the product of w_s and w_topo. Accordingly, the weight for the C-band was the complement of the L-band weight.

$$w\left(L\right)=w_{s,norm}·w_{topo}$$

(18)

$$w\left(C\right)=1-w\left(L\right)$$

(19)

The final estimate of GSV, GSV_fin, was obtained from the C- and L-band GSV estimates, GSV_mt,C and GSV_mt,L, with Equation (20):

$${GSV}_{fin}=\mathrm{w}\left(\mathrm{L}\right)·{GSV}_{mt,L}+\mathrm{w}\left(\mathrm{C}\right)·{GSV}_{mt,C}$$

(20)

4.4. Implementation of GSV Retrieval and Validation

To train the WCM, the field measurements of GSV were sorted and every other sample was selected. The remaining samples formed the test set. With this procedure, we ensured that both sets had the same statistical distribution of GSV. The test set was used to validate GSV estimated from both the trained and the calibrated WCM.

Equation (12) combined the GSV estimates from the Sentinel-1 images and from several combinations of the six ALOS-2 PALSAR-2 mosaics separately. After that, a final estimate of GSV was obtained for each of the two approaches using Equation (20). All retrievals were assessed against the field measurements of GSV in the test set. The retrieval error was quantified with

• The root mean square error (RMSE);
• The RMSE relative to the mean value of the GSV measurements in the test set;
• The bias, i.e., the systematic difference between estimated and reference GSV values;
• The coefficient of determination (R²).

5. Results

5.1. Correlation Between Radar Backscatter and GSV

Pearson’s coefficient of correlation between the plot-level GSV values and the backscatter observations was computed for each SAR image to understand the strength of the relationship between the two covariates (Figures S1 and S2). Statistics derived from the correlation coefficients are reported in

Table 5. Pearson’s coefficient of correlation between plot inventory GSV and Sentinel-1 backscatter. For each site and polarization, the mean, minimum, and maximum correlation values across the Sentinel-1 dataset are reported.

Site	VV			VH
	Mean	Min	Max	Mean	Min	Max
Finland N	0.50	0.38	0.73	0.52	0.42	0.65
Finland S	0.15	−0.13	0.40	0.14	−0.04	0.38
Catalonia	0.37	0.26	0.48	0.42	0.25	0.54
Romania	0.17	0.08	0.26	0.14	0.01	0.31

for the Sentinel-1 dataset and in

Table 6. Pearson’s coefficient of correlation between plot inventory GSV and ALOS-2 PALSAR-2 backscatter. For each site and polarization, the mean, minimum, and maximum correlation values across the ALOS-2 PALSAR-2 dataset are reported.

Site	HH			HV
	Mean	Min	Max	Mean	Min	Max
Finland N	0.31	0.24	0.42	0.51	0.44	0.62
Finland S	0.42	0.29	0.54	0.48	0.36	0.60
Catalonia	0.21	0.17	0.25	0.35	0.34	0.39
Romania	0.15	0.13	0.20	0.23	0.19	0.25

for the ALOS-2 PALSAR-2 dataset. Although this was not a rigorous way to assess the strength of the relationship because Pearson’s coefficient assumes collinearity between two variables while the backscatter at C- and L-bands loses sensitivity to GSV for increasing GSV, it provided indications on the theoretical performance of the GSV retrieval in space and time.

Except for images acquired in early spring over Finland N, the correlation between the Sentinel-1 backscatter and GSV was low (i.e., <0.3) to moderate (i.e., between 0.3 and 0.6) (Table 5 and Figure S1). While Finland N was characterized by the largest correlation, Finland S was characterized by the lowest correlation among all sites. This result was explained by the larger range of GSVs and the wet conditions prevalent at this site, as both factors negatively impact the correlation between biomass and SAR backscatter [19]. In addition, Finland S was the only site where different levels of correlation for ascending and descending passes (i.e., night-time and daytime) were observed (Figure S1). The night-time data were all characterized by correlation coefficients close to 0 whereas the correlation coefficients for the images acquired during daytime were higher. As for the two Finnish sites, the correlation coefficients in Catalonia were similar at VV and VH polarization (Table 5), with some seasonality occurring at VH polarization (Figure S1). On the contrary, in Romania, the correlation coefficient was low regardless of polarization because of the large proportion of plots with high GSV.

The correlation coefficients for the ALOS-2 PALSAR-2 dataset were either low or moderate (Table 6) with a few exceptions in Finland S (Figure S2). Compared to the C-band values, the correlation coefficients were higher at Finland S whereas they were lower in Finland N and Catalonia. The images used by JAXA to create the ALOS-2 PALSAR-2 mosaics over these sites were acquired in October. For the Finland N site, the weak correlation was a consequence of the air temperature being around 0 °C, this being consistent with previous studies [71]. In Catalonia, the weather was mostly dry and warm, thus not explaining the weaker correlation at the L-band compared to C-band observations. The visual analysis of the ALOS-2 PALSAR-2 mosaics revealed imperfections caused by sloped terrain, leading to strong variation in backscatter observations in relation to the GSV measurements. The low correlation between GSV and ALOS-2 PALSAR-2 backscatter for the Romanian site was explained by the overall high GSV.

5.2. Estimation of Parameters of GSV Structural Function

Figure 2 shows the regression of Equation (5) to the measurements of average canopy height from the ICESat-2 dataset and the corresponding GSV value for 377 NFI units in 26 countries. Data from Switzerland, Austria, Ireland, and Great Britain were discarded because for each of these countries, the correlation between canopy height and GSV values was smaller than 0.3. Topography in Austria and Switzerland and the fragmented forest landscape in all four countries except Austria likely contributed to observed weak correlation between forest variables.

5.3. Estimation of the Water Cloud Model Parameters

The coefficient α was obtained by allowing it to vary between a certain range of values and fitting each time the WCM to the observations of SAR backscatter and GSV. Based on experimental results reported in the literature, a plausible range of values between 0.1 dB and 2 dB was used [72,73]. For each value of α, the sum of the model residuals was computed for each SAR image at a given site. The value corresponding to the smallest of the cumulative residuals was selected as an estimate of α.

At the C-band, the procedure was strongly affected by the weak correlation between GSV and backscatter. The cumulative residuals showed no dependence on α except for the Finland N site where the sum of the residuals decreased as α increased. Since an estimate of 2 dB/m was considered plausible in previous studies [23], it was adopted regardless of the polarization, acquisition date, or site. At the L-band, we identified local minima in the cumulative residuals as a function of α, which were site-specific, being 0.8 dB/m in Finland S and Catalonia, and 1 dB/m in Finland N. Because of the very weak correlation between GSV and backscatter in Romania, the estimate of 0.5 dB/m was considered unreliable. As a first-order approximation, we finally assigned α the average of the site-specific values, i.e., 0.9 dB/m.

Figure 3 and Figure 4 show examples of measured and modelled backscatter as a function of canopy density for a pair of Sentinel-1 images covering the sites with the highest correlation between backscatter and GSV, i.e., Catalonia and Finland N. The relationship between the backscatter and canopy density differed depending on the local incidence angle. The sensitivity increased from steep to shallow angles as the thickness of the vegetation layer sensed by Sentinel-1 grew and the backscatter from the ground beneath the canopy diminished. Although the variability in the backscatter observations was substantial, the mean value per canopy density level presented a clear trend, which was captured by the WCM in Equation (8). Compared to the value representative of the backscatter for dense forests, the value of ${\widehat{\sigma }}_{veg}^0$ obtained from the regression of the WCM to the data was slightly higher.

To estimate the offset term ${SD}_{CD=100\%}$ in Equation (9), a model was fitted to canopy density and backscatter standard deviation per canopy density level. Figure 5 and Figure 6 show the relationship between canopy density and backscatter standard deviation for the same Sentinel-1 images used in the displays of Figure 3 and Figure 4, respectively. The standard deviation decreased monotonically with canopy density, and the use of a linear function was deemed reasonable.

For the ALOS-2 PALSAR-2 data, Figure 7 shows measured and modelled backscatter as a function of canopy density for the Finland N site where the correlation between GSV and SAR backscatter was the largest among all study sites. And Figure 8 shows the standard deviation of the backscatter modelled as a function of canopy density for the same site. The sensitivity of the backscatter to canopy density was strongest for the 30–40 deg incidence angle interval, regardless of polarization (2 dB at HH and 4 dB at HV polarization). The WCM reproduced the trend of the measurements. The correction applied to the estimate of ${\widehat{\sigma }}_{veg}^0$ was approximately 1 dB. This is a result of compensating for the standard deviation of the observations for a canopy density of 100%, which was indeed slightly larger than 0.5 dB (Figure 8), with the ENL-related term.

Figure 9 shows an example of trained and calibrated WCM for the same Sentinel-1 images over Catalonia and Finland N used in Figure 3 and Figure 4. Each panel also reports the mean value of the model residuals obtained with training and calibration. A comparison of all Sentinel-1 estimates of σ⁰_gr and σ⁰_veg from Catalonia and Finland N obtained with the two model-fitting procedures is illustrated in Figure 10.

The same analysis was undertaken for the ALOS-2 PALSAR-2 dataset over the study sites of Catalonia, Finland N, and Finland S. Figure 11 shows measured and modelled backscatter as a function of GSV for HH and HV polarizations for Catalonia, Finland N, and Finland S. The training and calibration of the WCM were tested on all ALOS-2 mosaic data and the estimates of σ⁰_gr and σ⁰_veg are displayed in Figure 12.

5.4. Retrieval of GSV

GSV retrieval was performed using the trained and the calibrated WCM in Equation (6). The retrievals were constrained by the maximum GSV values obtained with Equation (11) in which we applied site-specific values of (i) ICESat-2 maximum canopy height (Table 4), (ii) the structural model coefficients a and b in Figure 2, and (iii) the standard deviation of GSV values corresponding to the ICESat-2 h_max value.

Table 7. Site-specific maximum GSV derived from the plot-based measurements and modelled with Equation (11) using the ICESat-2 maximum canopy height in Table 4. Values in brackets represent the percentile at which the plot-based GSV measurements are equal to the modelled value of maximum GSV.

Site	Maximum GSV [m3/ha]
	Plot-Based	Modelled
Finland N	498	238 (96th)
Finland S	751	395 (94th)
Catalonia	480	504 (100th)
Romania	1677	1295 (99th)

compares our estimates of maximum GSV with the largest value from the dataset of plot-based measurements.

The performance of the retrieval using the Sentinel-1 dataset is illustrated with scatter plots comparing reference and estimated GSV and with retrieval statistics for the two sites of Catalonia and Finland N in Figure 13. For the remaining sites where the correlation between the Sentinel-1 backscatter and GSV was low, the scatter plots are included in the (Supplementary Material Figure S3).

The benefit of using multiple observations at the L-band is illustrated in Figure 14, which shows the retrieval results for the ALOS-2 PALSAR-2 mosaics of the same year of the field inventory data (2 SAR images) and for the 2015–2021 interval (12 SAR images). Results are shown for the Finland S site characterized by the highest correlation coefficients between SAR backscatter observations and GSV. The accuracy of the retrieval increased with the number of SAR observations; nonetheless, combining data acquired over several years might have been detrimental because the risk of not accounting for changes on the ground due to growth and losses would have biassed the GSV estimates. In the remainder of this section, we therefore assessed the L-band retrieval with ALOS-2 PALSAR-2 mosaics based on images acquired between the year prior to and after the field inventory. Results are displayed in Figure 15 for the sites of Catalonia, Finland N, and Finland S. The results for the Romanian site where correlation between SAR backscatter and GSV was low are illustrated in the (Supplementary Material Figure S4).

The weights used for merging the Sentinel-1 and PALSAR-2 GSV estimates differed depending on the site. The weighting factor for the L-band GSV estimates ranged between nearly 0 (Finland N), 0.23 (Catalonia), 0.54 (Romania), and 0.9 (Finland S). The scatter plots illustrating reference and estimated GSV are illustrated for each site in Figure 16.

6. Discussion

Estimating forest GSV using satellite data that correlate moderately or poorly with biomass suggests that any imperfection in the setup of the retrieval approach would cause substantial estimation errors. In the retrieval scheme proposed here, errors could be due to (i) the structural functions relating the forest variables in the WCM with GSV, (ii) the data used for training the structural models, (iii) the forest backscatter model, (iv) the procedure for estimating the parameters of this model, and (v) the procedure for combining GSV estimates from individual observations of the SAR backscatter. Assessment of any of these aspects requires a dataset of reference GSV measurements. In this study, such data consisted of plot-level measurements collected by field inventories at four sites. At the spatial scale of individual pixels, our datasets of SAR backscatter and GSV measurements were moderately or poorly correlated (Table 5 and Table 6). Plot-level data and pixel-level backscatter are known to correlate weakly. This effect was reinforced by the limited penetration of the C- and L-band signal into the canopy (see, e.g., Figure 9 and Figure 11) and, in part, by the pre-processing of the SAR data, which especially over sloped terrain accounted for the pixel area effect on the backscatter but not for scatterer-specific effects [74]. Under these circumstances, we acknowledge that the setup for assessing the performance of the method proposed here is sub-optimal. Nonetheless, our interest was to explain errors associated with a model calibration approach independent of GSV reference data. Hence, the plot-based datasets of GSV were assumed to be the most accurate data source that could be used for model calibration; i.e., the errors associated with forest inventories were not further addressed.

The training and calibration of the WCM relied on the same set of structural functions. Although the plot data could have been used to derive site-specific functions, our focus was on the WCM calibration so that the same structural functions had to be used to maintain consistency in the retrieval phase. While the relationship between canopy density and canopy height in Equation (4) was extensively documented, the experimental evidence for the power-law function relating canopy height and GSV was relatively limited. To overcome the lack of GSV measurements for large areas, we trained Equation (5) with sub-national average values that provided an almost gap-free coverage. The structural function in Equation (5) reproduced the relationship between canopy height and GSV in each of the four major European biomes. The two types of temperate forest biomes and the Mediterranean forest biome were characterized by a power-law relationship resembling a quadratic function (i.e., exponent close to 2). Boreal forests exhibited a more linear relationship with an exponent close to 1. The standard error (SE) of the regression indicated a low error across all biomes (between 9% and 20% of the average GSV) except for Mediterranean forest. There, SE of 45% relative to the average GSV was caused by the large spread of the observations, which was due to the variety of vegetation types in this environment.

The validity of the curves in Figure 2 can be questioned when considering that the structural function was trained on large-scale averages while our aim was to apply it at the pixel level of a SAR image. Our understanding was that the region of validity of the curves corresponded to the range of GSV values present in the training dataset (e.g., up to 400 m³/ha in temperate broadleaf and mixed forests). The structural function might instead produce inaccurate predictions of GSV for forest patches with GSV exceeding the largest average GSV in the dataset. As shown by the GSV statistics in Table 2, the largest of the GSV values in the dataset of sub-national statistics was considerably smaller than the maximum value of GSV at the plot level and in some cases comparable with the third quartile (Finland S) or with the second quartile (Romania). The impact of a presumably imperfect calibration of the structural function was further evidenced by the different values for the maximum GSV (Table 7). The estimate of the maximum GSV with the structural function for Catalonia was in line with the largest values reported in the plot-level GSV dataset. For both Finnish sites and the Romanian site, the maximum GSV predicted by the structural function instead corresponded only to a very high percentile in the reference plot dataset, being much smaller than the maximum measured GSV. For the two sites in the boreal zone, the discrepancy was attributed to the rather linear trend between height and GSV (Figure 2), which was caused by the lack of data points with a GSV above 200 m³/ha so that the high end of the function, where typically a non-linear correspondence occurs, could not be accurately constrained. Indeed, when fitting Equation (5) to the plot data, we obtained an estimate for the exponent close to 2 (not shown here), leading to a much steeper increase in GSV with height. For the Romanian site, the large difference between maximum GSV values (>400 m³/ha) was apparent since only 12 of the total 1306 plots had GSV larger than our model’s maximum estimate (Figure 2).

The SAR backscatter steadily increased with canopy density, regardless of the site, frequency, polarization, or incidence angle. The strongest sensitivity was observed at the L-band for HV-polarized data with an increase of 2.5 dB up to 4 dB (Figure 7). C-band VH-polarized data showed an increase of 2–3 dB (Figure 3 and Figure 4). The co-polarized data were characterized by weak sensitivity to canopy density, being mostly between 1 dB and 2 dB (Figure 3, Figure 4 and Figure 7). The overall poor-to-moderate association between SAR backscatter and canopy density did not affect the performance of the calibration. Indeed, the WCM in Equation (2) fitted to the observations of canopy density and backscatter was able to reproduce their trend (Figure 3 and Figure 4 at C-band and Figure 7 at L-band) and generated plausible estimates of σ⁰_gr and σ⁰_veg despite the large range of backscatter values for a given canopy density. The low ENL of the Sentinel-1 and ALOS-2 PALSAR-2 datasets (Table 3) partly explaining the large spread of the observations caused the strong correction to be applied to the initial estimate of σ⁰_veg to obtain its final value with Equation (9).

The relationship between canopy density and variability of the SAR backscatter measurements after correcting for their uncertainty (Figure 5 and Figure 6 at C-band and Figure 8 at L-band) was clearly linear, so that the correction factor to be applied to the initial estimate of σ⁰_veg was easily determined. This correction factor was between 0.2 and 0.5 dB at the C-band (Figure 5 and Figure 6) and around 0.5 dB at the L-band (Figure 8). For any SAR image, the resulting estimate of σ⁰_veg was among the largest measured SAR backscatter values within the study site, this being consistent with our assumption that σ⁰_veg is an ideal value of the backscatter for the largest possible stocks of wood volume in the area of interest.

The calibrated WCM was close to the trained WCM and both curves represented the trend in the observations (Figure 9 for C-band and Figure 11 for L-band). The mean of the model residuals confirmed the visual impression that the model fit obtained with the calibration tended to slightly overestimate the backscatter, regardless of frequency and polarization. The bias of the calibrated models was, however, substantially smaller than the range of measured backscatter values and less than the uncertainty characterizing the SAR backscatter data. In addition, the moderate correlation between backscatter and GSV did not allow us to conclude to which extent one curve is more representative of the true association between the two variables. Indeed, the trained WCM in Catalonia was characterized by somewhat unrealistic values for σ⁰_gr both at C- (Figure 9) and L-bands (Figure 11), because they were lower than any measurement, and by a steep increase in the modelled backscatter in correspondence with a moderate increment of GSV.

The comparison of all estimates of σ⁰_gr and σ⁰_veg from the training and the calibration showed strong agreement at both C- (Figure 10) and L-bands (Figure 12) except for the Catalonian site where a systematic offset for σ⁰_gr at the C-band for both polarizations and at the L-band at cross-polarization occurred. The calibration generated larger values than those obtained with the training of the WCM (Figure 11 and Figure 12), the latter one being unrealistically low.

With Sentinel-1 C-band data, the retrieval based on the calibrated WCM was comparable to the retrieval based on the trained WCM for the two sites of Finland N and Catalonia characterized by the highest correlation between backscatter and GSV (Figure 13). The retrieval error of 56% and 75% relative to the mean GSV was only minimally affected by the choice of the model-fitting procedure. The discrepancy between the retrievals for the two other sites was instead substantial and none of the model-fitting approaches could reproduce the plot-based GSVs in the test set. The overall large retrieval error was explained by the considerable variance of model predictions (Figure 13), a consequence of the weak sensitivity of the C-band backscatter to GSV (Figure 9). The accuracy was similar to previous studies comparing the estimates of forest biomass (GSV or AGB) at pixels with plot-level values. Nonetheless, the level of GSV was estimated correctly, with a tendency to underestimate the largest GSVs. This was a consequence of (i) the rather large estimate of σ⁰_veg and (ii) the inversion rule that set a maximum retrievable value, which often was smaller than the site-level value derived from the plot data.

With ALOS-2 PALSAR-2 data, the retrieval based on the trained WCM produced estimates of GSV that were only slightly overestimated in the central range of GSVs and underestimated for the highest GSVs (Figure 15). As in the C-band case, these were underestimated because of the limit imposed to the range of retrievable GSVs and because of the high value of σ⁰_veg. The performance of the retrieval from the calibrated WCM was comparable to the best-case scenario at the two Finnish sites (Figure 15). For Catalonia, slight underestimation occurred because of the slightly higher estimates of σ⁰_veg (Figure 12). This could be due to imperfect terrain flattening of the mosaic leading to uncompensated backscatter in sloped terrain. Translating variance of backscatter into variance of GSV estimations explained the large retrieval error between 63% and 87%. The error was marginally affected by the procedure used to estimate the model parameters, except for Catalonia. At the Romanian site, the retrievals were characterized by poor accuracy (Figure S4) due to the low correlation between forest backscatter and GSV (Table 6).

The weights used to merge the C- and L-band GSV estimates differed depending on the study site. For Catalonia, the weight was geared towards the C-band estimates, which was reasonable given the poorer performance of the retrieval with the L-band (Figure 15). The merged GSV estimates (Figure 16) were slightly more accurate than in the C-band-only case (Figure 13) with the RMSE decreasing from 79.1 m³/ha to 71.7 m³/ha and the bias decreasing from 4.6 m³/ha to 0.7 m³/ha for the calibrated WCM. For Finland N, the weighting excluded the L-band estimates although their accuracy was comparable to the C-band retrieval (Figure 13 and Figure 15). This was a shortcoming of how the normalization procedure in Equation (16) was implemented for this analysis. For Finland S, the L-band estimates were preferred to the C-band estimates, being consistent with the superior performance of the L-band retrieval and a consequence of the higher correlation between backscatter and GSV (Table 6). The error statistics of the merged estimates (Figure 16) were comparable to those obtained for the L-band estimates and were only minimally affected by the C-band retrieval. For Romania, the weights were similar, as a consequence of the poor correlation between GSV and backscatter at both C- and L-bands (Table 5 and Table 6) so that none of the GSV estimates prevailed.

The accuracy of the merged estimates of GSV from the trained and from the calibrated WCM was comparable (Figure 16). The closest agreement was obtained in Catalonia with a relative RMSE close to 65%, R² of approximately 0.7, and negligible bias for both model-fitting approaches. At Finland N, the retrieval based on the calibrated WCM was more accurate below 100 m³/ha, whereas the opposite occurred above 100 m³/ha with increasing underestimation regardless of the model-fitting approach. At Finland S, both sets of GSV estimates were characterized by both over- and underestimations in an equal manner. We explained this result as partial incapability of the WCM to reproduce the relationship between SAR backscatter and GSV. In Romania, the retrieval led to strong overestimation, which was expected given the lack of correlation between SAR backscatter and GSV. An investigation to explain the origin of the errors was beyond the scope of this study.

When putting these overall encouraging results in a broader context of biomass estimation with SAR backscatter data, our results show that the levels of GSV are well captured with the inversion of a Water Cloud Model and the combination of estimates from C- and L-bands. Except for the Romanian site, the estimates of GSV clustered along the identity line with systematic underprediction occurring for the highest GSVs only. Nonetheless, the estimates of GSV are characterized by substantial variance as shown by relative RMSE values between 53% and 79% of the average GSV and R² values between 0.61 and 0.78 (Figure 16). The weak sensitivity of the C- and L-band backscatter to GSV explains the large dispersion of the estimates.

Unlike studies that explained the systematic underprediction of high biomass as a consequence of backscatter saturation [75], we instead associate biases to either an incorrect value of the maximum retrievable GSV (see also [43,44]) or an overestimate of the parameter σ⁰_veg (see also [44]). As C- and L-band backscatter lose sensitivity for increasing AGB or GSV, measurements that are outside of the range of values predicted by a retrieval model might either not be inverted or associated with a biomass value within the range of the modelled backscatter. As a consequence, the high range of biomass is underestimated. By instead applying a constraint such as a maximum retrievable biomass, such backscatter measurements are here associated with a plausible biomass value (see, e.g., Figure 14). Because of the mostly weak sensitivity of the C- and L-band backscatter to biomass, the biomass estimates from a multi-temporal stack of observation are noisy and the combination of values suggested in this study acts as a filter so that (i) the impact of erroneous values and “hard coded” biomass on the final estimates is limited, and (ii) the statistical distribution of biomass is well reproduced.

7. Conclusions

This study compared the performance of forest GSV retrieval using two approaches to fit C- and L-band SAR backscatter data to forest GSV with a physically based WCM. Our key finding was that the performance of the WCM calibrated on auxiliary image data and independent from reference GSV data was comparable to a WCM trained on such reference data. To corroborate the results, the analysis was undertaken at four sites representing various forest conditions in Europe. Although the GSV estimates from the SAR data and from the plot inventory measurements were affected by large spread, the proposed procedure was able to correctly reproduce the spatial distribution of GSV except for very high GSV over sloped terrain. A calibrated WCM is therefore a viable approach when GSV reference data are poor or unavailable. This scenario applies to wide-area mapping, e.g., at the continental level, when trans-national reference measurements may be cumbersome to access due to data policy restrictions or difficult to interpret due to different measuring protocols and reporting.