BIOMASS ESTIMATION ALONG A CLIMATIC GRADIENT USING MULTI-FREQUENCY POLARIMETRIC RADAR VEGETATION INDEX

: Monitoring biomass changes in the Mediterranean ecosystems is important for better understanding their responses to climatic and anthropogenic changes. Multi-frequency SAR data (L-band PASLAR and C-band Sentinel-1) is investigated for developing a biomass estimation model along a climatic gradient in Israel. First, the relationships between biomass and transmissivity at various frequencies are discussed. Multi-frequency polarimetric radar vegetation index (MPRVI) is then proposed utilizing the ensemble average of degree of polarization and cross-polarized backscattering coefficients. After randomly partitioning the data for training set and testing set, the new MPRVI-based biomass model is evaluated. It shows a good agreement with reference biomass data with an r-square of 0.899 and a root-mean-squared error (RMSE) of 0.381 kg/m 2 and with a relative RMSE (RRMSE) of 10.8%.


INTRODUCTION
Mediterranean ecosystems along climatic gradients are characterized by high diversity of plant such as woodlands, dense shrublands, open shrublands, dwarf shrublands, and arid vegetation (Safriel, 2009;Shoshany, 2000). Mapping biomass for extensive regions along the Mediterranean-to-arid climatic gradient is instrumental for studying desertification and biodiversity changes and may play an essential role in understanding water availability, ecosystem changes, and their response to the global carbon cycle and climate change (Chang and Shoshany, 2017;Santi et al., 2017;Shoshany and Karnibad, 2015). Remote sensing is a primary information source for assessing such responses at wide regional extents (Shoshany, 2000). However, the extraction of biophysical information, including biomass, from remotely sensed data, is complex in Mediterranean areas due to their high geodiversity (high fragmentation and heterogeneity) (Santi et al., 2017;Sternberg and Shoshany, 2001). Optical sensor data have been widely used at different spatial resolutions; however, it provides information mainly for the upper vegetation layer because of its shallow penetration into the moderate-to-highly dense vegetation canopies (Guyot et al., 1989;Lu, 2006;Prabhakara et al., 2015;Shoshany and Karnibad, 2015;Chang et al., 2022). Since microwaves signal penetrate deeper into canopy layer than optical radiation, synthetic aperture radars (SAR) have been used for remote sensing of biomass, and it was proven that radar parameters such as backscattering coefficient and polarimetric coherence offer relatively high correlation with biomass (Le Toan et al., 1992;Neumann et al., 2012;Sarker et al., 2013). While radar backscatter models are well developed with input parameters including forest and shrubland biomass, the assessment of above-ground biomass (AGB) over a wide region in the Mediterranean is most challenging due to the high spatiotemporal variability of shrubs and background effects such as topography and surface roughness (Eisfelder et al., 2012;Lu, 2006a;Shoshany et al., 2000).
Polarimetric radar vegetation index (PRVI) was proposed by Chang et al. (2018) utilizing the degree of polarization (DOP) and the cross-polarized backscattering coefficients: (1) where DOP is the average of DOPH (degree of polarization for hh-polarization) and DOPV (degree of polarization vvpolarization) for full (or quad) polarization SAR data, and is the cross-polarized backscattering coefficient. PRVI shows a good performance for estimating biomass, by reducing the background effect and moderating the direct scattering effect (Chang et al., 2018). In addition, many studies have indicated that combining multi-frequencies (different wavelengths) shows a better performance than a single frequency for mapping wide range of biomass levels (Ferrazzoli et al., 1997;Jiao et al., 2010;Naidoo et al., 2015;Saatchi et al., 2007;Santi et al., 2017). Accordingly, this study aims to develop a biomass estimation model based on PRVI by using multi-frequency SAR data for wide regional mapping along a climatic gradient in the South-Eastern Mediterranean region. Taking into account the varying transmissivity characteristics of different frequencies, a technique for improving accuracy by combining multifrequency SAR data has been studied. Our research is focused on dual-polarization data due to wide availability and low cost of this data.

