Predictions of Biomass Change in a Hemi-Boreal Forest Based on Multi-Polarization L- and P-Band SAR Backscatter

Abstract Above-ground biomass change accumulated during four growth seasons in a hemi-boreal forest was predicted using airborne L- and P-band synthetic aperture radar (SAR) backscatter. The radar data were collected in the BioSAR 2007 and BioSAR 2010 campaigns over the Remningstorp test site in southern Sweden. Regression models for biomass change were developed from biomass maps created using airborne LiDAR data and field measurements. To facilitate training and prediction on image pairs acquired at different dates, a backscatter offset correction method for L-band data was developed and evaluated. The correction, based on the HV/VV backscatter ratio, facilitated predictions across image pairs almost identical to those obtained using data from the same image pair for both training and prediction. For P-band, previous positive results using an offset correction based on the HH/VV ratio were validated. The best L-band model achieved a root mean square error (RMSE) of 21 t/ha, and the best P-band model achieved an RMSE of 19 t/ha. Those accuracies are similar to that of the LiDAR-based biomass change of 18 t/ha. The limitation of using LiDAR-based data for training was considered. The findings demonstrate potential for improved biomass change predictions from L-band backscatter despite varying environmental conditions and calibration uncertainties.

To facilitate training and prediction on image pairs acquired at different dates, a backscatter offset correction method for L-band data was developed and evaluated. The correction, based on the HV/VV backscatter ratio, facilitated predictions across image pairs almost identical to those obtained using data from the same image pair for both training and prediction. For P-band, previous positive results using an offset correction based on the HH/VV ratio were validated. The best L-band model achieved a root mean square error (RMSE) of 21 t/ha, and the best P-band model achieved an RMSE of 19 t/ha. Those accuracies are similar to that of the LiDAR-based biomass change of 18 t/ha. The limitation of using LiDARbased data for training was considered. The findings demonstrate potential for improved biomass change predictions from L-band backscatter despite varying environmental conditions and calibration uncertainties. RÉSUMÉ Les changements de la biomasse accumul ee au-dessus du sol dans une forêt semi-bor eale ont et e pr edits au cours de quatre saisons de croissance a l'aide de la r etrodiffusion d'un radar a synth ese d'ouverture (RSO) a eroport ee en bande L et P. Les donn ees radar ont et e recueillies durant les missions BioSAR 2007 et BioSAR 2010 sur le site test de Remningstorp, dans le sud de la Su ede. Des mod eles de r egression pour le changement de biomasse ont et e d evelopp es a partir de cartes de biomasse cr e ees a l'aide de donn ees LiDAR a eroport ees et de mesures sur le terrain. Afin de faciliter l'entraînement et la pr ediction a partir de paires d'images acquises a diff erentes dates, une m ethode de correction du d ecalage de la r etrodiffusion pour les donn ees en bande L a et e elabor ee et evalu ee. La correction, bas ee sur le rapport de r etrodiffusion HV/VV, a facilit e l'obtention pour des paires d'images de pr edictions similaires a celles obtenues pour la même paire d'images et ce tant pour l'entraînement que la pr ediction. Pour la bande P, des r esultats ant erieurs positifs, utilisant une correction de d ecalage bas ee sur le rapport HH/VV, ont et e valid es. Une erreur moyenne de la carr ee racine (RMSE) de 21 t/ha et de 19 t/ha a et e respectivement obtenue pour le meilleur mod ele en bande L et pour celui en bande P. Ces pr ecisions du changement de biomasse etaient semblables a celle des donn ees LiDAR (18 t/ha). La limitation de l'utilisation des donn ees LiDAR pour l'entraînement a aussi et e prise en consid eration. Les r esultats d emontrent un potentiel d'am elioration des pr evisions du changement de la biomasse a partir de la r etrodiffusion de la bande L, malgr e des conditions environnementales variables et les incertitudes d' etalonnage.

