Improving Estimation Accuracy of Growing Stock by Multi-Frequency SAR and Multi-Spectral Data over Iran’s Heterogeneously-Structured Broadleaf Hyrcanian Forests

Abstract

Via providing various ecosystem services, the old-growth Hyrcanian forests play a crucial role in the environment and anthropogenic aspects of Iran and beyond. The amount of growing stock volume (GSV) is a forest biophysical parameter with great importance in issues like economy, environmental protection, and adaptation to climate change. Thus, accurate and unbiased estimation of GSV is also crucial to be pursued across the Hyrcanian. Our goal was to investigate the potential of ALOS-2 and Sentinel-1’s polarimetric features in combination with Sentinel-2 multi-spectral features for the GSV estimation in a portion of heterogeneously-structured and mountainous Hyrcanian forests. We used five different kernels by the support vector regression (nu-SVR) for the GSV estimation. Because each kernel differently models the parameters, we separately selected features for each kernel by a binary genetic algorithm (GA). We simultaneously optimized R² and RMSE in a suggested GA fitness function. We calculated R², RMSE to evaluate the models. We additionally calculated the standard deviation of validation metrics to estimate the model’s stability. Also for models over-fitting or under-fitting analysis, we used mean difference (MD) index. The results suggested the use of polynomial kernel as the final model. Despite multiple methodical challenges raised from the composition and structure of the study site, we conclude that the combined use of polarimetric features (both dual and full) with spectral bands and indices can improve the GSV estimation over mixed broadleaf forests. This was partially supported by the use of proposed evaluation criterion within the GA, which helped to avoid the curse of dimensionality for the applied SVR and lowest over estimation or under estimation.

1. Introduction

Hyrcanian forests are known as remnants of the Pleistocene era that survived the frost period [1]. These forests are located in regions of northern Iran and part of Caucasus, and embrace a high species and structural diversity of uneven-aged mountainous broadleaf forests distributed across a high altitudinal gradient [1,2]. Recently, portions of these forests were inscribed in the list of UNESCO World Natural Heritages [3]. The growing stock volume (GSV) is one of the important allometric biophysical forest attributes. It is closely related to other forest quantities such as height and aboveground biomass and is of great importance in the forest ecology, management, and carbon storage [4,5]. Tree-level biomass is conventionally derived by using species-specific allometric relations and wood density from ground-based measurements. However, the high cost, time, and the limited geographical coverage prohibitively challenge these methods. On the other hand, remote sensing data from spaceborne SAR and multispectral sensors with proper radiometric and spatial resolution and sufficient time intervals of data acquisition from the desired areas have been proven to provide important proxies in forestry research [6,7]. Due to mentioned historical and environmental reasons, development and implementation of remote sensing-assisted methods serve the overarching aim of monitoring and sustainable management of Hyrcanian forests.

Amongst the recent attempts for GSV estimation by state-of-the-art multispectral data, Chrysafis et al. (2017) [8] estimated the GSV by blending the Sentinel-2 and the Landsat data with the Random Forest (RF) model, and concluded that near-infrared and the red edge domains greatly affect GSV estimation. In addition, Mura et al. (2018) [9] estimated the GSV by the Sentinel-2, Landsat, and the Rapideye sensors and concluded that, beside the near-infrared and red edge regions, the SWIR region is also effective because of its sensitivity to the water content in the canopy. In both studies Sentinel-2 data were suggested to excess others in performance.

The radio detection and ranging (radar) sensors considerably contributed to solving the limitations of optical sensors including their inability to penetrate the canopy and less sensitivity to the vertical canopy structure. The ability of radar data to estimate biophysical forest characteristics is also less susceptible to weather conditions and acquisition time, which eases monitoring of mountainous forests that mostly occur in humid, cloudy, and foggy areas. However, changes in radar wavelength and type of polarization results in differences in both analytical workflow and the achieved estimation performance. Moreover, higher trunk volume leads to underestimation of actual GSV values due to the saturation in the scattering form dense canopy, which improves with an increase in wavelength [10,11,12]. Gao et al. (2018) [13] estimated GSV by the Dual polarized ALOS-1 data and reported a higher potential and later saturation of L-band HV cross-polarization channel than HH co-polarized channel. By using multi-temporal dual polarized ALOS-1 sensor data, Antropov et al. (2013) [14] concluded that the multi-temporal method was superior in prediction (with HH co-polarized channel performing better for mature trees), yet the saturation happens in the high GSV values. In addition, Chowdhury et al. (2014) [15] estimated the GSV by the ALOS-1 multi-temporal full polarimetric data, from which covariance and coherency matrices, as well as the phase difference between HH & VV channel and the coherency between HH and VV channels, were extracted. They showed that full polarimetry data has a high ability for GSV estimation.

