Assessment of Forest above Ground Biomass Estimation Using Multi-Temporal C-band Sentinel-1 and Polarimetric L-band PALSAR-2 Data

Abstract

Synthetic Aperture Radar (SAR), as an active sensor transmitting long wavelengths, has the advantages of working day and night and without rain or cloud disturbance. It is further able to sense the geometric structure of forests more than passive optical sensors, making it a valuable tool for mapping forest Above Ground Biomass (AGB). This paper studies the ability of the single- and multi-temporal C-band Sentinel-1 and polarimetric L-band PALSAR-2 data to estimate live AGB based on ground truth data collected in New England, USA in 2017. Comparisons of results using the Simple Water Cloud Model (SWCM) on both VH and VV polarizations show that C-band reaches saturation much faster than the L-band due to its limited forest canopy penetration. The exhaustive search multiple linear regression model over the many polarimetric parameters from PALSAR-2 data shows that the combination of polarimetric parameters could slightly improve the AGB estimation, with an adjusted R² as high as 0.43 and RMSE of around 70 Mg/ha when decomposed Pv component and Alpha angle are used. Additionally, the single- and multi-temporal C-band Sentinel-1 data are compared, which demonstrates that the multi-temporal Sentinel-1 significantly improves the AGB estimation, but still has a much lower adjusted R² due to the limitations of the short wavelength. Finally, a site-level comparison between paired control and treatment sites shows that the L-band aligns better with the ground truth than the C-band, showing the high potential of the models to be applied to relative biomass change detection.

1. Introduction