Study area (Israel)
The South-Eastern Mediterranean region (Israel) represents a typical Mediterranean-to-arid transition zone from a sub-humid Mediterranean climate to a desert, where the mean annual rainfall ranges from above 600 mm to below 100 mm ( Figure 1) (Shoshany et al., 1996;Sternberg and Shoshany, 2001). Vegetation patterns consist of four life forms in this area: woodlands (mainly Pinus halepensis), taller than 5 m with average biomass of 8 kg/m 2 ; shrubs (dominant species are Pistacia lentiscus and Rhamnus palaestinus), with a hight range of 0.5 m -5 m and average biomass of 4.5 kg/m 2 ; dwarf shrubs (such as Sarcopoterium spinosum and Thymelaea hirsuta), with heights less than 0.5 m and an average biomass of 0.7 kg/m 2 ; herbaceous plants, with maximum biomass of 0.8 kg/m 2 (Svoray et al., 2001). These vegetation types vary due to precipitation and magnitude and frequency of human disturbances. According to the precipitation variation across the gradient zones and the intensity of human disturbance, the shrubland biomass varies approximately from 0.5 kg/m 2 to 5.0 kg/m 2 .

Figure 1. The study area with field experiments sites is marked
with a rectangular box. The contour indicates mean annual precipitations (mm/year). Table I shows the remote sensing data used in this study: dualpolarization L-band ALOS-PALSAR2, dual-polarization Cband Sentinel-1, and mean annual precipitation (MAP). The SRTM (1 arc-second) elevation data was applied to SAR data for terrain correction for PALSAR2 and Sentinel-1. Landsat 8 data were acquired for land classification. All the remote sensing data were co-registered in 10-meter spatial resolution (WGS84/UTM) by ENVI 5.5  Table 1. Remote sensing data.

Reference biomass
Numerous allometric models have been suggested to provide an adequate methodological solution for biomass estimation of shrubs and trees based on their geometric properties (Chave et al., 2014;Conti et al., 2019;He et al., 2018;Le Toan et al., 2011;Pereira et al., 1995;Sternberg and Shoshany, 2001). The average shrub biomass allocation of all species combined was 9.1~27.6% of leaf biomass and 72.4~90.9% of the woody biomass (Sternberg and Shoshany, 2001). Shrub stems (or trunks) and branches are therefore the dominant component of shrub' volume. Here we used the allometric model presented by Pereira (Pereira et al., 1995): ( 2) where h is the shrub's height and d is its diameter.
Forty-seven ROIs (regions of interest) of shrublands in the study area, 100m by 100m of each ROI, were selected representing local areas of low disturbance. Their biomass was calculated by combining high-resolution orthophoto images and field measurements (Chang and Shoshany, 2017). The average of tall tree biomass allocations of most species was determined mainly by 72.6% of trunks, followed by 23.8% of branches and 3.4% of foliage (Vargas-Larreta et al., 2017). For estimating tree biomass in this paper, we utilized the following model of He et al. (2018) using the trunk diameter as an independent variable: ( 3) where DBH is a diameter at breast height for a tree.
Ten ROIs with relatively homogeneous tree density in the study area were selected for utilizing this method. For each ROI, five (or six) trees were measured within 50m by 50m between February and April 2019.