Multi-sensor remote sensing approaches are highly capable for forest applications. In our study, we percept GSV estimation over mountainous broadleaf forests from a slightly different perspective. The biophysical characteristics of forest can be studied in multi-spectral approach by focusing on the biochemical aspects such as chlorophyll and in SAR approach on radar wave penetration in the canopy [16,17]. Mauya et al. (2019) [18] estimated GSV by the ALOS-2’s global mosaic, the Sentinel-1 and the Sentinel-2 sensors data. They concluded that using SAR data alone was unlikely to provide a good estimation ability for GSV, while a combined use of Sentinel-1 and Sentinel-2 data were advantageous. In Iran’s Hyrcanian region, recent investigations of the ability of remotely sensed data and methods include Vafaei et al. (2018) [19] who estimated biomass by the ALOS-2 full polarimetric and Sentinel-2A data. The Sentinel-2 returned the moderate accuracy, whereas the ALOS-2 individually led to minimum estimation accuracy.

The overall objective of this research was to estimate the GSV in Hyrcanian uneven-aged mountainous broadleaf forests based on leveraging a broad range of possibilities in optical and radar data processing. To this aim, we used ALOS-2 full polarimetric, Sentinel-1 dual polarimetric and Sentinel-2 multi-spectral data. We only concentrated on combined use of polarimetric and spectral features. In addition, we compared five different kernels in support vector regression (SVR) for GSV estimation. We additionally applied a heuristic feature selection by binary genetic algorithm, in which we simultaneously optimized the root mean square error (RMSE) and coefficient of determination (R²) for each kernel separately. The workflow and findings of this study are mainly significant due to multiple challenges associated with the structure and composition of our test site, including the severe topography, limited field samples, highly mixed tree species in various ages, and complex slope-aspect structure.

2. Materials and Methods

2.1. Study Area

The study area is located in Guilan province’s Nav-Asalem region in northern Iran (see Figure 1). The altitude varies between 100 to 2000 m above sea level, and the slope ranges from 0 to 73 degrees. The climate is temperate and cold, average annual rainfall is 1200 mm per year and the mean temperature is 12.4-degrees Celsius. The forest comprises uneven-aged mountainous broadleaf stands dominated by oriental Beech (Fagus orientalis Lipsky.) and hornbeam (Carpinus Betulus L.), accompanied by other broadleaf shrub and tree species.

2.2. Field Data Inventory

Field data were collected during July 2014 to March 2015. The in situ measurements were conducted in 148 circular plots with 18 m radius each and distributed in a randomly positioned square grid, as prescribed by the technical bureau of the Iranian Forests, Rangelands and Watershed Management Organization (FRWO) [20]. In each plot, the diameter at breast height (DBH; i.e., 1.5 m above ground surface) and the species type were recorded for all trees. For a sample of trees, the dominant species was recorded in each plot. Finally, the GSV values were calculated in each stand by means of DBH and lookup tables for each species separately.

Plot center locations were collected by global positioning system (GPS) in the WGS84 coordinate system. The original values of GSV in each plot were multiplied by 10 in order to be transformed from m³ per 0.1 hectares to m³ per hectare. The inventoried GSV ranged between 98.60 and 385.8 m³ ha⁻¹ (see

Table 1. Summary of descriptive statistics of GSV.

Descriptor	Value
Mean(m3 ha−1)	247.86
Minimum(m3 ha−1)	98.60
Maximum(m3 ha−1)	385.8
Standard deviation(m3 ha−1)	54.63
Number of plots	148

).

2.3. Remotely Sensed Data

2.3.1. Sentinel-2 Data

Sentinel-2 satellite carries a multi-spectral sensor in 13 bands from 400 nm to 2400 nm with spatial resolutions of 10, 20, and 60 m [21]. The 60-m resolution bands include B1 (430 nm), B9 (940 nm) and B10 (1340 nm) that are mainly used for atmospheric correction. The Sentinel-2 data in Level1C was acquired on 21 June 2016 and downloaded from European Space Agency (ESA) repository. After performing the atmospheric and geometric correction incorporating SRTM 1 aSec DEM [22], the spectral features were extracted from Sentinel-2 data. Since the acquired Level1C data represents the top-of-atmosphere (TOA) reflectance, we used Sen2Cor algorithm for atmospheric correction to calculate Level 2A data representing the bottom-of-atmosphere (BOA) reflectance. The sen2cor algorithm is an image-based correction that performs image correction based on lookup tables extracted from RadTran algorithm [23], with the main benefit that it does not need local meteorological data for correction. We only used the original bands and vegetation indices (see

Table 2. Sentinel-2 spectral features.

Feature	Used Bands	Resolution	Feature	Used Bands	Resolution
B2, B3, B4, B8	Original band	10 m	GEMI [25]	B8, B4	10 m
SAVI [26]	B8, B4	10 m	ARVI [27]	B2, B4, B8	10 m
TSAVI [28]	B8, B4	10 m	NDVI [8]	B4, B8	10 m
MSAVI [29]	B8, B4	10 m	B5, B6, B7, B8a B11,B12	Original bands	20 m
MSAVI2 [29]	B8, B4	10 m	NDI45 [30]	B4, B5	20 m
DVI [8]	B8, B4	10 m	MTCI [31]	B4, B5, B6	20 m
RVI [32]	B8, B4	10 m	MCARI [33]	B3, B4, B5	20 m
PVI [8]	B8, B4	10 m	REIP [34]	B4, B5, B6, B7	20 m
IPVI [35]	B8, B4	10 m	S2REP [34]	B4, B5, B6, B7	20 m
WDVI [36]	B8, B4	10 m	IRECI [34,37]	B4, B5, B6, B7	20 m
TNDVI [32]	B8, B4	10 m	PSSRa [38]	B4, B7	20 m
GNDVI [39]	B3, B8	10 m

B2 (490), B3 (560), B4 (665), B5 (705), B6 (740), B7 (783), B8 (842), B8a (865), B11 (1610), B12 (2100), Unit = nm.

) for GSV estimation. Due to the different spatial resolutions of the Sentinel-2 imagery, we used the 5 × 5 local window for bands and indices with 10 m spatial resolution, whereas 3 × 3 local window was used for bands and indices with 20 m spatial resolution. The mean values of data extracted in each local window over each field plot inventory were applied for GSV estimation. The entire process was performed in SNAP V6 software [24].