Introduction
Increased concentrations of greenhouse gases in the atmosphere cause global warming, and the largest contributor to this process is carbon dioxide (Stocker et al. 2013). The biosphere acts as a net carbon sink, and changes in its efficiency in storing carbon need to be better understood (Canadell et al. 2007). Forests constitute a large part of the terrestrial biosphere and their carbon storage is proportional to the biomass they contain. This calls for large-scale mapping of forest biomass. Mapping of forest biomass is also of interest in assessing fire, storm, insect, and disease damages, and to support decisions in commercial forest management. Synthetic aperture radar (SAR) systems have the benefit of being able to form images of the terrestrial landscape regardless of clouds, precipitation, and sun illumination conditions, which can inhibit LiDAR or optical systems, making it easier to reliably get full coverage data even over large areas. The relatively long wavelength of L-and P-band SAR systems compared to sensors operating with shorter wavelengths or in the optical domain is an advantage because of the signal penetration into the canopy, which facilitates estimation of above-ground biomass (AGB) (Kasischke et al. 1997).
Advanced SAR methods like InSAR (Solberg et al. 2014;Persson and Fransson 2017) or tomography (Ho Tong Minh et al. 2016) have been used to predict AGB or AGB change, but place high demands on the SAR data in terms of temporal coherence that not all SAR products can meet. In particular, L-band data often do not possess the temporal coherence needed for repeat-pass interferometry over forests (Hamadi et al. 2017;Monteith and Ulander 2018). While data from formation flying satellites ought to overcome this limitation, no such L-band systems are currently in operation. Backscatter-based algorithms for AGB or AGB change prediction do not have this inherent coherency requirement, and can be applied to both Land P-band data. In this context, simple backscatterbased predictions are more feasible for large-scale mapping of carbon stock dynamics on a national or continental level. Additionally, backscatter algorithms are simpler to implement for non SAR experts, since they can be applied to standard SAR products without the use of specialized SAR software.
Previous studies have used both L-and P-band SAR data to predict AGB, or the highly correlated variable stem volume, in boreal and hemi-boreal forests (Rignot et al. 1994;Kurvonen et al. 1999;Saatchi and Moghaddam 2000;Rauste 2005;Sandberg et al. 2011;Neumann et al. 2012;Santoro et al. 2015; Schlund and Davidson 2018;Cartus et al. 2019;Santoro et al. 2019). L-band data have also been used to detect clear-cuts and storm damage (Fransson et al. 2007;Santoro et al. 2012), and to estimate AGB growth from multi-satellite SAR data (Balzter et al. 2003). However, algorithms based on backscatter are known to be subject to saturation (loss of sensitivity) when biomass increases beyond a certain level. The onset of sensitivity loss is dependent on the wavelength and forest type, and occurs at 40-100 t/ha and 100-200 t/ha for L-and P-band, respectively, in conifer-dominated forests (Dobson et al. 1992;Le Toan et al. 1992;Rignot et al. 1994;Imhoff 1995;Fransson 1999). While L-and P-band data have previously been compared in their ability to estimate or predict AGB, this study evaluated both bands for prediction of AGB change. The data used in this study are highly comparable between bands as it was acquired in the same geometry over the same area, with almost all acquisitions made on the same days for both bands.
This study extends the work of Sandberg et al. (2014), who analyzed the P-band data from the BioSAR 2007 and BioSAR 2010 campaigns, but omitted the L-band data from the analysis. While the coming BIOMASS mission will operate in P-band, in light of both the scarcity of available P-band data over Europe, and North and Central America in the present and foreseeable future, following the restrictions imposed by the US Department of Defense (Carreiras et al. 2017), and the present and planned L-band missions including ALOS-2, ALOS-4, NISAR, ROSE-L, and SAOCOM (The CEOS database: Missions, instruments, measurements and datasets 2020), it was of interest to expand the analysis to include the L-band data from the same campaigns. Currently, only ALOS-2 PALSAR-2 can, to a certain extent, provide continual global L-band backscatter data sets. In Huuva et al. (2017), single polarization AGB change models were investigated for the L-and P-band data from the campaigns, and it was found that using models based on HV backscatter explained most of the variation in AGB change for both L-and P-band.
In Sandberg et al. (2014), LiDAR-based AGB maps were used to train P-band SAR models. These AGB maps were, however, created using a different method for the two acquisition years. For 2007, the AGB map was created in two steps, by first creating a stem volume map based on field and LiDAR data, which was then converted to an AGB map using biomass expansion factors. For 2010, the AGB map was instead created directly from AGB estimates from field and LiDAR data. In contrast, this study used LiDAR-based AGB maps created according to a unified and consistent method for both acquisition years to support the analysis and modeling. Unlike Sandberg et al. (2014), individual tree biomass models were applied to the field data for both years, eliminating the need to use biomass expansion factors when creating the AGB maps. This is expected to result in more accurate AGB values (Petersson et al. 2012). In contrast to the analysis in Sandberg et al. (2014), change predictions were back-transformed to facilitate the comparison of models based on different AGB change transformations. Additionally, while in that study the largest field measured changes used for evaluation were limited in magnitude to 50 t/ha, the evaluation data in the present study were expanded by including additional evaluation plots which were clear-cut between the acquisition years, thereby permitting evaluation of the models on much larger AGB changes, reaching 300 t/ ha in magnitude, which is interesting in light of the previously discussed saturation effects.
The intensity of SAR backscatter from a certain location is not only dependent on the amount of biomass and SAR parameters like wavelength, incidence angle, and polarization, but varies with other variables such as soil, stem, and canopy moisture, surface roughness and topography, and forest structure variables such as tree species, stem density, and spatial distribution of trees (Lucas et al. 2010;Kasischke et al. 2011). Moisture can vary significantly with not only seasons but even time of day due to precipitation and diurnal cycles, making it especially challenging when predicting AGB change. The non-AGB related change in backscatter may vary within a given SAR image pair due to, e.g. local precipitation in a region of one of the two images, resulting in wetter soils compared to other areas. Ideally, a non-constant offset correction across scenes would therefore be required. However, a first order approximation can be achieved by estimating the non-AGB related change with a constant offset for a given image pair. While such an offset in backscatter can, when modeling AGB changes from one SAR image pair, be corrected for by a constant term in the model, the offset will generally not be the same for another image pair, thus hindering the applicability of a trained model to other image pairs. To mitigate this effect of moisture and radiometric calibration uncertainty a backscatter change offset correction method was proposed for P-band data by Sandberg et al. (2014). The correction is based on the HH/VV backscatter ratio, found in P-band data to be relatively insensitive to moisture (Soja et al. 2013), and aims to correct the backscatter so that areas with small changes in AGB also have small changes in backscatter. While the correction method was successful in facilitating AGB change prediction from P-band backscatter data, no such correction method has been developed for L-band data.
The objectives of this study are to develop and evaluate a backscatter offset correction method for Lband SAR data, and to validate the existing correction method for P-band data. In addition, we compare the potential of predicting AGB changes in a hemi-boreal forest using L-and P-band backscatter, given that a suitable offset correction method can be applied to both bands.