MULTI-FREQUENCY RADAR VEGETATION INDEX
The microwave transmissivity (horizontal co-polarization) relationships with biomass are shown in Figure 2 (Entekhabi and Moghaddam, 2007;Ulaby et al., 1986). In this study, the Cband and L-band SAR data have been used as primary imagery source data for biomass estimation (Shimada et al., 2009;Torres et al., 2012). Multivariate linear regression analysis is frequently used by the remote sensing community for bio-physical estimation, using multi-frequency and multi-polarization SAR data (Berninger et al., 2018;Eisfelder et al., 2012;Lu, 2006b). However, it requires the reparameterization of regression coefficients for each site, thus it is very limited for wide areas mapping. The following generalized PRVI presented by Chang et al., (2018) is suggested here as an alternative: where a0 and a1 are constant coefficients, and B is the biomass in kg/m 2 .
While the PRVI was developed based on full-polarimetric data, the dual-polarization based PRVI can be obtained by replacing the DOP (average of DOPH and DOPV) parameters with one of the DOPH and DOPV. The average of DOP and DOPH are alternative to each other with little difference, and the r-square between the two is 0.952, as obtained from the previous study using the full polarimetric PALSAR data (Chang et al., 2018). Therefore, the PRVI for dual-polarization has the following form, under the assumption of the homogeneously distributed target (Chang et al., 2018;Ulaby et al., 1992Ulaby et al., , 1986): where the subscript pp denotes vv-or hh-polarization, the subscript pq denotes the cross-polarization, and the dualpolairzation combination would be VV+VH or HH+HV.
Consideration of the frequency-dependent transmissivity, a combined PRVI based biomass model using multiple frequency, can be expressed by the extended form of equation (4): where Tk indicates the transmissivity at frequency band k, PRVIk for the polarimetric radar vegetation index, and b0 and bk are constant coefficients.
The Tk can be replaced by simple equations derived from the relationships between biomass and transmissivity ( Figure 2) at P-, L-, C-, and X-bands (Ulaby et al., 1986): Equation (6) and (7) would offer better accuracy with containing the proportion of biomass due to each different penetration depth for each frequency, but they are too complex to be used in developing a biomass model.
The DOP parameter of PRVI is relevant with volume scattering (multiple scattering) than biomass itself (Chang and Shoshany, 2017;Chang et al., 2018). Therefore, equation (6) can be expressed as separating an integral form of DOP and cross-pol backscattering coefficients under the assumption that multiple scattering, according to the wavelengths, is similar. In the multiplication between the transmissivity T and the crosspolarized backscattering coefficient ( ), the transmissivity T can be ignored since the rate of the change in T with the increase in biomass is much smaller than the changes in ' ' with the increase of biomass.
Therefore, equation (7) can be simplified by a new parameter, multi-frequency polarimetric radar vegetation index (MPRVI) which represent multiplying the ensemble average of DOP and cross polarized backscattering coefficient: where the subscript Q is H or V for HH+HV or VV+VH dual polarization, respectively, and 〈 〉 denotes the ensemble average for multi-frequency bands.
Utilizing C-band Sentinel-1 and L-band PALSAR2 data, the MPRVI is calculated and assessed with reference to biomass from the previous research (Chang et al., 2018). Figure 3 shows that correlation between the MPRVI (equation 8) and the multiplying PRVI with the characteristic of transmissivity (equation 7) obtains high correlation (r 2 : 0.908). Such a high correlation indicates that the proposed MPRVI incorporates the characteristics of multi-frequencies while taking into account transmissivities. It can replace thus the complicated non-linear equation (7) with the simple form of equation (8) for biomass estimation. Therefore, the biomass model utilizing multifrequency SAR data can be expressed by the MPRVI, in a modified form of equation (4): where c0 and c1 are constant coefficients.

Figure 2.
Canopy transmissivity for sample microwave bands (horizontal co-polarization) corresponding to biomass (Ulaby et al., 1986), and structure parameters for woody plants (bottom left for tree and bottom right for shrub).

RESULTS AND DISCUSSION
For each ROI, the average of the values extracted for 5 by 5 pixels (50 m by 50 m) was used as radar parameters (backscatters, DOPs, and PRVIs) from the PALSAR and Sentinel-1 data. For a preliminary analysis, the correlations between radar parameters and biomass were analysed by using dual-polarization C-band (DC) Sentinel-1 data, and dualpolarization L-band (DL) PALSAR2 data: ( : 0.066), ( : 0.534), ( : 0.604), ( : 0.570), ( : 0.639), and ( : 0.573). The L-band cross polarized backscattering coefficient ( ) and C-band DOP ( ) showed relatively better correlations with biomass than other parameters. PRVI and MPRVI were calculated and analysed with reference biomass: ( : 0.613), ( : 0.656), and MPRVI ( : 0.739). The PRVI values showed better performance than other radar parameters in both C-band and L-band data. Above all, the MPRVI performed better than single frequency PRVI. These results show that synergetic use of multi-frequency could improve biomass estimation.