2.3.2. Sentinel-1 Data

Sentinel-1 is equipped with a synthetic aperture radar (SAR) antenna that scans the Earth in dual polarization mode in VV and VH channels in the C band. Our study area was scanned in dual polarization mode on 22nd of July 2017 and data was downloaded from the ESA repository. The data was scanned in single look complex (SLC) mode using the interferometric wide (IW) method. The dataset was primarily calibrated. Since the data was in dual polarimetry mode, the calibration was performed in the complex calibration manner to preserve the phase information. For extracting polarimetric information, we first extracted the C2 matrix (non-coherent covariance matrix) from the calibrated data. For dual polarimetric data, we used the Equation (1) to calculate the C2 matrix [40]. Data was scanned by IW method and was thus required to be debursted by the Sentinel-1 TOPS deburst method. For squared pixels, we applied the multi-looking operator with five looks in the azimuth direction and one look in the range direction to reach the 17-meter resolution, i.e., the spatial resolution close to the field sample plot size. For protecting polarimetric information, we applied the polarimetric refined Lee filter (7 × 7 local window) on the debursted C2 matrix to decrease the speckle noise [39]. The range-doppler terrain correction method with the SRTM 1Sec HGT DEM [22] was used to georeference the C2 matrix. Finally, the dual-polarimetric H-A-Alpha decomposition was applied on the georeferenced C2 matrix [40]. Consequently, we used a 3 × 3 local window due to the available 17 m spatial resolution to extract information on sample plot level. The entire analysis was performed in SNAP V6 software [24]. The extracted features are summarized in

Table 3. Sentinel-1 features.

Feature	Elements
C2 matrix	C11, C12_real, C12_image, C22
H/A/Alpha	Entropy, Anisotropy, Alpha

.

$$C_2=\left[\begin{array}{cc}\langle {\left|S_{VV}\right|}^2\rangle & \langle S_{VV}{S}_{VH}^{*}\rangle \\ \langle S_{VH}{S}_{VV}^{*}\rangle & \langle {\left|S_{VH}\right|}^2\rangle \end{array}\right]$$

(1)

2.3.3. ALOS-2 PalSar Data

The ALOS-2 satellite features a fully polarimetric synthetic aperture radar (SAR) antenna. This sensor scans the Earth in HH, HV, VH, and VV polarization channel in L band. We downloaded the data acquired on 10th of June 2015 over our study site from the repository provided by the Japan Aerospace Exploration Agency (JAXA). This data was scanned in SLC mode with Stripmap-2 method. The data preprocessing was performed solely for the polarimetric features. For preserving phase information, the raw data was initially calibrated in the complex calibration manner. The T3 coherency matrix was extracted from the calibrated data (Equation (2)) [41]. One may note the negative effects caused by the severe topography on scattering from the scatterers and consequently on the polarimetric information [42]. One crucial effect of topography is the rotation of radar wave in the line of sight because of the slope in the azimuth direction. This phenomenon is known as polarization orientation angle (POA) and is compensated by applying Equations (3)–(5) [41,42,43]. The multi-looking operator with six looks in azimuth direction and four looks in range direction was applied to the T3 matrix for the squared pixel and speckle noise reduction. Finally, we reached the 15-m spatial resolution that was close to those of field sample plots. As already mentioned for Sentinel-1 data processing, we used refined Lee speckle filter to simultaneously reduce speckle noise and preserve polarimetric information [41,44,45]. The filter was applied with 7 × 7 local window, 3 × 3 target window, and sigma of 0.9 to the T3 matrix. Following these steps, the T3 matrix was consequently georeferenced by SRTM 1Sec HGT DEM [22] using the range-doppler terrain correction method. In this research, we only extracted the polarimetric features (see

Table 4. ALOS-2 polarimeteric features.

Feature	Elements
H/A/Alpha [41]	Anisotropy, Alpha, Entropy, Beta, Delta, Gamma, Lambanda Alpha (1,2,3), Lambanda (1,2,3)
Yamaguchi [46]	Surface, Double, Volume, Helix
Van zyl [47]	Surface, Double, Volume
Cloude [41]	Surface, Double, Volume
Generalized Freeman-Durden [48]	Surface, Double, Volume
Touzi [49]	Psi, Tau, Phi, Alpha Psi (1,2,3),Tau (1,2,3), Phi (1,2,3), Alpha (1,2,3)
RVI [41]	RVI
SPAN [41]	SPAN
Pedestal height [41]	Pedestal height

). We used a 3 × 3 local window to extract the information on sample plot level. Similar to the preceding steps, the processes were performed in SNAP V6 [24].