Mapping forest Above Ground Biomass (AGB) using Earth observations has become an integral part of monitoring and assessment programs. For many years optical data have become embedded as main sources of information driving forest decision support tools—for example, time series Landsat and MODIS phenological indicators within tools like Global Forest Watch and Forewarn (https://www.globalforestwatch.org/ and https://forwarn.forestthreats.org/). These metrics have tended to focus on forest vs. not forest spectral differences delineated in time series and degradation or anomaly alerts based on deviation compared to historical patterns. While optical-based mapping techniques have made significant contributions toward operational forest assessment, less effort has been placed on using radar observations within operational decision support tools due to less systematic coverage available and more complex data structures [1].

Synthetic Aperture Radar (SAR), as an active observation technique, can transmit longer electromagnetic wavelengths from 1 mm to 1 m and receive the scattered waves after interacting with targets. SAR data are often noted as having value because of their all-weather, day-and-night capability, the ability to penetrate clouds, and sensitivity to structural attributes (surface roughness, orientation, and moisture). Due to backscatter saturation issues of shorter wavelength X- and C-bands [2,3], relatively longer wavelength SAR systems such as L- and P-band have been promoted for characterizing AGB in more dense or intact forested landscapes for many years [4,5,6]. P-band shows a strong ability to accurately map AGB; however, no spaceborne P-band instrument is readily available. Since the launch of the JERS-1, the science community has used L-band to map forested landscapes such as the Global Rain Forest Mapping Project [7]. The launch of ALOS-1 PALSAR-1 by JAXA in 2007 enabled additional observations to be used for mapping forest AGB [8,9,10,11] and this thematic focus continues with ALOS-2 PALSAR-2 launched in 2014. With the planned launch of the open access NASA ISRO SAR (NISAR) Mission in a few years, the availability of L-band data will continue to grow. At L-band (1.25 GHz; 23 cm wavelength), radar waves penetrate into a forest canopy through gaps and areas with relatively low density of scatterers (sparse medium effect) and scatter back from various components of the forest such as the crown (leaves and branches) and stems [12]. The co- and cross-polarized radar measurements contain different information about the orientation and structure of forest canopy and tree stems within the resolution cell of radar images [13] making these observations particularly useful. Hence, this wavelength domain has particular strengths for forest structure mapping.

The use of dual pol L-band, such as stamps or open access global JAXA mosaics, have provided opportunities for mapping forest AGB and several studies have noted saturation thresholds. For example, Bouvet et al. [14] used 2010 PALSAR-2 mosaics to estimate woody AGB in African savannahs and woodlands using Bayesian methods. This work notes limitations of estimation AGB over 85 Mg/ha. Similar AGB thresholds (<100 Mg/ha) were described by Ningthoujam et al. [15] who used HH and HV polarizations of PALSAR-1 for woody biomass estimation in tropical deciduous mixed forest. Mitchard et al. [16] used co-polarized L-band data to map AGB across several sites of woody vegetation up to 150 Mg ha⁻¹ before the relationships became noisy from saturation. Hame et al. [17] found L-band effective in mapping AGB in Laos with dense tropical forest but noted many of those sites had undergone degradation. Outcomes from these and other studies have highlighted how dual polarizations have limitations and saturation thresholds for mapping forested AGB.

The full polarimetric SAR signature provides additional parameters other than the single polarization channels. Many decomposition methods have been developed to take advantage of this additional information [18,19,20,21]. The phase difference and polarimetric coherence have proved to be very useful for AGB estimation [3,22]. For example, Antropov et al. [23] recently found that the combination of the co-polarization coherence and HV had higher accuracy in AGB estimation than the single or dual polarizations. However, a traditional limitation is the lack of operational full polarimetric satellite observations.

The recent launches of Sentinel-1A/B now provide open-access and operational C-band observations to the community. This has created new opportunities to use different satellite microwave observations now that availability is growing. Historically, many studies have found that blending different radar wavelengths improves forest AGB mapping. For example, Saatchi et al. [24] used an empirical approach to map biomass and fuel loads in Yellowstone using L- and C-band radar. In this research application we build on a lineage of radar AGB mapping and relative biomass change detection by evaluating the use of multitemporal Sentinel-1 and quad pol PALSAR-2 for mapping AGB in the northeast USA. Objectives include testing backscattering amplitude, polarimetric parameters derived by recently developed decomposition methods, use of multi-temporal Sentinel-1, and relative biomass change detection over control and treatment sites. A variable selection framework and a semi-physical scattering model approach were applied to comprehensively evaluate these data for AGB mapping in a complex landscape.

2. Study Area and Data Collection

2.1. AGB Data

A field campaign was carried out in northern New England, USA (NNE) to obtain a robust set of calibration and validation measurements (

Table 1. Forest, sites, and plots.

Forest Property Full Name	Forest Short Name	Site	Sample Plots	Description
Bartlett Experimental Forest	BAR	07_C	9	White pine-northern hardwood reference
		07_T	9	White pine-northern hardwood thinning
		22_C	9	Northern hardwood reference
		22_T	9	Northern hardwood thinning
		35_C	9	Northern hardwood reference
		35_T	8	Northern hardwood thinning
Bear Camp River Property	BEAR	T	9	Oak-pine thinning
Harmon Preserve	HAR	C	9	Pitch pine barrens reference
		T	9	Pitch pine barrens thinned and burned
Kingman Farm	KING	C1	9	Oak-pine reference
		T1	9	Oak-pine thinned
Lord Farm	LORD	LORD	9	Oak-pine reference
Lovell River	LOVE	C1	9	Oak-pine reference
		CC2	7	Oak-pine clearcut
		T1	6	Oak-pine thinned
Massabesic Experimental Forest	MAS	C2	9	White pine-hemlock reference
		HD	9	White pine high-density thinning
		LD	9	White pine low-density thinning
		RP	9	Red pine plantation
Mendums’ Pond	MEND	C	9	Oak-pine reference
		T	9	Oak-pine thinning
Jones Forest	MIL	C	9	Oak-pine reference
		T	9	Oak-pine thinning
Moore Fields	MOOR	C1	9	Oak-pine reference
		T1	9	Oak-pine thinning
Pine River State Forest	PR	C2	9	White pine reference
		T1	9	White pine tornado

). Most of the natural forested landscape in the region is broadly classified as Eastern Temperate Forest with Level 3 ecoregions of Atlantic Maritime Highlands, Northeastern Coastal Zone, and Acadian Plains and Hills according to the U.S. EPA Ecoregions [25]. Within this region are flat coastal plains and extreme topography in the mountainous areas, providing a range of conditions. The tree species composition is highly variable, and mixed-species stands are the rule rather than the exception; common forest types include northern hardwoods (dominated by Acer saccharum, A. rubrum, Fagus grandifolia, and Betula allegheniensis), oak‒pine‒hemlock (dominated by Quercus rubra, Pinus strobus, and Tsuga canadensis), and spruce‒fir (dominated by Picea rubens and Abies balsamea), with occasional stands of other types such as red pine (Pinus resinosa). The climate in the region is classified as continental wet warm (Dfb: Koeppen classification) with warm, humid summer conditions and cold winters. Strategic field plots were coordinated with regional environmental and government agencies, private industry lands, and university partnerships. A total of 27 sites were selected within 10 different properties (Table 1). In most cases, these sites were selected based on the occurrence of a known disturbance (natural or anthropogenic), with sufficient spatial extent to allow installation of a grid of measurement plots with a buffer of at least 60 m to the edge of the disturbance, and an adjacent site reflecting forest structure similar and composition similar to that which had been present prior to the disturbance. For example, at the Bartlett Experimental Forest, site 07_T is a recently thinned northern hardwood forest, while 07_C is an adjacent “control” or reference condition. The Bear Camp River property was too small to install both disturbed and reference sites, so a reference site was installed at Lord Farm only a few km away. The Massabesic Experimental Forest included thinnings to low and high density (LD and HD respectively), as well as a red pine plantation (RP) and a reference site (C). Each site includes a systematic grid of nine plots wherever space and the required buffer allowed, with plot centers spaced 30 m apart; at each plot center, a variable-radius plot (also known as a horizontal point sample) was installed, using a basal area factor of 4 m²/ha [26]. For all live and standing dead trees, the species and diameter at breast height (DBH, cm) were recorded; a variety of other measurements were also taken but are not reported here. Aboveground biomass (AGB) for each tree was calculated using the allometric equations of Chojnacky et al. [27].The variable radius plots represent a statistical subsample of the trees within each site, which can be scaled to represent population-level estimates using standard techniques [26]; biomass estimates from variable-radius plots have shown correlations comparable to, or in some cases better than, fixed-area plots, with airborne lidar-derived metrics in this study region [28]. At each plot within the paired disturbed/treatment sites (T) and matched reference control sites (C), georeferenced measurements were taken following best practices [26]. Figure 1 illustrates an example SAR pixels matched with plots in the Massabesic Experimental Forest (MAS), which has four sites of varying biomass: a control site (C2), high-density cut (HD), low-density cut (LD), and Red Pine planation (RP).

2.2. Sentinel-1

Dual polarization (VV+VH) Sentinel-1 Interferometric Wide (IW) Ground Range Detection (GRD) data were obtained from the Alaska Satellite Facility (ASF) for this application with its coverage (Figure 1). Sentinel-1A carries a C-band imager at 5.405 GHz with an incidence angle between 20° and 45°. The platform follows a Sun-synchronous, near-polar, circular orbit at a height of 693 km. These C-band radar data were first transformed into Beta Naught to assist in radiometric terrain correction and thermal noise removal using the look-up table (LUT) within the metadata. Radiometric distortion caused by the terrain slope in this region was corrected using the method proposed by Small [29] to achieve Gamma Naught γ°, which is less dependent on the incidence angle than the Sigma Naught σ°. The inherent speckle noise was filtered using a Boxcar filter method with a 5-by-5 window size. Finally, with the assistance of the precise orbit and 30 m Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM), all data were geocoded using the range-doppler approach with the output resolution of 15 by 15 m (to match a nominal application scale in this case). All radar preprocessing was executed using open access routines (https://github.com/Applied-GeoSolutions/gips and https://github.com/Applied-GeoSolutions/gippy) based in Python, C, and GDAL and are available for sharing (please contact the authors).

2.3. PALSAR-2

Fine beam quad polarimetric single look complex PALSAR-2 data were obtained from JAXA over the field sites (Figure 1), and it is also processed in the same way as Sentinel-1 data except that the output of the PALSAR-2 is the 3 × 3 coherency matrix. Thus, many polarimetric parameters could be derived from the coherency matrix using the polarimetric decomposition and other polarimetric related methods. We evaluated parameters from two decomposition methods including both eigen and model-based approaches developed by Touzi [21] and Huang et al. [30], respectively. In addition, select geometric parameters such as the geometric randomness and shape factor to describe the structure of the forest canopy developed by Huang et al. [31] are also used in this paper. All available parameters (

Table 2. SAR parameters from Sentinel-1 and PALSAR-2, in which the yellow highlighted parameters are the backscatter amplitudes that have been extensively investigated previously.

Sensors	Parameters		Description
Sentinel-1	VH		VH channel
	VV		VV channel
PALSAR-2	VH		VH channel
	VV		VV channel
	HH		HH channel
	Touzi Decomposition [19]	Lambda 1	First dominant scattering amplitude
		Lambda 2	Second dominant scattering amplitude
		Lambda 3	Third dominant scattering amplitude
		Alpha 1	Dominant polarization angle
		Alpha m	Mean polarization angle
		Phi 1	Phase difference
	Eigen Decomposition [16]	Entropy	Scattering Entropy
		Alpha	Scattering Mechanism
		A	Anisotropy
	SVSM [28]	RVI	Geometric randomness
		Sigma	Shape factor
	SAVSM [27]	Ps	Surface Scattering
		Pd	Double-bounce scattering
		Pv	Volume scattering
		Beta_Ps	Beta parameter related to surface scattering
		Alpha_Pd	Alpha parameter related to double-bounce

) that will be investigated in this paper.

3. Methods

In this section, two modeling frameworks (Figure 2) are explored: a semi-physical simple water cloud model is used first for the investigation of the relationship between AGB and the backscattering amplitudes (VH and VV) for both Sentinel-1 and PALSAR-2 data (Framework 1), and a multiple linear regression using an exhaustive search algorithm is chosen to obtain the best model for AGB estimation only for PALSAR-2 based on a number of polarimetric parameters from both old and recent developed decompositions.

3.1. Simple Water Cloud Model

The Water Cloud Model (WCM), proposed by Attema and Ulaby [32], is a semi-empirical model assuming that the vegetation consists of a collection of spherical water droplets that are held in place structurally by dry matter. The primary assumption of the WCM is based on the fact that the dielectric constant of dry vegetation matter is much smaller than that of the water content of vegetation, and more than 99% of the volume is composed of air in the vegetation canopy [32]. Therefore, such a model was developed assuming that the canopy “cloud” called the water cloud contains identical water droplets randomly distributed within the canopy. In this paper, the Simple Water Cloud Model (SWCM) of Cartus et al. [9] is adopted for the initial exploration, which is written as

$${\gamma }_{veg}^0(B_{df})=\frac{{\gamma }_{df}^0-{\gamma }_{gr}^0{e}^{-\delta B_{df}}}{1-e^{-\delta B_{df}}}$$

(1)

where $B_{df}$ is the biomass, ${\gamma }_{df}^0$ and ${\gamma }_{gr}^0$ are the pure scattering from the forest canopy and ground respectively; $\delta$ represents the extinction coefficient, and ${\gamma }_{veg}^0$ is the measured backscattering amplitudes of the co- or cross-polarization. Note that ${\gamma }_{df}^0$, ${\gamma }_{gr}^0$ and $\delta$ are unknown parameters that need to be determined by the measured AGB and ${\gamma }_{veg}^0$. To do this, a non-linear regression is perform using the measured ${\gamma }_{veg}^0$ plot-level AGB and under the assistance of the Python scipy curve_fit function. Finally, the biomass will be written as

$$AGB=-\frac{1}{\delta }ln\left(\frac{{\gamma }_{df}^0-{\gamma }_{veg}^0}{{\gamma }_{df}^0-{\gamma }_{gr}^0}\right)$$

(2)

3.2. Exhaustive Search Multiple Linear Regression

We framed the comparison of the many polarimetric parameters provided by the PALSAR-2 in this paper as a variable selection problem. There are many methods to carry out variable selection. We select an exhaustive search variable selection of a multiple linear regression model because of the ease of interpretability and ability to select a parsimonious model. The exhaustive search algorithm implementation is from the package leaps [33] in R [34], which enumerates all possible candidate models with a given number of parameters and ranks them using several information criteria. In this study, polarimetric parameters are compared by the adjusted R² of the model in which they occur and by the frequency that the parameter occurs in all best model of a given size (1 to 19). Our emphasis here is on the identification of useful predictors among the set of candidate variables, rather than formal development of a final predictive model.

4. Results

4.1. Simple Water Cloud Model

A saturation level comparison between Sentinel-1 and PALSAR-2 was executed to get a sense of capabilities and scattering properties. For this comparison only the VH and VV polarizations are investigated since the Sentinel-1 only has the VH and VV polarizations to make Sentinel-1 and PALSAR-2 comparable. The semi-physical SWCM is used to investigate the relationship between polarizations from both C- and L-band and AGB. Near simultaneous PALSAR-2 data acquired on 15 May 2017 (day of year 135) and C-band Sentinel-1 data acquired on 14 May 2017 (day of year 134), which has only one day difference, were used in this subtask to avoid weather effects such as rain and snow. The fitting lines of VV (red line) (Figure 3) show a much flatter range than those of the VH for both Sentinel-1 and PALSAR-2, indicating that VH is much more sensitive to AGB than VV in this landscape. In addition, comparisons of the VH fitting line between Sentinel-1 and PALSAR-2 show that the C-band Sentinel-1 data reaches saturation at around 50 Mg/ha, which is orders of magnitude faster than PALSAR-2, with a saturation point around 150 Mg/ha. As expected, this is due to the shorter C-band and much higher attenuation coefficient than that of PALSAR-2, and it is almost 6 times more than that of the PALSAR-2. Therefore, the L-band backscattering amplitude data are more useful than C-band over higher AGB estimation in this landscape.

4.2. AGB Estimation from Polarimetric PALSAR-2 Data

In the above section, the VH polarization shows much higher potential than the VV polarization from the L-band to estimate the AGB. In addition to the VH and VV polarizations, other parameters (Table 2) from the polarimetric PALSAR-2 data are also investigated to further study the capability of the quad polarization data on the AGB estimation. Figure 4 illustrates (left) all one-parameter models sorted by strength (adjusted R²) and the right one is adjusted R² of best models of size n (1 to 19) overlaid with indicated matrix of parameters included in each best model.

Among the single parameter models, the Pv derived from the model-based decomposition developed by Huang et al. [30] recently from the L-band PALSAR 2 has the highest adjusted R² among all one-parameter models with its value greater than 0.3 and RMSE of 76.72 Mg/ha. The VH has an adjusted R² that is slightly lower than Pv. This is because the measured VH polarization represents the multiple scattering other than just the volume scattering from the forest canopy, instead, the Pv derived from the decomposition represents the volume scattering after removing the contribution from the surface and double-bounce scattering. The second-best parameter is the Lambda3, which represents the third dominant scattering amplitude reflecting either the double-bounce or the surface scattering amplitude over the forest areas. In terms of the best mode with two parameters, R² of the combination of Alpha from the eigen based decomposition and the Pv reaches the highest with its value of approximately 0.4 and RMSE of 70.46 Mg/ha (

Table 3. RMSE of best models of size 1 to 5.

Model ID	RMSE (Mg/ha)	Parameters
1	76.72	Pv
2	70.46	Alpha + Pv
3	69.00	Alpha + Pv + VH
4	69.68	Alpha + Lambda1 + RVI + VH
5	68.90	Alpha + Anisotropy + Lambda2 + Lambda3 + Pv

). This indicates that the Alpha value, which represents the angle between the surface and double bounce scattering, contributes a lot to the AGB estimation in addition to the Pv. This is likely due to how the AGB adopted in this research application depends on the trunk size and tree species, while the Alpha angle also represents the scattering from the trunk to some degree. As expected, when more parameters are added to the regression, the adjusted R² first increases, later remains stable, and finally slightly decreases, which is due to the overfitting caused by the intercorrelation among those polarimetric parameters.

The selection of a final best model took into consideration model interpretability, parsimony and issues of multicollinearity. Given the many parameters are derived from varying decompositions of the same data, several parameters are expected to be highly correlated (Figure 5). High correlation among model parameters results in variance inflation meaning that significant parameters appear insignificant and well as bias and overfitting. To avoid the overfitting issue, we note that there is little increase in the adjusted R² value when more than two parameters are used. Furthermore, the parameters Alpha and Pv that make up the best two-parameter model are the two parameters occurring most frequently in all best models of a given size. Therefore, for the AGB estimation, Model 2 is used here with the equation

$$AGB=e^{0.066Alpha+0.289Pv+3.106}$$

(3)

Both Alpha and Pv are significant, with p-values less than 1 × 10⁻⁶.

4.3. AGB Estimation from the Single and Multi-Temporal Sentinel-1 Data

As noted in Section 4.1, the C-band Sentinel-1 data have a severe saturation problem on day 134. In this section, we attempt to test the ability of the time dimension on the AGB estimation. Figure 6 indicates that R² using the data on different days shows a much high fluctuation ranging from 0.0 to 0.25. This might be caused by the weather conditions or the satellite orbits geometries [16,35,36]. In addition, comparison of the estimated AGB using the single day and seasonal mean over the MAS forest (Figure 7), shows that many coherent noise exists in a single day AGB map, but the seasonal mean is much smoother. Therefore, instead of using the single day data, the temporal mean gamma naught shown in red line in Figure 6 is used for the AGB estimation. Compared with the estimation using the data on day 134, the R² of the temporal mean is significantly improved, with its value as high as around 0.15. However, even though the estimation is significantly improved, it is still too low to be used for the AGB estimation. From this perspective, we might conclude that the C-band Sentinel-1 backscattering amplitude is difficult to be used for an operational use of the AGB estimation.

4.4. Relative Biomass Difference Detection at Site Level

In this section, we compare the models using a site-level comparison where we focus on ability to detect a relative difference in biomass instead of ability to estimate exact biomass. We test whether paired sites (control, treatment) that show a significant difference in field measured AGB also show a significant difference in model predicted AGB, and that the difference is in the same direction. The results of the comparison of sites where both S1 and P2 extent are available are shown in Figure 8. Of the paired sites, HAR, KING, MIL, and MOOR, only HAR shows a significant difference (5% level) in field-measured AGB, while both P2_ALPHA_Pv and seasonal mean S1 show significant differences between treatment and control sites for all four pairs. Despite the significant level of difference, all directions of differences are consistent with field measurements where treatment sites have less biomass than control sites. Formal comparison tests were not carried out for MAS and LOVE sites. However, boxplots and AGB maps (Figure 9) show that the P2 model is more consistent than the S1 model at following the trends between sites at MAS while the S1 model appears to be a better fit to the relative differences between LOVE sites. Some of this difference may be contributed to forest type. MAS is predominately conifer while LOVE is predominately deciduous (see Figure 10.)

5. Discussion

Taking the fact that the water occupies substantial volume in the vegetation canopy, the SWCM is applied as a “water cloud”, which is straightforward and widely accepted. In addition to the forest AGB estimation, it also has been extensively used for both the vegetation indices (leaf area index, biomass, etc.) and ground soil moisture estimation over agricultural areas [37,38,39]. However, this model might produce uncertain results when inverted especially when the observations fall outside the range covered by the model curve, which will result in either infinite or negative estimates [11,40]. Taking the P2 VH polarization (Figure 3c) in our case as an example, the maximum and minimum SAR backscatter of SWCM are −11.47 dB and −16.22 dB, which corresponds to the scattering from the pure forest and bare soil, respectively. According to Equation (2), any the backscatter greater than −11.47 dB or less than −16.22 will not be inverted because it results in a negative or infinite value, which is not consistent with reality. In this research application, the SWCM is used simply to investigate the saturation issues of both C- and L-band, in which the C-band saturated much faster than the L-band.

In terms of the polarimetric parameters, the measured HH, VH, and VV polarizations have been extensively investigated previously on the AGB estimation using C-, L- and P-band data. VH is shown to be more useful for the AGB estimation than VV and HH due to less influence by soil moisture, which is also consistent with our findings in Section 4.2. However, the azimuthal terrain slope will increase the measured VH [41], which makes its AGB estimation very problematic over mountain areas. Instead, some parameters derived from the recent developed methods from polarimetric SAR would be insensitive to this effect and present “purer” vegetation scattering components from the vegetation canopy than VH itself. One example is the decomposed volume scattering component (Pv), which is the best single parameter model rather than VH. To fully understand the capability of the polarimetric SAR, a total of 19 parameters, including old and newly developed polarimetric parameters, are tested in our paper. In Section 4.2, an improvement of the AGB estimation is indeed observed using the polarimetric parameters with the R² of around 0.4 and the RMSE of around 70 Mg/ha, which is consistent with other studies [4,11,42]. However, the R² is still very low and RMSE is also high, which complicates the operational use. From this perspective, we might infer that the polarimetric parameters could only slightly improve the AGB estimation when all biomass (low and high) are considered. This conclusion is still meaningful because to map the low AGB, the polarimetric parameters might be better than the VH alone, while to map the high AGB, the polarimetric parameters are not sufficient and we must seek other methods such as InSAR, PolINSAR, or TomoSAR [43].

The primary reason of the low accuracy is that our models consider both low and high biomass, while the L-band backscatter might only be suitable to map young, sparse forest with low biomass [11]. Note that when only the low biomass less than 200 Mg/ha is considered, the RMSE reduces significantly with its value of 31.18 Mg/ha. However, most of those plots are from one or two sites, which are highly spatial correlated and might result in a nonreliable result. Therefore, the low biomass estimation is not investigated in this paper. In addition, the tree species might be another reason to low the accuracy as [42] obtained an R² of around 0.5 over the pure Slash Pine, which is only slightly better than our results. Kellndorfer et al. (2003) [42] investigated the capacity of the C-band, L-band, and P-band on the retrieval of many biometric parameters in the plot level over Slash Pine, in which the L-band cross-polarization achieves the adjusted squared correlation of 0.52, which is close to our results (0.43). The difference lies in the fact that they focus on Slash Pine trees, which are pure forest, while most of the sites in our study area consist of mixed forest. Our results are consistent with the results [44] for the mixed species and complex landscapes in southern Sweden. Peregon and Yamagata (2013) [11] investigate the ALOS/PALSAR-1 L-band data to estimate the AGB in Western Siberia, stating that the saturation level may depend on the tree species and forest type as well as the ground surface type and weather conditions. When all the SAR backscatter are used, the R² is between 0.35 and 049, which is also consistent with our results. In this paper, the effects of tree species are not considered since most of the sites in our study area are a mixture of deciduous and coniferous trees (Figure 10).

To our knowledge, multi-temporal SAR data have rarely been used for forest AGB estimation. While some studies show low correlation (less than 0.1) between the AGB and backscatter of C-band [42], our results show on some days the correlation achieves as high as 0.25. Although it is still low due to the penetration issue, it still tells us on what days we could obtain a high-accuracy AGB map when longer wavelength SAR data (such as L- and P-band) are available. In addition, we found that the aggregation mean of the multi-temporal data could obtain a reliable result and reduce the coherent speckle noise. Although we did not obtain a high-accuracy estimation of the AGB map using both C- and L-band, the statistic T-testing of the AGB over the control and treatment sites from best models of L-band polarimetric parameters and of the aggregated multi-temporal C-band data, show that those models have high potential to detect the relative biomass change.

To improve the AGB estimation for operational use, in terms of an SAR system, the much longer wavelength P-band SAR systems, which have a higher penetration depth than the L-band, would be a better option. In addition, the fusion of the SAR parameters and the optical sensors might be another choice for biomass estimation. Laurin et al. [44] show that the dual polarization Sentinel-1 data, together with optical multispectral data, could improve the AGB estimation accuracy, and the regression achieves an R² of 0.7 in central Italy when an upper bound of 400 Mg/ha is considered. Most importantly, the multi-temporal method not only reduces the coherent speckle noise, but also gives a reliable estimation, as discussed in this paper. Unfortunately, we did not obtain high accuracy due to the limited penetration of the short wavelength (5.3 cm) C-band in this paper, even though multi-temporal data were used. Using long wavelength (23 cm) L-band PALSAR-1 data, the improvement of AGB estimation in central Finland was observed in [23] when the multitemporal aggregation of selected PolSAR scenes was performed.

6. Conclusions

This research application build on a lineage of SAR AGB estimated by evaluating multitemporal C-band dual-polarization Sentinel-1 and L-band polarimetric PALSAR-2 data in a northeastern forested landscape. Results derived from polarimetric PALSAR-2 data are consistent with other research results, but the accuracy is still too low to use operationally. This also demonstrates that the polarimetric parameters could only slightly improve the AGB estimation. Using the multi-temporal Sentinel-1 data, the highest R² is around 0.25 on day 250, which is still not very low, making it unsuited to operational use. A qualitative comparison between the control and treatment sites in all forests are analyzed, showing that the estimated AGB by L-band is more consistent with the ground truth than that of the C-band, except for sites at LOVE. In future, to improve the AGB estimation, the fusion of the long wavelength SAR (L- or P-band) and optical data, which provide more information in terms of structure and spectrum, will be further investigated.