Test site
The remote sensing data consisted of SAR and LiDAR data collected during the airborne SAR campaigns BioSAR 2007 and BioSAR 2010 conducted at the Remningstorp estate (Figure 1), in southern Sweden (58 30 0 N, 13 40 0 E). In addition, parts of the in situ data were also collected within the campaigns. Full descriptions of the campaigns can be found in the BioSAR reports by Hajnsek et al. (2008) and Ulander et al. (2011). The test site is located within the hemiboreal zone, which covers the transition between the boreal and temperate zones. This zone is characterized by mixtures of coniferous and deciduous species. The test site consists of managed forest dominated by Norway spruce (Picea abies (L.) H.Karst.) and Scots pine (Pinus sylvestris L.) with some birch (Betula spp.). In a Swedish context, the climate and geology are favorable for forest growth. Two thirds of the forest within the estate are located on till (i.e. a mixture of glacial debris) with a field layer consisting of different herbs, blueberry (Vaccinium myrtillus L.), and narrow-leaf grass (e.g. Deschampsia flexuosa (L.) Trin.). In denser old spruce stands the field layer is absent. The remaining forest grows on peatland, which is dominated by Scots pine. Regardless of tree cover, the main soil type on till is brown earth. Extensive ditching of peatlands on the estate over the years has been conducted to increase the productive forest area. However, some of the peatlands are frequently saturated (Ahlberg and Kardell 1997). The above-ground biomass range is 0-400 t/ha, at plot level (with 10 m radius). The ground slopes on the test site are generally small, with an average slope of 4 and a 95th percentile of 11.5 , as calculated from the 2 m Â 2 m resolution LiDAR-based digital terrain model (DTM). The ground surface topology lies between 120 m and 145 m above sea level.

In situ data
Three sets of field plot data were inventoried at the Remningstorp test site. The first set consisted of a systematic grid with 849 circular plots with 10 m radius and 40 m plot spacing. The inventory was started before, and completed after the vegetation season in 2005. In order to match the SAR data, the forest variable estimates were forecast to 2007, using established models in the Heureka forestry decision support system (L€ amås and Eriksson 2003;Wikstr€ om et al. 2011). The second data set, part of the BioSAR 2010 campaign, was inventoried after the growth season 2010, and consists of 216 circular plots with 10 m radius and 200 m spacing Ulander et al. 2011). The third data set consisted of ten 80 m Â 80 m plots, distributed across the Remningstorp estate, all but one fully contained within homogeneous stands. All ten large plots were surveyed just before the growth season of 2007. Three of the large plots were clear-cut between 2007 and 2010, and the remaining seven large plots were inventoried again after the growth season of 2010, thereby matching the SAR data without forecasting. The AGB was set to 0 t/ha in 2010 for the three large plots that were clear-cut between the acquisitions. One third of plot 9 was not clear-cut, since two of its corners are located in another, retained, stand. Therefore, only data from the central part of this plot, contained in the clear-cut, was used to represent the whole plot, as if the whole plot were clear-cut. Details of these large plots are given in Table 1. For all data sets, all trees with a diameter at breast height (DBH, measured at 1.3 m above ground) of more than 4 cm were calipered, species determined, and positioned using a real time kinematic global positioning system. For the data sets consisting of circular plots, the height was measured on about 10% of the trees, randomly sampled with probability proportional to the basal area, using a hypsometer. For the 80 m Â 80 m plots, the height of every calipered tree was measured. AGB was estimated using the Heureka system (Wikstr€ om et al. 2011), which used established models for height by S€ oderberg (1986) and AGB by Marklund (1988). Only plots found to be completely contained in a single stand, and with a Lorey's mean tree height (i.e. basal area weighted mean tree height) above 3.7 m (corresponding to DBH greater than 4 cm) were used in the study. This resulted in sample sizes of 763 and 208 plots from 2005 and 2010, respectively. In addition to the three field data sets, a stand delineation map over the area was used to allocate training and validation sets in the evaluation of regression models. Moreover, tree species information from a forest management plan of the estate was used in creating a LiDAR-based biomass map of the test site, as explained below. The average size of a stand in the test site was 2.5 ha, and the area covered by all remote sensing data sets contained 158 stands.
Since the 10 m radius field plot grids for the two years were not the same, and the number of shared plots was small, changes in AGB could not be directly measured in situ with a reasonable amount of field plots. Instead, AGB maps were created for 2007 and 2010 using the field plots and LiDAR data. The AGB maps are further described below. The 80 m Â 80 m plots were set aside for validation of the results, and were thus not used in model selection or training.