$$S=\left[\begin{array}{cc}{S}_{HH}& S_{HV}\\ S_{VH}& S_{VV}\end{array}\right].\overrightarrow{k}=\frac{1}{\sqrt{2}}\left[\begin{array}{c}{S}_{HH}+S_{VV}\\ S_{HH}-S_{VV}\\ 2S_{HV}\end{array}\right] .T_3=\overrightarrow{k}\overrightarrow{k}{}^{*T}$$

(2)

$$\delta =\frac{1}{4}\left[{\tan }^{-1}\left(\frac{-4Re(\langle \left(S_{HH}-S_{VV}\right)S_{HV}^{*}\rangle )}{-\langle {\left|S_{HH}-S_{VV}\right|}^2\rangle +4\langle {\left|S_{HV}\right|}^2\rangle }\right)+\pi \right]$$

(3)

$$V=\frac{1}{2}\left[\begin{array}{ccc}1+\cos 2\delta & \sqrt{2}\sin 2\delta & 1-\cos 2\delta \\ -\sqrt{2}\sin 2\delta & 2\cos 2\delta & \sqrt{2}\sin 2\delta \\ 1-\cos 2\delta & -\sqrt{2}\sin 2\delta & 1+\cos 2\delta \end{array}\right]$$

(4)

$$T_3^{POA}=V{T}_3{V}^T$$

(5)

In the equations S is the stocke’s matrix, S_HH to S_VV are the four polarizations, k is the Pauli target vector for calculating T₃ matrix, and δ is the POA angle.

2.4. Modeling by Machine Learning

2.4.1. Support Vector Regression

Support vector machine (SVM) is known as powerful, flexible, and robust-to-noise machine learning method. The main goal of SVM is to minimize the sum of squares error (SSE) in Equation (6) during the model training procedure. The support vector regression (SVR) is a regression version of SVM [50,51]. There are two main methods for SVR, ε and υ regression. In ε regression, an insensitive tube defines by ε parameter, but this parameter is user-defined and the procedure affects model accuracy by only involving noisy observation in the big values of ε or the non-important observation data in the small values of ε, which results in non-generalized model that cannot model different behaviors of the phenomena [50,52,53]. In the υ regression method, the υ parameter defines the fraction of support vectors number involved in the modeling procedure [52]. In this method, the ε is automatically calculated in the algorithm and better controls support vectors. Therefore, both caveats of the ε method including challenge of the selection of proper ε value and its effect on accuracy could be resolved [52]. Therefore, we used the υ regression. The kernel based solution was used for accommodating the non-linear and complex behaviors. The used kernel passes the features into a new feature space [53], following which the model is trained in the new feature problem [51]. Whereas there are various kernels for SVR, the choice proper kernel for each problem depends on the nature of the problem [49]. We tested five different kernels (see

Table 5. Applied SVR kernels [54].

Kernel	Hyper-Parameter
Laplace	Sigma, Cost
RBF	Sigma, Cost
Polynomial	Degree, Offset, Scale, Cost
Sigmoid	Scale, Offset, Cost
Bessel	Sigma, Order, Degree, Cost

). We tuned the hyper-parameters by repeated cross-validation based on RMSE (see Section 2.4.3). Besides the kernel’s hyper-parameter, the cost parameter was also tuned with each kernel to maintain the model’s flexibility [51].

$$\min \left(\frac{1}{2}\Vert w{\Vert }^2+C\left(\nu \epsilon +\frac{1}{l}{\displaystyle \sum }_{i=1}^l\left({\xi }_i+{\xi }_i^{*}\right)\right)\right)$$

(6)

In Equation (6), C regularization constant, υ parameter between 0 and 1, ε value and ${\xi }_i^{*}$ is slack variable. We analyzed the results of three different experiments. In the first experiment, the GSV was modeled using a single sensor approach. In the second experiment, we modeled the GSV using all the features obtained through the multi-sensor approach. Finally, in the third experiment, we model the GSV using the GA selected features.

2.4.2. Feature Selection