Model fitting
Two unknown coefficients (c0 and c1) of equation (9) were obtained by using the regression analysis with allometric biomass data based on the following logarithmic form: A simple color composition map was shown in Figure 4 (left) with an RGB image combination (R: MPRVI, G: NDVI, B: HHL-VVC), where NDVI was acquired from Landsat data. Green color indicates the vegetated areas and blue shows bare surface regions, and red indicates dense forest regions or urban and rocky regions with high backscattering coefficients. A MPRVI-based biomass map was generated by equation (10) as shown in Figure 4 (right). The MPRVI was still affected by background effects from rough bare surfaces and rocks, while it was found to be instrumental in reducing the effect of surface roughness compared with other radar parameters. The biomass extracted from the MPRVI based model ranges from 1.1 to 1.8 kg/m 2 for shrublands where the annual rainfall is above 400 mm/year, and from 0.4 to 1.4 kg/m 2 for shrubs in the semi-arid region with 200-400 mm/year. The average biomass for forest regions (mostly Pinus halepensis) exceeds 3.0 kg/m 2 . Southern Israel (Negev desert, reddish color in RGB image) is located mostly in the arid region having very low biomass; however, the biomass obtained in practice shows moderate biomass because the MPRVI was overestimated by high crosspolarized backscattering coefficients for very rough surfaces and rocks in the region. For these reasons, the biomass model was not applied in southern Israel (Negev desert) with NDVI values lower than 0.2.

Model validation
For validation, data set was partitioned randomly into two subsets of 50% for training and testing, using regression leaner toolbox of MATLAB. Figure 5 shows the response plot of the MPRVI based biomass model in log scale. The average of square error was 0.031 for 57 data points. In general, the error level was significant at the relatively low biomass region (semiarid to arid region). One of the reasons would be the high geodiversity of shrubs (high fragmentation and heterogeneity), and another reason could be the topography effect (the study area is hilly). Figure 6 shows the validation of the MPRVI based model with reference to biomass data. Since there were more measurements of shrubs than tall trees, many of the estimated biomass values were less than 2kg/m 2 . Overall, it shows a good agreement with reference biomass including shrubs and tall tree (r-square: 0.899, root-mean-squared error (RMSE): 0.381 kg/m 2 , and relative RMSE (RRMSE): 10.8%). Most other investigations of remote sensing biomass estimation have tested biomass ranges up to 600 Mg/ha, primarily in forests and dense woodlands, whereas semi-arid shrublands' biomass estimation and mapping is most challenging due to their relatively low biomass, ranging between 0 to 5 kg/m 2 (50 Mg/ha) (Eisfelder et al., 2012;Vaglio Laurin et al., 2021). Recently, advanced SAR technology such as PolInSAR, timeseries, and integration techniques have been applied to multifrequencies data to improve biomass estimation (Kraatz et al., 2021;ME and Kumar, 2021;Vaglio Laurin et al., 2021). Future BIOMASS mission (satellite-based P-band SAR data) which is planned to be launched by the European Space Agency in 2022 (Le Toan et al., 2011) would be instrumental for further testing and implementing our methodology.

CONCLUSION
Based on the concept of polarimetric radar vegetation index (PRVI), the multi-frequency polarimetric radar vegetation index (MPRVI) was developed and applied to the semi-arid Mediterranean ecosystem of Israel. Key radar parameters (i.e., the degree of polarization and the cross-polarized backscattering coefficient) for biomass estimation are integrated into the MPRVI as a simple form for biomass estimation by using dual polarization PALSAR2 (L-band) and Sentinel-1 (C-band) data. MPRVI based biomass model shows a high r-square value of 0.899 with reference to biomass. Our investigation may contribute to improving a wide range of biomass mapping. Further research is needed to apply this model to a wide range of biomass data. Integration techniques by using multi-sources such as optical sensors (multi-spectral), Lidar, Interferometry SAR (InSAR) may improve biomass estimation in a wide range of biomass (e.g. Chang and Shoshany, 2016;Cutler et al., 2012;Vaglio Laurin et al., 2021).