SAR data
The SAR scenes used are fully polarimetric in L-and P-band. For 2007, the SAR data were collected using the German E-SAR system by the German Aerospace Center (DLR), while the 2010 SAR data were collected using the French SETHI system by ONERA. The two systems and hence the SAR data for the two years are not identical. While the acquisitions in 2010 were all conducted on the same day, the SAR data for 2007 were collected during four dates ranging from early spring with recently thawed snow cover and wet soil, to late spring and dry soil. Soil moisture was continuously measured using an Aquaflex sensor installed 15 cm below ground in a forested area within the test site. In addition, air temperature and precipitation were recorded at a weather station located in the eastern part of the estate (Table 2). At the time of SAR acquisition, no precipitation was observed (Hajnsek et al. 2008;Ulander et al. 2011).
To eliminate the influence of different imaging geometries between the bands, only acquisitions with heading (along-track direction) 200 were chosen for the analysis. All of the selected acquisitions were made from the same flight track, i.e. the incidence angle for a given location on ground is the same in all images (regardless of sensor). The acquisition dates are given in Table 2. For both bands, one image per acquisition date in 2007 was included. The first and last 2007 acquisition dates are the same for L-and Pband, with about two months passing between them, while the intermediate date differs by two days between the bands. For each band, both scenes in the chosen heading from 2010 were included in the analysis to allow both images to be different across image pairs. However, the differences between the scenes from the same date are expected to be small. Table 3 lists the image pairs used in the subsequent analysis. The 2007 data, from E-SAR, were geometrically calibrated using corner reflectors deployed at the test site.  The reflectors were also used to validate the radiometric calibration of the data. For the 2010 data, from SETHI, the corner reflectors were used for both geometric and radiometric calibration. The radiometric calibration uncertainty was estimated to be ±1 dB for both SAR systems. The delivered geocoded standard products of the two SAR systems were used for the analysis. A full description of the SAR data can be found in the BioSAR reports (Hajnsek et al. 2008;Ulander et al. 2011). A brief summary of the main characteristics of the data used is presented here.
The 2007 SAR data, from E-SAR, were acquired with a center frequency of 1,300 MHz and a bandwidth of 94 MHz for L-band, and a center frequency of 350 MHz and a bandwidth of 70 MHz for P-band. The center frequencies correspond to wavelengths of 23 cm and 86 cm for L-and P-band, respectively. The single look resolutions of the data were 2.1 m Â 0.9 m for L-band and 2.1 m Â 1.6 m for P-band (slant range Â azimuth). The L-band data were multilooked by DLR 8 times, and the P-band data 4 times, both only in azimuth, with a 50% overlap. The resulting resolution for both bands was 2.1 m Â 4.0 m.
The 2010 SAR data, from SETHI, were collected with a center frequency of 1,325 MHz and a bandwidth of 150 MHz for L-band, and a center frequency of 360 MHz and a bandwidth of 166 MHz for P-band, with some notches and gaps in the latter spectrum to avoid radio frequency interference. The center frequencies correspond to wavelengths of 23 cm and 83 cm for L-and P-band, respectively. The single look resolutions of the L-and P-band data were 0.9 m Â 0.9 m and 0.8 m Â 0.8 m (slant range Â azimuth), respectively.
The SAR images were radiometrically corrected for local ground slope and incidence angle induced variations in the intensity and transformed to c 0 using a 2 m Â 2 m resolution DTM. From the DTM, the surface normal for each 2 m Â 2 m pixel was obtained. The angle w between the surface normal and the slant range image plane normal was then used to obtain the average radar cross section per unit ground area, r 0 , according to where b 0 is the average radar cross section per image pixel area (Ulander 1996). After this, a first order correction for variations in incidence angle was applied to obtain where h is the local incidence angle (Ulaby et al. 1982). The radiometrically corrected backscatter c 0 will henceforth be referred to as backscatter. The resulting backscatter images were averaged with a 50 m Â 50 m filter and resampled to backscatter maps with 50 m Â 50 m pixels following the methodology in Sandberg et al. (2014). The resampling has the effect of suppressing most of the speckle in the SAR data, as a pixel of 50 m Â 50 m contains between 700 and 4,000 resolution cells. This is also expected to make the differences in original resolution between the sensors negligible. The selected pixel size is on the same order of magnitude as the 80 m Â 80 m plots used for evaluation of AGB change predictions, while still permitting the models to capture biomass variations on a sub-stand scale.