Feature selection is an essential task in working with hyper-dimensional data, in which the best subset from all features is selected based on lower cost of model training and higher accuracy. This solution is essential in the presence of a large number of extracted features and small sample size to avoid the prohibitive effects of the curse of dimensionality on model performance and interpretability [55]. Therefore, a proper method for feature selection cares for the low training cost while avoiding the curse of dimensionality in the modeling [56]. The genetic algorithm (GA) is an evolutionary method for optimization. The method is inspired by natural selection in the real world, with the main idea of incorporating a wide range of suitable solutions that leads to selection of an optimum solution [57]. We used binary GA for feature selection. First, a population of chromosomes with 0 and 1 genes was randomly generated in which each chromosome is considered as a solution. The goodness of each solution was determined by the fitness function. The solutions with higher goodness had a higher chance to create a new solution. The new solution was then created by crossover procedure from two randomly selected solution with a high fitness value. The mutation operator changed the random genes with little probability and was used to maintain genetic diversity from one generation of a population of genetic algorithm chromosomes to the next. Finally, the algorithm was iterated until the fitness value remained unchanged [58]. Here, we selected the features for each kernel individually by cross-validation (see Section 2.4.3). The fitness function pursued the goal to simultaneously RMSE and R² values (see Figure 2). For this purpose, the fitness function as defined in Equations (7)–(9) was proposed. In this fitness function, the highest fitness value could potentially occur in R² equal to 1 and RMSE equal to 0, i.e., a fitness value equal to 190.

$$\alpha =90-{\tan }^{-1}\left(\frac{RMSE}{100\times R^2}\right)$$

(7)

$$d=\sqrt{RMS{E}^2+{\left(100\times R^2\right)}^2}$$

(8)

$$fitness=\alpha +d$$

(9)

2.4.3. Validation

The validation of model performance is the essential step in modeling [20]. Due to the most common case of limited sample size, using the k-fold cross-validation was necessary to avoid noisy and unstable predictions [18,50]. The RMSE, R² (Equations (10)–(11)) were used as error diagnostics based on their proven usefulness in assessing prediction of forest parameters [18,59]. In addition, the model stability in the presence of limited field sample size was checked by the repeated k-fold cross-validation. The model was assessed as more stable if the standard deviation (Equation (13)) of each of the diagnostic metrics was closer to zero [60]. Also we used MD (see Equation (12) for models over-estimation or under-estimation analysis. The higher positive values indicates the over-fitting and higher negative values indicates the under-fitting problem [61]. We used the 5-fold cross validation with 20 repetitions for more stable feature selection in the GA (Section 2.4.2), kernel’s hyper-parameter tuning (Section 2.4.1) and finally, the model validation. The statistical modeling and the feature selections were implemented in open source domain in R using the kernlab, caret, and GA libraries.

$$RMSE=\sqrt{{\displaystyle \sum }_{i=1}^n{\left(y_i^{obs}-y_i^{pred}\right)}^2/n}$$

(10)

$$R^2=\sqrt{{\displaystyle \sum }_{i=1}^n{\left(y_i^{pred}-\overline{y}\right)}^2/{\displaystyle \sum }_{i=1}^n{\left(y_i^{obs}-\overline{y}\right)}^2}$$

(11)

$$MD=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left(y_i^{pred}-y_i^{obs}\right)$$

(12)

$$SD=\sqrt{{\displaystyle \sum }_{i=1}^m{\left(P_i-\overline{P}\right)}^2/m}$$

(13)

In the Equations, the $y_i^{obs}$ is the observed GSV in plot i, $y_i^{pred}$ is the predicted GSV, n is the number of samples in fold, $\overline{y}$ is the GSV’s average in fold, $P_i$ is the validation metric in each fold, $\overline{P}$ is the metric’s average and m is the number by which the model was trained with different data.

3. Results

3.1. Single Sensor Models

3.1.1. Sentinel-2

Here, the 31 features extracted from the Sentinel-2 data were used for the GSV modeling (see Table 2). As previously mentioned in Section 2.3, the SVR model was employed using five different Kernels. The kernel’s hyper-parameters were tuned with repeated cross-validation. The results can be seen in

Table 6. Sentinel-2 modeling with 31 features.

Kernel	RMSE	R2	MD	RMSE_SD	R2_SD
Laplace	31.793	0.673	−0.003	2.919	0.071
RBF	31.644	0.673	−1.194	2.639	0.072
Sigmoid	31.693	0.674	−2.07	2.642	0.074
Polynomial	31.553	0.677	−0.62	2.630	0.071
Bessel	31.536	0.676	−1.38	2.662	0.072

.

As can be seen in Table 6, from the perspective of performance, the Bessel kernel with RMSE = 31.536 and R² =0.676 provided a better performance compared to other kernels. The polynomial kernel presented a very close performance to the Bessel Kernel. From the perspective of model stability (see Section 2.3.3), the Polynomial kernel with RMSE_SD = 2.630 and R²_SD = 0.071, provided more stable results than other kernels. The RBF kernel was the second more stable kernel (Table 6). The MD index shows that the all models are under estimate the GSV.

By considering both the models performance, stability and MD analysis the polynomial kernel exceeded other kernels, while all models generally showed close performance. Finally, no significant difference was observed amongst the kernels in the GSV modeling using the features extracted from Sentinel-2 (see Figure 3, Figure 4 and Figure 5).

3.1.2. ALOS-2

Here, 45 features were extracted from ALOS-2 data (see Table 4). The results are shown in

Table 7. ALOS-2 modeling with 45 features.

Kernel	RMSE	R2	MD	RMSE_SD	R2_SD
Laplace	46.959	0.306	6.93	6.508	0.124
RBF	46.588	0.318	6.99	6.510	0.123
Sigmoid	46.086	0.327	7.09	6.543	0.124
Polynomial	46.691	0.311	7.1	6.583	0.121
Bessel	46.576	0.316	6.97	6.516	0.122

.

As shown in Table 7, the sigmoid kernel with provided a better performance compared to other kernels by returning RMSE = 46.086 and R² = 0.327. This was closely followed by the RBF kernel with RMSE = 46.588 and R² = 0.318. Concerning the stability, the RBF kernel was more stable than other kernels (RMSE_SD = 6.510 and R²_SD = 0.123), followed by Bessel kernel. The Laplace kernel showed similar stability to that of Bessel kernel, yet a considerably lower performance compared with other kernels (see Table 7). The MD index shows that all of models significantly over estimate the GSV.

By considering both performance and stability and MD analysis, the RBF provided better results compared to the other kernels. However, all kernels performed roughly close to each other, with no significant difference between the kernels in the GSV modeling using the ALOS-2 features (see Figure 3, Figure 4 and Figure 5).

3.1.3. Sentinel-1

In this experiment, seven features, extracted from Sentinel-1 data solely were used for the GSV modeling (see Table 3). The results are show in

Table 8. Sentinel-1 modeling with seven features.

Kernel	RMSE	R2	MD	RMSE_SD	R2_SD
Laplace	46.815	0.299	4.536	6.360	0.117
RBF	47.093	0.294	3.232	6.349	0.111
Sigmoid	47.173	0.296	4.549	6.594	0.118
Polynomial	47.473	0.287	4.517	6.348	0.110
Bessel	46.879	0.297	4.596	6.302	0.113

.

As can be seen in Table 8, the Laplace kernel exceeded all other applied kernels with RMSE = 46.815 and R² = 0.299 again, closely followed by the Bessel kernel with RMSE = 46.879 and R² = 0.297. Moreover, the Bessel kernel provided more stable results than other kernels (RMSE_SD = 6.302 and R²_SD = 0.113). The MD index shows that all of models significantly over estimate the GSV after ALOS2. By considering both performance and stability, the Bessel kernel was slightly better compared to the other kernels.

3.2. Multi-Sensor Models

Here, we employed 83 features, obtained from Sentinel-2, ALOS-2, and Sentinel-1 data (see Table 2, Table 3 and Table 4). The results are shown in

Table 9. Performances achieved by modeling with 83 features.

Kernel	RMSE	R2	MD	RMSE_SD	R2_SD
Laplace	33.053	0.655	−0.538	3.431	0.071
RBF	32.701	0.654	−1.12	3.052	0.073
Sigmoid	32.636	0.659	−1.96	2.905	0.067
Polynomial	32.682	0.659	−2.19	2.884	0.070
Bessel	40.961	0.628	0.69	5.068	0.073

.

The sigmoid kernel performed best (RMSE = 32.636 and R² = 0.659), followed by polynomial kernel (RMSE = 32.682 and R² = 0.659). In addition the sigmoid kernel resulted in the most stable model with RMSE_SD = 2.905 and R²_SD = 0.067. The MD index shows that all of models except Bessel moderately under estimate the GSV.

As can be observed in Figure 3, Figure 4 and Figure 5, stacking the entire features into one feature set did not improve the results. This is mainly due to the joint effect caused by the curse of dimensionality and the limited sample size, and suggests the necessity of feature selection as an essential task.

3.3. Multi-Sensor Approach with Selected Feature

The

Table 10. Modeling optimized feature sets for each kernel.

Kernel	RMSE	R2	MD	RMSE_SD	R2_SD	Number of Features
Laplace	29.815	0.715	0.288	3.338	0.064	18
RBF	30.132	0.704	−0.069	2.595	0.060	19
Sigmoid	30.527	0.7	−0.43	2.903	0.058	16
Polynomial	30.034	0.708	−0.041	2.844	0.058	12
Bessel	30.707	0.696	−0.667	2.901	0.061	18

summarizes the results of models built with features optimized using the binary genetic algorithm as described in Section 2.4.2.

The Laplace kernel provided a superior predictive accuracy to other kernels. This was followed by the polynomial kernel. However, the RBF kernel proved to return the most stable predictions, followed by the polynomial kernel.

All in all, the highest accuracies for GSV prediction were returned by the polynomial kernel.

Table 11. Selected features for kernels by GA.

Kernel	Sensors	Selected Features by GA
Laplace	ALOS-2	Entropy, Anisotropy, Lambda2, Lambda3, Psi3, Alpha3
	Sentinel-1	Entropy, Anisotropy, Alpha
	Sentinel-2	B8, B5, B7, B8a, B11, B12, gndvi, pssra, tsavi
RBF	ALOS-2	Beta, Lambda, Freeman_dbl, Cloude_dbl, VanZyl_vol, Alpha2, Psi3, Alpha3
	Sentinel-1	Entropy
	Sentinel-2	B4, B5, B8a, B11, B12, gndvi, mcari, msavi2, pssra, tsavi
Sigmoid	ALOS-2	VanZyl_dbl, Psi1, Phi2, Psi3
	Sentinel-1	Entropy, Anisotropy
	Sentinel-2	B4, B5, B11, B12, gndvi, mcari, msavi2, pssra, s2rep, tsavi
Polynomial	ALOS-2	VanZyl_vol, Psi3, Alpha3
	Sentinel-1	Entropy
	Sentinel-2	B4, B5, B12, gemi, gndvi, mcari, pssra, tsavi
Bessel	ALOS-2	Entropy, Beta, VanZyl_dbl, Psi1, Phi1, Psi3, Phi3
	Sentinel-1	Entropy
	Sentinel-2	B4, B8, B5, B12, gemi, gndvi, mcari, pssra, s2rep, tsavi

summarizes the range of various features selected for the GSV modeling using each kernel. The selected features mainly include the ALOS-2 polarimetric decomposition, Sentinel-1 H-A-Alpha decomposition, and Sentinel-2 spectral features including NIR, red edge, SWIR bands and vegetation indices. The proper feature selection which resulted in the parsimonious feature sets enables achieving both predictive accuracy and model stability and robustness. Also, the MD index for all of the models compared to other experiments approach zero, especially for the polynomial kernel (MD = −0.041).

As summarized in Figure 3, Figure 4 and Figure 5 the performance of multi-sensor approach with selected features was slightly better than that achieved by Sentinel-2, which suggest that including selected polarimetry features along with the multi-spectral information improves the GSV across mountainous broadleaf stands. Also the over estimation or under estimation problem was solved in this experiment. For polynomial kernel the MD is equal to −0.041 and it is best MD index between all of experiments.

4. Discussion

The main objective of this research was to estimate the GSV in heterogeneously-structured and mountainous Hyrcanian forests in northern Iran. Predictive models and continuous monitoring of forest are especially vital in Iran’s Hyrcanian forests because of the current ongoing rate of degradation and their crucial role in Iran’s Hyrcanian forest ecosystem. Thus, investigating the potential of multi-frequency SAR, multi-spectral optical data for generating the reliable predictive models are essential in this study area. In this research, the GSV modeling was carried out using the polarimetric and multi-spectral features and their combination using three different approaches. The SVR with five various kernels was used for modeling because of non-linear and complex relation of features with the GSV particularly in multi-sensor approach. Also, there were several challenges in this procedure including the severe topography, limited field samples, highly mixed tree species in various ages, and complex slope-aspect structure.

In the first experiment, the GSV was modeled using the Sentinel-2 multi-spectral features. B11 and B12 bands from SWIR spectral region showed a high contribution because of their sensitivity to the water content in the canopy, that confirmed in research by Chrysafis et al. (2017) [8] and Mura et al. (2018) [9]. We also found the NIR and the red edge bands as good features in modeling which is due to their sensitivity to the chlorophyll and pigments of tree leaf. Relevant works of Chrysafis et al. (2017) [8], Mura et al. (2018) [9], and Chrysafis et al. (2019) [62] also described these features as influential features for the GSV estimation in their research. Basically, the vegetation indices that use the NIR and red-edge spectral bands have an effective contribution in the GSV estimation. This is mainly because of that the VI uses the combination of spectral bands and reaches the information that cannot be extracted from the single spectral band, as Chrysafis et al. (2017) [8] found this in their research. The Mura et al. (2018) [9] and Chrysafis et al. (2019) [62] also found that VIs have an effective contribution to GSV estimation.

The performed analysis indicated that the GSV modeling with ALOS-2 polarimetric features produced non-satisfactory results. However, increasing the number of polarimetric features can improve the results. The two main challenges for GSV modeling with ALOS-2 data in this study area are the harsh topography and the limited number of field inventory samples. As shown in Figure 6, the study area embrace harsh mountainous forest stands with high elevational gradient and complex slope-aspect structure. In this type of topography, phenomena like shadow and layover occur in SAR data. Particularly, no scanned data from the study areas were found for shadowing. Nevertheless, the radiometric values were significantly affected in layover behavior. Another effect is the POA (see Section 2.3.3) which was compensated in the preprocessing steps and thus results in no notable problem for the GSV modeling [43]. Yet the effect of slope in range direction, known as angular variation effect (AVE), still remains. In areas with forest cover and longer radar wavelength, the AVE introduces a notable effect. As such the double scattering mechanism increases in the slopes that face the radar incoming wave, whereas it decreases in opposite slopes. In addition, the volume scattering mechanism increases since the severe topography causes an increase in cross-pol scattering. In a severe mountainous area like our research site, the proposed correction methods do not completely refine the results [42,63]. In addition, a further effect is the effective scattering area due to the non-homomorphic imagining of radar. This effect can be maximized in the mountainous areas and thus dramatically affect the geo-referencing of the applied SAR data. This effect causes the higher radiometric values in front slopes facing the incoming radar wave and lower radiometric values in opposite slopes [42,43,63].

The second challenge for modeling GSV is caused by limited number of field inventory samples that is commonly the case across heterogeneously-structured uneven-aged Hyrcanian forests. Apart from statistical and model-related problems, this particularly hindered us from building species-specific GSV models and thus considering differences in scattering from tree species or different age classes. The collective effects of the above mentioned issues led to the rather mediocre results obtained by ALOS-2 data in our study (see Figure 3, Figure 4 and Figure 5). The H-A-Alpha decomposition features have an average ability in the GSV estimation. Vafaei et al. (2018) [19] found that these features have an average ability for biomass estimation in Hyrcanian forests too. The Touzi decomposition’s features are also effective, Sharifi et al. (2015) [64] found that the Tozi decomposition describes the target asymmetry in the forested area and can use for biomass estimation in Hyrcanian forest. The Van Zyle, Cloude, and Freeman decomposition features based on volume or double mechanism are effective in GSV modeling because of L-band wavelength penetration in the canopy and interaction with the tree’s stem and ground. So they are reasonable features for the GSV modeling as Kumar et al. (2012) [65] found this behavior. The achieved performance and model stability are comparable to similar research. Compared with the only available relevant case study from the area, Vafaei et al. (2018) [19] used limited polarimetry features. Unlike their study, we extracted a large number of polarimetric features in our research. Especially comparing to Vafaei et al. (2018) [19] with ALOS-2 in the same region, his conclusion is that ALOS-2 has a very weak ability for biomass estimation (R² = 16%). Since the biomass and GSV have near relation [5], we greatly improved the prediction ability for GSV with different target decomposition features (R² = 32%).

The performed analysis indicated that the modeling with Sentinel-1 polarimetric features produced non-satisfactory results for GSV modeling. The Sentinel-1 features are not frequently used in forestry research because of the shorter wavelength (C-band) of the data. In this wavelength, the scattering only occurs in the upper part of the canopy. In our study, we applied the dual polarimetric H-A-Alpha decomposition. Our result has better performance (R² = 29%) than Mauya et al. (2019) [18] who only used the scattering coefficient in GRD mode and achieved R² = 18%. Therefore, we can conclude that using the Sentinel-1 dual polarimetry H-A-Alpha decomposition is useful for GSV estimation. As can be seen in Figure 3, Figure 4 and Figure 5, the results obtained by the Sentinel-1 data were inferior than those of Sentinel-2 and ALOS-2 data. In addition to the underlying reasons as referred earlier, one may also note the lower penetration rate of the C-band Sentinel-1 compared with the L-band ALOS-2. This will result in a lower sensitivity to biophysical parameters such as GSV.

Finally, the performed analysis on the multi-sensor approach indicated that the GSV modeling with all of the features did not enhance the results compared to the modeling by Sentinel-2. This is mainly due to the curse of dimensionality effect and the limited sample size. To deal with this problem, we modeled the GSV by multi-sensor approach using the selected features. The results indicated that the modeling with the selected features provide better results than the modeling by Sentinel-2. The features were selected for each kernel separately proper feature selection procedure provides a good and stable GSV modeling which is not suffering from the curse of dimensionality anymore. Also, the over-estimation or under-estimation was reduced greatly by our feature selection method (see Figure 5). As a result, the simultaneous integration of multi-spectral and radar data features produce satisfactory and stable models for GSV modeling compare to the Mauya et al. (2019) [18], Vafaei et al. (2018) [19], and Sharifi et al. (2015) [64] works. Due to the general novelty of in-depth and state-of-the-art remote sensing analysis over Hyrcanian forests. In future research, issues such as the effects of different topographic corrections for polarimetry data, forest species type mapping, novel SVR kernels, application of feature fusion, and texture features will be explored.

5. Conclusions

This study generated empirical evidence on the use of ALOS-2, Sentinel-2, and Sentinel-1remotely sensing data for GSV estimation in Iran’s Hyrcanian forests on a landscape scale. We used nu-SVR with five kernels for non-linear behavior modeling. The result showed that in single sensor approach only the Sentinel-2 returned better results, while shortcomings were observed when applying ALOS-2 and Sentinel-1 across heterogeneously-structured and mountainous Hyrcanian forests with limited field samples. In multi-sensor approach, the curse of dimensionality caused inferior results compared with the case of Sentinel-2. In addition, we tested selecting features in multi-sensor approach individually for each kernel by GA based on simultaneous optimization of diagnostic measures in a proposed fitness equation. This contributed to both higher performance and more model stability with the lowest underestimation.