LiDAR data
LiDAR acquisitions for both years were made using helicopter-mounted TopEye systems (Mk II in 2007, Mk III in 2010. The 2007 scanning was performed on April 24 with densities of 30-40 pulses/m 2 . In 2010, the scanning took place on August 29 with an average density of 69 pulses/m 2 . Further details can be found in the BioSAR reports (Hajnsek et al. 2008;Ulander et al. 2011). In processing the point clouds, point densities were locally found to be significantly lower than the averages, and to avoid artifacts in the All images were acquired from the same flight track, in the same imaging geometry.
AGB modeling from the varying point density, the point clouds were thinned to 10 points/m 2 to obtain an even point density throughout the test site.

LiDAR-based biomass maps
The field data from the two inventory years were not from the same set of plots, necessitating an indirect approach to estimate the AGB change from 2007 to 2010. Therefore, AGB maps were created using data from the LiDAR measurements. Multiple linear regression models were developed to relate in situ AGB values with LiDAR percentiles and density metrics. The biomass modeling included investigation of numerous transformations of both the response and explanatory variables, and inspection of scatter plots, residual plots, and covariances. As the AGB state predictions for each year were to be used to compute the AGB change, the determination of a common model that could describe both of the data sets well was prioritized. A log-log model was selected that fulfilled the requirements of linearity, normality and homoscedasticity for both LiDAR data sets. The relationship between AGB and LiDAR metrics was similar for spruce dominated plots for both years, while pine and deciduous dominated plots exhibited relationships that were different both from each other and from spruce dominated plots, and differed for the different years. The regression model used indicator and interaction variables to handle the different relations due to tree species (pine and deciduous) and acquisition year. A plot was defined as pine or deciduous dominated if the fraction of biomass from pine or deciduous trees was 0.7 or higher. Accordingly, a single overall model was fit to the LiDAR and field data for both acquisition years. The model was given by where i denotes an observation, i is an error term, b 0 is a common intercept term, b ik ¼ 0 for spruce dominated plots for both years, and k ¼ 1 for pine dominated plots scanned in 2007, k ¼ 2 for deciduous dominated plots scanned in 2007, k ¼ 3 for pine dominated plots scanned in 2010, and k ¼ 4 for deciduous dominated plots scanned in 2010. P 50% and P 90% are the 50th and 90th height percentiles of points above 1.37 m, and dns is the canopy density, computed as the fraction of all points in an observation that are above this height. The accuracy of the LiDAR-based AGB model was estimated by calculating the root mean square error (RMSE) of predictions by applying leave-one-out cross-validation on the plots used to create them. The RMSE of the AGB model was 32.4 t/ha, corresponding to 23.1% of the mean AGB derived from the field measurements of both years. The coefficient of determination, R 2 , was found to be 0.78. By applying the model to LiDAR metrics over the test site, and stand level dominant species information from the forest management plan, AGB maps for each year were produced with a 17.5 m Â 17.5 m pixel size to match the pixel area with that of the field plots (with a radius of 10 m). The maps were subsequently resampled through spatial averaging to 50 m Â 50 m pixels coregistered with the backscatter maps. The resampling has the effect of further reducing the prediction variance, as the resampled pixel value is essentially a mean value of about eight predictions of the model relating AGB to LiDAR metrics.
The change prediction accuracy of the AGB maps was evaluated by computing the average pixel differences between 2007 and 2010 for all pixels inside the ten 80 m Â 80 m plots, and comparing these values to the AGB differences derived from the field measurements from these plots. In the calculation, the value of a pixel partially inside a plot was weighted by the fraction of its area covered by the plot polygon. The change prediction RMSE of the LiDAR-based AGB maps was determined to be 18.1 t/ha (5.5% of the range of field measured AGB change in the ten evaluation plots) with a bias of 2.0 t/ha. Although the uncertainties of both AGB maps are combined when the difference between them is evaluated, the lower RMSE is not surprising, since the evaluation plots are much larger than the 10 m radius plots used to estimate the uncertainty in the AGB model itself.

Backscatter offset correction
The backscatter intensity varies due to environmental conditions such as moisture variations. For both Land P-band, the SAR data used in this study showed significant changes in backscatter for areas with little change in AGB between the acquisition years. Scatterplots of the differences from all image pairs of both bands are shown in Figure 2.
Most of the data points are clustered around relatively low positive values of AGB change, indicating natural growth, while larger AGB changes present in the data are from clear-cuts in stands with relatively high AGB, resulting in a tail of large negative changes in AGB. Aside from a significant spread in backscatter change for the same AGB change, we can observe that many points close to zero AGB change are not centered around zero in backscatter change, but are offset by 1-2 dB for both bands. While all image pairs have a positive offset, its magnitude differs between different image pairs, as witnessed by the horizontal striping in the plots. These differences in offset pose a problem in constructing a linear model of AGB change that would generalize across different image pairs. To mitigate these offsets, Sandberg et al. (2014) proposed a backscatter change offset correction method for P-band data. The correction was based on the HH/VV backscatter ratio, found to be correlated with biomass but relatively insensitive to changes in moisture (Soja et al. 2013;Sandberg et al. 2014). To find a suitable method for correcting L-band data, the correction method was generalized by analyzing different candidate polarization ratios.
The backscatter offset correction procedure is as follows. For each pair of scenes, find the set of 50 m Â 50 m pixels, X, that have a change in the chosen polarization ratio DR, on a decibel scale, within some threshold t. Then, choose a polarization, calculate the mean backscatter level of the pixels in X for each year, multiply each full image in the pair by the mean of the two X-means, and divide it by its own X-mean. The resulting images have the same mean backscatter in the low DR pixels for each year. For an image pair I 1 and I 2 , Equations 4 and 5 give the corrected images I 1 0 and I 2 0 , when applied separately to each of the polarizations of I 1 and I 2 : where q 1 and q 2 are the mean backscatter of I 1 and I 2 over areas with DR within t, that is, Figure 3 shows the changes in backscatter polarization ratios HH/VV, HV/VV, and HV/HH for all Land P-band pairs. At P-band, the offset for the ratio HH/VV (Figure 3d) was clearly smaller than the offset of the single polarization channels themselves (Figure 2d-f). Despite backscatter offsets due to moisture changes and possible radiometric calibration errors, the HH/VV ratio could therefore be used to detect areas with relatively low changes in AGB between SAR acquisitions, and the images from the two acquisitions could then be corrected so that the average backscatter in these areas was equal in each image. The HV/VV ratio also showed a dependence on AGB change (Figure 3e), but with a slightly negative offset, albeit significantly smaller than the offset in backscatter of the individual polarization channels. The HV/HH ratio showed no apparent dependence on AGB change and a slight negative offset.
At L-band, the HH/VV ratio was weakly correlated with AGB change, and it furthermore showed a positive change in backscatter irrespective of AGB change ( Figure 3a). As such, it was not suitable for correcting L-band data in the manner described above. Instead, the HV/VV ratio of the L-band data showed the desired characteristics of sensitivity to AGB change paired with an insensitivity for changes external to AGB change (Figure 3b). The HV/HH showed a dependence on AGB change, but also an average negative offset of about 2 dB for pixels with no change in AGB. In this study, the L-band data were therefore offset corrected using the HV/VV ratio, while the Pband data were corrected using the HH/VV ratio as in Sandberg et al. (2014).
A suitable threshold t for backscatter correction was assessed by calculating the change in mean backscatter for pixels with a small change in the AGB maps, i.e. pixels with a change in AGB on a natural logarithmic scale less than 0.1. Figure 4 shows how this remaining backscatter offset depends on the threshold used for Land P-band. The two correction methods reduced the backscatter offset for areas with small changes in AGB significantly. Many of the estimated offsets for the different channels were not within the radiometric calibration uncertainty of ±1 dB reported by DLR and ONERA. In the worst case, this could lead to a calibration offset with a magnitude of ±2 dB. For both bands, the remaining offset was greatly reduced when applying a correction using even a small threshold. Note that uncorrected data were plotted as being corrected using a threshold of zero in Figure 4. Then, after initial fluctuations as the threshold was increased, a minimum magnitude of the offset was achieved. As the threshold was further increased, the offset started to grow again as a larger portion of the images was used for the correction. Finally, as the full dynamic range of the images was within the threshold, the remaining offset became constant.
Backscatter correction using the HV/VV-ratio at a threshold of 2 dB reduced the offset in L-band backscatter from at most 2.85 dB to within 0.12 dB for all image pairs and polarizations (Figure 4). The difference in offset between image pairs for a given polarization was also greatly reduced from at most 1.77 dB to 0.12 dB or less for all polarizations. Backscatter correction using the HH/VV ratio at a threshold of 3 dB reduced the offset in P-band data from a maximum of 2.16 dB to within 0.19 dB for all image pairs and polarizations. For P-band, the largest difference in offset between image pairs was also reduced from about 1.13 dB to 0.11 dB or less for all polarizations. This reduction in offset difference is really the important metric, since the aim was to train a model on data from one pair of images, and predict AGB change on data from another image pair. The difference in offset should be small so that the constant term in the model is applicable across image pairs. Based on this evaluation, backscatter offset correction using thresholds of 2 dB and 3 dB was performed on the L-and P-band SAR data, respectively, before the model selection. It was noted that the correction was not sensitive to the precise value of the threshold, and that values within ±1 dB gave very similar correction results.
Scatterplots of offset corrected backscatter maps from all image pairs using the chosen thresholds are shown in Figure 5. Comparing to the scatterplots of uncorrected backscatter change in Figure 2, we can see that the corrections have succeeded in reducing the offset so that the scatterplots in Figure 5 are centered at the origin, and that areas of low backscatter change were roughly centered on zero change in AGB. The horizontal striping of the plots was also reduced, as the differences in offset between different image pairs were reduced.

Modeling of AGB change using backscatter
The LiDAR-based AGB maps were used in the multiple linear regression analysis to select and train backscatter change models. The evaluation plots and non-forest areas were masked out, and only pixels available in the AGB maps for both 2007 and 2010 were used. The data set then consisted of 1,355 50 m Â 50 m pixels for each image. The incidence angle for the SAR data used in the analysis ranged from 28 to 50 . The dependence of backscatter change on AGB change is not linear, and therefore a number of transformations of both AGB and backscatter were investigated in the model selection process. The transformations are an attempt to empirically find a suitable model that is linearly correlated with the LiDAR-based AGB maps. To mitigate overfitting, a cross-validation procedure was used in the model selection. The regression models investigated are generally given by where the response variable y is the AGB change on the linear, square root, or logarithmic scales defined in Equations 9, 10, and 11, respectively: DAGB ≝ AGB after À AGB before (9) The b k are the regression coefficients and is an error term. The x ij in Equation 8 are backscatter change in amplitude, power or decibel (dB) units defined according to Equations 12, 13, and 14, respectively, for a given polarization.
Different backscatter change measures were not mixed in the same model, meaning the same transformation was applied to all explanatory variables for a certain model. The described response and explanatory variable transformations result in nine different transformation combinations for the left and right hand side of Equation 8. All submodels of these with two to five non-zero regression coefficients b k were included as candidate models. Since the intercept b 0 is always included, this means all models with one to four explanatory variables were investigated. The different combinations of response variable and explanatory variable transformations, and the different possible choices of explanatory variables give a total of 504 possible models.
As a first step of two in the model selection, the models with the highest adjusted R 2 out of all models with the same number of explanatory variables and the same response variable transformation were chosen. This was done separately for each response variable transformation, since adjusted R 2 values cannot be compared across different response variable transformations. Additionally, a model was not chosen if its adjusted R 2 was not increased compared to a model with fewer variables with the same response variable transformation.
As the same model was to be used for all image pairs, we were interested in models that simultaneously fit data from all pairs. Because of this, the first step of model selection was done on the backscatter offset corrected data, with data from all image pairs concatenated into one combined data set on which the model parameters were estimated. The backscatter offset correction adjusts the relation between backscatter change and AGB change to be similar across image pairs.

Cross-validation of backscatter-based AGB change models
In the second step of model selection, an approach that was denoted leave-one-stand-out (LOSO) crossvalidation was used to validate the AGB change predictions. Using the stand delineation map, each model was trained on data from all stands but one, and then evaluated on the left-out stand. This was repeated until each stand had been used once for evaluation. The LOSO cross-validation reduced the computational cost compared to leave-one-out cross-validation at the pixel level, and aggregated alternatives such as spatial leave-one-out cross-validation, as proposed by Le Rest et al. (2014). The latter has been used to overcome spatial autocorrelation and was demonstrated in the modeling of ecological processes.
Since the models should be applicable between image pairs, the cross-validation was simultaneously conducted not only across stands, but also across image pairs. This showed the potential of using a model trained in one area on one image pair to predict AGB changes in another area from another image pair.
The procedure was repeated for all image pairs and the average performance of each model over all image pair combinations was estimated to produce an overall pixel level root-mean-square deviation (RMSD), and a stand level RMSD for each model. RMSD was used to denote comparisons with LiDAR-generated references, while RMSE was used to denote comparisons with field-based references. RMSD is calculated in the same manner as RMSE, the only difference being the used reference data set.
RMSD values were also calculated by applying LOSO cross-validation, but without switching image pairs between training and validation. The difference between these in-pair RMSDs (best-case conditions) and the across-pair RMSDs (showing temporal stability) was the accuracy loss (in terms of RMSD) when training a model on one pair of images, and predicting on another pair of images. Additionally, in-pair and across-pair biases were estimated through calculating the mean residuals for all image pair combinations. The biases were estimated through comparisons with both the LiDAR-based AGB change and the field measured AGB change in the evaluation plots.
To be able to compare models with different transformations of the response variable, the predictions on logarithmic and square root form were back-transformed and corrected for back-transformation bias prior to computing RMSDs. The predictions on logarithmic form were bias corrected using a ratio estimator according to Snowdon (1991), while the predictions on the square root scale were bias corrected following the method described in Gregoire et al. (2008).

Comparisons with LiDAR estimated AGB change
Summaries of the validation results are presented in Tables 4 and 5 (L-and P-band, respectively). These accuracies were used as criteria for the model selection. For both bands, for a given number of variables, the best performance when compared on the pixel level to the LiDAR-based AGB maps was achieved for the AGB change modeled on the square root scale and backscatter change modeled on the logarithmic (dB) scale. For these models, at L-band, the standard deviations of RMSDs from different image pair combinations were within 1 t/ha. At P-band, the corresponding values were within 0.5 t/ha. At L-band, there was a decrease in RMSD of about 4 t/ha and 6 t/ha for the pixel and stand level, respectively, when adding a second explanatory variable, but further increasing the number of variables reduced the predictive performance, indicating overfitting. At P-band, each added explanatory variable decreased the RMSDs of these backscatter models (with AGB change modeled on the square root scale and backscatter change modeled on the logarithmic scale). However, RMSDs did not decrease more than 0.1 t/ha and 0.3 t/ha for the pixel and stand level, respectively, when more than two explanatory variables were added. Based on these results, the models highlighted in Tables 4 and 5 were selected to be evaluated through comparison with the field measured AGB change in the ten 80 m Â 80 m field plots.
In Table 6 (upper half), the mean biases as compared to the LiDAR-based AGB maps, standard deviation of biases, and p values for the selected models are shown. The p value indicates whether the variance in bias across pairs (i.e. when training a model on one image pair, and using the model to make predictions on another image pair) is significantly larger than in pairs. The models are identified by a letter indicating the band and a figure indicating the number of explanatory variables. In Table 6, it can be seen that the mean biases for 8 Pixel-RMSD and Stand-RMSD are averages over all combinations of training and evaluation pairs. Columns 5-7 give the minimum, maximum, and standard deviation of the pixel level RMSDs from all combinations. The models which are highlighted were also evaluated by comparison to field measured AGB change. Definitions of the AGB and c 0 change measures are given in Equations 9-14. 3 Pixel-RMSD and Stand-RMSD are averages over all combinations of training and evaluation pairs. Columns 5-7 give the minimum, maximum, and standard deviation of the pixel level RMSDs from all combinations. The models which are highlighted were also evaluated by comparison to field measured AGB change. Definitions of the AGB and c 0 change measures are given in Equations 9-14.
L-and P-band were positive and less than 3 t/ha. Furthermore, the mean biases decreased with the number of explanatory variables. It is also worth noting that the mean biases were lower for the L-band models compared to P-band. The mean biases across pairs and in pairs differed by less than 0.1 t/ha for a given model. For P-band the variances in bias across pairs were significantly larger (p<0.05) than in pairs.

Comparisons with field estimated AGB change
Details of the models with the best predictive performance with respect to the LiDAR-based AGB maps (RMSD) for both bands are given in Table 7, where in-pair and across-pair accuracies are presented. To also evaluate the predictive performance of the models and avoid the uncertainties related to using the LiDAR-based estimates, the predictions were compared with the field estimates, and listed as RMSE in Table 7. The prediction accuracies (across-pair) from SAR (18.6-25.8 t/ha, Table 7) were similar to those from LiDAR (18.1 t/ha). The L-band model L2 (highest validation accuracy) had an RMSE of 20.8 t/ha, while the corresponding P-band model, P4, had an RMSE of 19.0 t/ha. Compared to the range of field measured AGB change in the ten evaluation plots, the RMSEs correspond to 6.3% and 5.7%, for L-and P-band, respectively. The chosen models indicate that HV was the polarization most sensitive to AGB changes, and HH-and VV-derived variables contributed only in Pband models with three or more explanatory variables, resulting in very small improvements to the RMSEs. However, the prediction accuracies across image pairs depend on the offset corrections, which cannot be performed with access only to the HV polarization.
The mean biases compared to the field measurements in the evaluation plots, standard deviation of biases and p values for the selected models are shown in Table 6 (lower half). Overall, the magnitude of the mean biases for L-and P-band were less than 2.3 t/ha. The lowest mean biases were achieved for the P-band models with more than one explanatory variable.
To compare the model prediction accuracies from different sensors and bands, the change predictions of the SAR models L2 and P4 and the LiDAR-based AGB maps for the large plots are shown in Figure 6 together with field measured AGB changes. There is no apparent overall band-related difference in the predictions, and neither a clear bias or a change in variance for the largest changes, originating from the clear-cut plots 9, 12, and 16.
Evaluation of the SAR-based predictions by comparison to the LiDAR-based AGB maps captures a more complete set of AGB change values than present in the field plots, despite the uncertainty of using LiDAR-based estimates as reference. Figure 7 shows the AGB change prediction maps of models L2 and P4 together with the LiDAR-based AGB change map. The SAR maps depict across-pair LOSO cross-validated predictions.
The training-prediction pairs used for Figure 7 are L a -L f , and P a -P f , for L-and P-band, respectively. These combinations represent worst cases in terms of backscatter offset, in that the 2007 acquisitions across the pairs have the largest time difference out of all combinations. Change maps based on other combinations are, however, similar. The predictions both with L-and P-band data overall capture the general trends in LiDAR predicted AGB change quite well, albeit with larger variance locally than the LiDAR-based change map. Table 6. Validation and prediction biases of selected models for L-and P-band (models with the lowest pixel level RMSD for a given number of variables   Figure 6. Field measured AGB change and predicted AGB change based on LiDAR, L-, and P-band SAR data. The SAR-based predictions are obtained using the models with the lowest RMSD on the pixel level, L2 and P4 in Table 7. For these models, one prediction per image pair is shown. The dominant tree species for each plot is reported above the plot ID. The predictions were temporally robust, with differences in RMSD between in-pair and across-pair predictions within 0.1 t/ha (Table 7). Correspondingly, a similar robustness was apparent for the RMSEs, which were calculated from the predictions on the field evaluation plots ( Table 7). The differences in RMSE were within 0.1 t/ha for L-band, and within 0.2 t/ha for Pband. Moreover, the mean biases (Table 6) across pairs and in pairs differed by less than 0.4 t/ha for a given model. The variances in bias across pairs were not significantly larger (p > 0.05) than in pairs.

Discussion
The relatively long wavelength of the radar systems used in this study is a strength when estimating forest biomass via SAR backscatter, due to long wavelengths penetrating deeper into the forest canopy and being sensitive to forest elements down to the order of a wavelength. Consequently, because the wavelength is on the order of decimeters for L-band, compared to meters for P-band, we may further expect L-band backscatter to be less sensitive to biomass than Pband backscatter. However, the results of this study did not show a large difference in performance between the two bands for predicting hemi-boreal forest AGB change. Instead, the predictions with L-band data achieved RMSD and mean bias values slightly lower, but close to those with P-band data, when compared to the LiDAR-based AGB maps. Yet, when compared to direct change measurements on field plots, the roles were reversed, and the predictions at L-band were slightly worse than at P-band, but still achieving RMSEs close to those of the LiDAR-based AGB maps used for estimating the parameters of the models.
The discrepancies between model predictions on the LiDAR-based AGB maps and the evaluation plots illustrate a limitation of the study in comparing the performance of the bands. As the parameters of the SAR models were estimated on LiDAR-based AGB maps, and the models achieve accuracies close to the accuracy of the training data, it may be that their performance is limited by the quality of the LiDAR-based AGB maps. The prediction accuracy achieved using SAR data is generally not better than when using LiDAR data for AGB or AGB change modeling, but as the SAR data used in the study are high resolution airborne acquisitions, the two data sources might be close to each other in information content. Higher accuracies for the SAR models would possibly have been achieved if field measurements of AGB change could have been used to train the models. However, one must often rely on LiDAR-based predictions due Figure 7. Biomass change maps, from left to right, based on LiDAR, L-, and P-band SAR data. The SAR-based predictions were obtained using the models with the lowest RMSD on the pixel level, L2 and P4 in Table 7. The black squares show the locations of the evaluation plots. The axes are labeled in meters, in UTM zone 33 N, datum WGS84 projection.
to the lack of field plot inventory for largearea mapping.
The change prediction results presented in Figure 6 show no decrease in accuracy even for the largest changes of around 300 t/ha. This contradicts the saturation effect noted in earlier studies of biomass state prediction (Dobson et al. 1992;Rignot et al. 1994;Imhoff 1995;Fransson 1999), which indicates that it could be difficult to measure changes below about 100 or 200 t/ha. This result should be further evaluated on other data sets. However, some previous studies have indicated that the saturation for high biomass was not caused by the high biomass itself, but from different structural changes in the forest as the biomass increases, with some of them increasing the backscatter, and others decreasing it Joshi et al. 2017). Another important reason for saturation is due to ground topography (Soja et al. 2013), which is relatively modest at the Remningstorp test site and thus increases the point of saturation. As almost all pixels with a biomass change below 100 t/ ha in the test site are from a handful of stands, the fact that the models can capture unexpectedly large changes may be due to the specific structure of the stands or area in question.
The offset correction and change prediction results were obtained in a hemi-boreal forest that is relatively flat in topography, and while the results are promising, they may not generalize to forests growing on more undulating terrain in the boreal region. This type of forest is relatively sparse, and a non-negligible ground contribution is expected to disturb the signal more with increasing topography. Insights could also be gained from comparing the performance of the offset correction methods for different tree species. While the correction methods may not be suitable for other forest types, studies evaluating the offset correction methods in temperate and tropical forests could give further insights about the underlying processes. In high biomass tropical forests, it is possible that the HV/VV-based correction could be better suited for also correcting P-band data, as the SAR signal does not penetrate as far into the dense canopy, making the interaction more like that of L-band scattering in less dense forests. For this test site, the P-band HV/ VV ratio showed a slight offset in backscatter, but has some dependence on biomass change (Figure 3), and could therefore give a reasonable offset correction for P-band.
The most important contribution of this study was the development of a backscatter offset method for Lband data, which facilitated the prediction of biomass change between images acquired in different moisture conditions or with radiometric calibration differences between them. Specifically, the proposed correction enabled predictions across image pairs with similar accuracies as those obtained within an image pair. An important finding was that a correction based on the HH/VV polarization ratio, as previously developed for P-band data, was not suitable for L-band in the same forest and environmental conditions. Instead, it was found that for the L-band data, the HV/VV ratio had the qualities needed for detecting low change areas, and the offset correction method developed was based on this ratio. The offset corrections made predictions across image pairs almost as good as those within each pair. Yet, the remaining variance in predictions on the evaluation plots from different image pairs indicate a more complex relationship between moisture and backscatter intensity than that of a simple offset. Understanding these relationships considering the underlying backscatter mechanisms might enable even better correction methods and more precise SAR-based AGB change predictions. The correction methods developed and applied in this study remove an offset due to moisture or calibration error that is constant over the whole scene, but cannot remove local variations.

Conclusions
We have evaluated the performance of a method for predicting AGB change from changes in P-and Lband SAR backscatter data in hemi-boreal forests. An important element of the method is a backscatter-offset correction, which is based on a polarization ratio and does not require any reference ground data. By training regression models for each band on AGB maps created from field surveys and LiDAR data, we showed that both L-and P-band backscatter could predict AGB changes with errors close to those of the LiDAR-based AGB maps used for training the models. The prediction accuracies of L-and P-band models differed only slightly, although this result may have been due to limitations of the LiDAR-based AGB maps used to train the models.
The P-band results were in line with earlier studies, verifying the results of Sandberg et al. (2014). The method and evaluation for L-band presented in this paper is new and the results show almost the same prediction errors as for P-band. Furthermore, it was empirically found that the HH/VV-based backscatteroffset correction used for P-band was not suitable for correcting L-band data, and instead an HV/VV-based correction performed better. The reason for this is yet to be fully understood, but as P-band SAR is known to interact more strongly with the ground and stems of the forest, while L-band SAR backscatter is mainly from the upper parts of the canopy (Tebaldini and Rocca 2012), it is possible that the HH/VV-based correction mainly corrects for ground moisture and that L-band data instead need to be corrected for canopy moisture.
The proposed backscatter offset correction method for L-band data was shown to facilitate prediction results across different image pairs almost identical to those obtained when training and predicting on the same image pair. Finally, while the HV/VV-based offset correction worked well for the L-band data in this data set, more studies are needed to see if this generalizes to other data sets and regions. Provided that it does, and that future L-band SAR missions provide HV and VV data, the method can alleviate and improve the mapping of AGB change using L-band SAR data. The findings suggest a potential for large-area mapping of forest biomass change using, e.g. the upcoming satellite SAR missions ALOS-4 and NISAR (L-band), and BIOMASS (P-band), despite varying environmental conditions and calibration uncertainties.