Forest Growing Stock Volume Estimation in Subtropical Mountain Areas Using PALSAR-2 L-Band PolSAR Data

Abstract

Forest growing stock volume (GSV) extraction using synthetic aperture radar (SAR) images has been widely used in climate change research. However, the relationships between forest GSV and polarimetric SAR (PolSAR) data in the mountain region of central China remain unknown. Moreover, it is challenging to estimate GSV due to the complex topography of the region. In this paper, we estimated the forest GSV from advanced land observing satellite-2 (ALOS-2) phased array-type L-band synthetic aperture radar (PALSAR-2) full polarimetric SAR data based on ground truth data collected in Youxian County, Central China in 2016. An integrated three-stage (polarization orientation angle, POA; effective scattering area, ESA; and angular variation effect, AVE) correction method was used to reduce the negative impact of topography on the backscatter coefficient. In the AVE correction stage, a strategy for fine terrain correction was attempted to obtain the optimum correction parameters for different polarization channels. The elements on the diagonal of covariance matrix were used to develop forest GSV prediction models through five single-variable models and a multi-variable model. The results showed that the integrated three-stage terrain correction reduced the negative influence of topography and improved the sensitivity between the forest GSV and backscatter coefficients. In the three stages, the POA compensation was limited in its ability to reduce the impact of complex terrain, the ESA correction was more effective in low-local incidence angles area than high-local incidence angles, and the effect of the AVE correction was opposite to the ESA correction. The data acquired on 14 July 2016 was most suitable for GSV estimation in this study area due to its correlation with GSV, which was the strongest at HH, HV, and VV polarizations. The correlation coefficient values were 0.489, 0.643, and 0.473, respectively, which were improved by 0.363, 0.373, and 0.366 in comparison to before terrain correction. In the five single-variable models, the fitting performance of the Water-Cloud analysis model was the best, and the correlation coefficient R² value was 0.612. The constructed multi-variable model produced a better inversion result, with a root mean square error (RMSE) of 70.965 m³/ha, which was improved by 22.08% in comparison to the single-variable models. Finally, the space distribution map of forest GSV was established using the multi-variable model. The range of estimated forest GSV was 0 to 450 m³/ha, and the mean value was 135.759 m³/ha. The study expands the application potential of PolSAR data in complex topographic areas; thus, it is helpful and valuable for the estimation of large-scale forest parameters.

1. Introduction

Forest carbon stocks are essential to our understanding of global climate change, and can be represented through extracting forest parameters [1]. The forest growing stock volume (GSV) is a key forest variable in the context of forest management and monitoring. Also, the forest GSV is referred to as the total volume (m³/ha) of the boles or stems of all living trees, and can be converted into above-ground biomass (AGB) by its density factor [2]. Therefore, the accurate quantification of forest biomass or GSV is essential for understanding the spatial distribution of carbon in vegetation areas, which can also provide effective predictions for the change trend of carbon stock [3]. In particular, large-scale forest GSV retrieval has become a research hotspot in recent years.

At present, many methods have been reported for estimating forest GSV. Traditional forest inventory approaches rely upon ground surveys by manually collecting the forest parameters of a single tree at sample plots [4]. However, the large amount of time required, labor intensity, and cost limit its application on a larger scale [5]. Remote sensing technology provides a possible solution to overcome such limitations, in particular the spaceborne remote sensing technique, which plays an important role in forest monitoring and management [6]. Optical remote sensing datasets (e.g., Landsat Thematic Mapper (TM) and Moderate Resolution Imaging Spectroradiometer, MODIS) can be used to estimate forest GSV [7], mainly by analyzing the relationship between forest parameters and vegetation indices (e.g., enhanced vegetation index (EVI), normalized difference vegetation index (NDVI) and perpendicular vegetation index (PVI)) [8,9,10,11]. However, the retrieved GSV values using optical remote sensing data are usually troubled with saturation effects, especially in the high carbon stock forests [12]. Another problem is the impact of cloud cover on image collection, constraining its application to moist regions (e.g., the tropical region) [13]. Light detection and ranging (LiDAR) data provides high accurate forest parameters for GSV estimation [14,15,16]. However, due to space discontinuous problems and complex data processing, LiDAR-derived GSV estimates usually can only be obtained over limited areas [17]. HHHSynthetic aperture radar (SAR) enables imaging in all-weather conditions and with continuous temporal coverage [18]. Now, it has been successfully applied in various fields [19,20,21,22]. Especially, the long-wavelength SAR has a wide potential in forestry applications [23,24,25]. Currently, the SAR techniques that have been utilized for the retrieval of forest parameters mainly are polarimetric SAR (PolSAR) [23,26], interferometric SAR (InSAR) [27,28], polarimetric interferometric SAR (PolInSAR) [29], polarization coherence tomography (PCT) [30,31,32], and tomography SAR (TomoSAR) [33,34,35]. In this paper, the full polarimetric SAR technique will be further used for retrieving forest GSV in subtropical mountain areas.

Radar polarimetry is the technique of acquiring, processing, and analyzing the polarization state of an electromagnetic field [36]. Forest characteristic information about the geometrical structure and geophysical properties can be obtained by analyzing polarimetric SAR signatures [37]. In an earlier study, the relationship between polarimetric signatures and forest GSV or AGB was studied by using high-frequency SAR data (e.g., X and C-band). Due to the low penetration, the short wavelengths interacted primarily with the forest canopy, and were suitable for low-carbon stock areas [38]. Lower frequency SAR data have stronger penetrating capability and can interact with various components of vegetation, and have been discovered to be more preferable than higher frequencies in high-carbon stock forests [39,40]. The phased array-type L-band synthetic aperture radar (PALSAR-2) can capture images in quad polarization modes, which provides an opportunity to study forest parameters using multi-polarization and multi-temporal data [41,42,43].

The L-band backscatter coefficient is sensitive to the biophysical parameters of forest [44,45,46]. However, the sensitivity is affected by many factors (e.g., radar polarization, forest structure, environment conditions, and topography) [39]. Among these factors, the complex terrain conditions can affect full polarimetric SAR data regarding both azimuth and distance, which is mainly reflected in the following three aspects. (1) The azimuthal slope causes a change in the polarization state, which leads to polarization orientation angles (POA) offsetting [47]. (2) Local terrain undulations cause a change in the effective scattering area (ESA), which leads to the change of actual backscatter [48]. (3) In vegetation-covered areas, the local terrain causes variation in penetration depth and scattering mechanisms, which are reflected in the angular variation effect (AVE) [49]. Zhao et al. [50] showed that the correction of these three aspects (POA compensation, ESA correction, and AVE correction) could reduce the topographic effect. However, in AVE correction, the critical correction factor n is obtained only according to the impact of the entire forest area, without considering the impact of different forest cover types. Meanwhile, the range of n values in different polarization channels in subtropical mountain areas need to be further explored.

The main purposes of this study are to (1) understand the role of terrain correction in forest GSV estimation; and (2) investigate the potential of PALSAR-2 L-band full polarimetric data for the retrieval of forest GSV in the subtropical mountain regions.

2. Materials

2.1. Study Area

The work was carried out in the Youxian county in Hunan of central China (27°05′ to 27°24′ N, 113°35′ to 113°55′ E, see Figure 1). It is a field site ground for forest research at Central South University of Forestry and Technology. The topography varies between 60–1386 m. The slope ranges from 0° to 84°. The climate type is a subtropical monsoon humid climate. The annual mean temperature is 17.8 °C. The average annual rainfall is 1410.8 mm, and most of rainfall occurs in the summer. The dominating forest type is coniferous forest, including fir and pine. In addition, there are some other vegetation types, such as bamboo and camphorwood.

2.2. Field Inventory Data

Field data collection was conducted from June to July 2016, with the help of the Central South University of Forestry and Technology, and the Chinese Academy of Forestry. A total of 60 forest plots with the size of 30 m × 30 m were surveyed for the experiment (Figure 1). The center of each plot was located by using a global positioning satellite (GPS) receiver, and the location (latitude and longitude) of the central point was recorded. These plots were independent from one another to avoid the spatial autocorrelation. Within each plot, the diameter at 1.3 m above the ground of each individual tree was measured by the diameter at breast height (DBH) ruler, and the tree height measurement was performed by the laser altimeter. The GSV was calculated using the method presented by Fang et al. [51,52], and the GSV of 60 plots ranged from 6.88 m³/ha to 434.42 m³/ha, with an average value of 194.75 m³/ha (

Table 1. Main biophysical properties of 60 plots in the study area. DBH: diameter at breast height.

	Range	Mean
DBH	4.06 to 30.10 cm	17.84 cm
Height	4.60 to 20.20 m	13.24 m
Number of Stems	30 to 350	96
Growing Stock Volume	6.88 to 434.42 m3/ha	194.75 m3/ha

). By random sampling, the plots were divided for the training (n = 44) and validation (n = 16) of models into two groups.

2.3. Polarimetric SAR Data and Pre-Processing

Full polarimetric (HH, HV, VH, and VV polarizations) L-band SAR data over this experiment site were acquired by the Japanese Aerospace Exploration Agency (JAXA) using the PALSAR-2. A total of seven scenes data were ordered as L1.1 level with the single-look complex (SLC) format in slant range geometry, and were acquired from June to October 2016 at approximately 4:22 local time. The central location of these datasets is approximately 27.18° N–113.68° E, and the dimensions are 69 km in azimuth and 25.8 km in range. The incidence angle ranges from 37.8° to 40.1°. The azimuth resolution is 2.97 m, and the range resolution is 2.86 m.

The basic data pre-processing steps were applied to reduce the geometric and radiometric distortions and speckle effects. The radiometric calibration of these data was first performed [53]. The coherency matrix [T₃] was generated and converted into the covariance matrix [C₃], which could represent the full polarimetric data. Then, these data were multi-looked with 7 × 10 in the azimuth and range directions. A Lee filter with a 3 × 3 window was applied to reduce the speckle effects. Geocoding was performed using shuttle radar topography mission (SRTM) elevation data (30-m spatial resolution). Finally, the SAR images were re-sampled to 30-m spatial resolution. PolSARpro software (Version 5.1.3, European Space Agency, Paris, France) was used to pre-process the SAR data. Gamma software was used to perform geocoding and resampling.

2.4. Ancillary Data

The ancillary data used in this study mainly include the SRTM digital elevation model (DEM) (https://earthexplorer.usgs.gov/) and land-use data product (http://www.dsac.cn/). The SRTM DEM (Figure 2a) has a 30-m resolution and was created by the National Geospatial Intelligence Agency and Jet Propulsion Laboratory. We used it to assist the SAR dataset geocoding. Besides, based on the SAR imaging geometry, terrain correction factors (i.e., projection angle and local incidence angle) could also be obtained from the DEM data. The land-use classification data (Figure 2b) with a spatial resolution of 30 m was provided by the Geographical Information Monitoring Cloud Platform at the same time as the PALSAR-2 dataset acquisition. According to the secondary classification of land use, the forest was divided into four types: woodland, shrubbery, sparse woodland (S-Woodland), and other forest (O-Forest), which would be used to assist the terrain correction factors.

3. Methodology

The processing flow chart presented in Figure 3 illustrates the framework of analysis steps. We firstly carried out basic pre-processing for the SLC level 1.1 datasets, including radiometric calibration, multi-looking, filtering, and geocoding. Secondly, POA and ESA correction were performed. Then, we calculated the values of n for different polarization channels by analyzing the correlation coefficients between the local incidence angles and backscatter coefficients. Thirdly, AVE correction was performed for different forest cover types, and then the results were spliced for later analysis. Fourthly, the correlations were analyzed between multi-temporal PolSAR data backscatter and forest GSV of all the sample plots, and the optimal data was selected for GSV estimation. Finally, the estimation models were constructed and compared by using the primary diagonal elements of the covariance matrix. Then, we estimated and mapped the GSV for the whole experiment area.

3.1. Terrain Correction

3.1.1. Polarization Orientation Angle Correction

The azimuth slope was the main factor that caused the polarization ellipse to rotate, and then affected the polarization state of the electromagnetic wave. To compensate for the impacts of the azimuth slope, the polarization orientation angles could be obtained by the circular polarization algorithm [54], as shown in Equation (1):

$$\eta =\frac{1}{4}\left[\mathit{arctan}\left(\frac{2\mathit{Re}\left[T_{23}\right]}{T_{22}-T_{33}}\right)+\pi \right]$$

(1)

where $\eta$ is the shift angle, and $T_{22}$, $T_{23}$, $T_{33}$ are the corresponding elements of matrix [T]. After acquiring the shift angle, a new rotated polarimetric covariance matrix ($C_{POA}$) can be formed by Equation (2):

$$C_{POA}=\left[U_{3(\eta )}\right]\left[C_3\right]{\left[U_{3(\eta )}\right]}^T$$

(2)

where $C_3$ denotes a polarimetric covariance matrix that represents multi-looked PolSAR data, and $U_{3(\eta )}$ is a rotation matrix.

$$\left[U_{3(\eta )}\right]=\frac{1}{2}\left[\begin{array}{ccc}1+\mathit{cos}2\eta & \sqrt{2}\mathit{sin}2\eta & 1-\mathit{cos}2\eta \\ -\sqrt{2}\mathit{sin}2\eta & 2\mathit{cos}2\eta & \sqrt{2}\mathit{sin}2\eta \\ 1-\mathit{cos}2\eta & -\sqrt{2}\mathit{sin}2\eta & 1+\mathit{cos}2\eta \end{array}\right]$$

(3)

3.1.2. Effective Scattering Area Correction

The ratio between the radar cross-section and reference area was usually used to express the backscatter coefficient [48]. The theoretical reference area was the pixel area, which did not change with topographic fluctuation. However, in most practical applications, the reference area was defined to be ground area (i.e., the effective area) that was affected by topographic fluctuation. The relation between the effective area and the theoretical reference area is shown as Equation (4) [55]:

$$A_{\sigma }=A_{\beta }/\mathit{cos}\phi$$

(4)

where $A_{\beta }$ and $A_{\sigma }$ represent the theoretical reference area and the effective area, respectively. $\phi$ is the projection angle, and $\mathit{cos}\phi$ is the correction factor for this step.

Then, a general correction equation for ${\sigma }_0$ could be obtained, which was the product of ${\beta }_0$ and $\mathit{cos}\phi$. For full PolSAR data, we corrected each element in the polarimetric covariance matrix using the same correction factor, and the equation could be written as follows:

$$C_{ESA}=\left[C_3\right]\times \mathit{cos}\phi$$

(5)

Here, $C_3$ is the polarimetric covariance matrix that has been POA corrected and geocoded. The projection angle $\phi$ is complementary to the smallest angle between the surface normal and the image plane, and can be obtained by DEM and orbit information.

3.1.3. Angular Variation Effect Correction

Since the local scattering mechanisms within the forest structure vary with the local incidence angles, further AVE correction was needed after the ESA correction. A simple cosine model was derived to reduce the angular effect [56]:

$${\sigma }_{corr}^0\left({\theta }_{loc}\right)={\sigma }_0\times {\left(\frac{\mathit{cos}{\theta }_{ref}}{\mathit{cos}{\theta }_{loc}}\right)}^n$$

(6)

where ${\sigma }_{corr}^0$ represents the terrain corrected backscatter, ${\theta }_{loc}$ denotes the local incidence angle, ${\sigma }_0$ denotes the uncorrected radar backscatter coefficient, ${\theta }_{ref}$ is the radar incidence angle, and n is a parameter that needs further discussion.

In a similar way, the polarimetric covariance matrix can be corrected by a $3\times 3$ correction coefficient matrix $\left[K_3\right]$ [50], and the expression is given as:

$$C_{AVE}=\left[C_3\right]\times \left[K_3\right]$$

(7)

where $C_3$ is the polarimetric covariance matrix that has been corrected by the previous two steps, and $K_3$ is the correction coefficient matrix.

According to Equation (6), the local scattering mechanisms are mainly affected by the local incidence angles. Therefore, the optimal correction factor n of different polarizations can be obtained through calculating and evaluating the correlation results between the local incidence angle and the corrected backscatter coefficient, which is:

$$n_{p,q}(z)=\mathit{arg}\mathit{min}\left\{\left|\rho \left({\theta }_{loc},C_{i,j}\right)\right|\right\}$$

(8)

where $\rho$ denotes the correlation between two parameters, C_i,j represents the elements of the corrected covariance matrix, z represents the different types of forest cover, and p and q represent different polarization channels. According to Equation (7), only the values of n corresponding to the primary diagonal elements of the covariance matrix need to be obtained. Considering the impact of forest characteristics on terrain correction, we try to calculate n corresponding to different types of forests in this paper in order to effectively reduce the topographic effect. The initial ranges of the n values are from zero to two. In addition, considering the computational complexity and accuracy of the n value, we set the interval of n to 0.01. The optimal n is determined by the absolute value of the correlation.

3.2. Retrieval of GSV

We performed terrain correction on all the full PolSAR data, and then selected the most relevant data for forest GSV estimation through time series analysis. The elements on the diagonal of the covariance matrix (corresponding to backscattering intensity of different polarization channels) were used as variables of estimation models. A few studies reported that individual backscattering measurements could be used to extract forest biological parameters [57,58]. Therefore, five single-variable models (Equations (9) to (13)) were first fitted to analyze the relationship between the single variables and the GSV. Among the five models, model (e) was derived from the parameterization of the improved Water-Cloud model, which we named the Water-Cloud analysis model.

(a) Linear function:

$${\sigma }^0={\beta }_1+{\beta }_{2}GSV$$

(9)

(b) Logarithmic function:

$${\sigma }^0={\beta }_1+{\beta }_2\mathit{ln}\left(GSV\right)$$

(10)

(c) Quadratic function:

$${\sigma }^0={\beta }_1+{\beta }_2\mathit{ln}\left(GSV\right)+{\beta }_3{\left(\mathit{ln}\left(GSV\right)\right)}^2$$

(11)

(d) Exponential function:

$${\sigma }^0={\beta }_1+{\beta }_{2}sqrt\left(GSV\right)$$

(12)

(e) Water-Cloud analysis function:

$${\sigma }^0={\beta }_1+{\beta }_2{e}^{\left({\beta }_{3}GSV\right)}$$

(13)

In addition, we also constructed a multi-variable regression model using three elements (HH, HV, and VV backscatter) on the diagonal of the covariance matrix, and compared it with the above five models to find a suitable model for forest GSV estimation and mapping.

$$\mathit{ln}\left(GSV\right)=a+b_1{\sigma }_{HH}^0+b_2{\left({\sigma }_{HH}^0\right)}^2+c_1{\sigma }_{HV}^0+c_2{\left({\sigma }_{HV}^0\right)}^2+d_1{\sigma }_{VV}^0+d_2{\left({\sigma }_{VV}^0\right)}^2$$

(14)

4. Results

4.1. Acquisition of Terrain Correction Factors

Before implementing terrain correction, the correction factors of each correction stage should be obtained. These correction factors could be divided into two categories: angular factors and parameter n, where the angular factors include the POA shift angle, projection angle, incidence angle, and local incidence angle. All of the angular factors are shown in Figure 4. Here, it is worth noting that the projection angle and local incidence angle of a single pixel are not complementary, especially in the terrain undulating regions. In order to show the correction parameters and correction effects, we chose one scene of data as an example to display the results. Here, the data acquired on 14 July 2016 was randomly selected.

Based on Equation (8), the distribution of correlation coefficients at different polarization channels can be obtained with a 0.01 interval. As shown in Figure 5, the different polarization channels are labeled with solid lines in different colors: HH polarization in red, HV polarization in blue, and VV polarization in green. The black dotted lines represent the position corresponding to the optimal n value. In order to effectively reduce the topographic effect, we have obtained the distribution of the correlation coefficients for different forest cover types, i.e., woodland (Figure 5a), shrubbery (Figure 5b), sparse woodland (Figure 5c), and other forest (Figure 5d).

From these figures, we can see that the variation trend of the correlation coefficients of different polarization channels is consistent for different forest cover types, increasing with the increase of parameter n values. After obtaining the correlation coefficient distribution, we can easily extract the optimum values of n by using Equation (8) in all four forest cover types.

In addition to the data acquired on 14 July 2016, the optimum n values of the remaining six scene data (16 June 2016, 30 June 2016, 25 August 2016, 22 September 2016, and 6 October 2016) are also calculated, and the results are shown in

Table 2. Results of the optimum n values. S-Woodland: sparse woodland, O-Forest: other forest.

	Data	Woodland	Shrubbery	S-Woodland	O-Forest
HH	16 June 2016	1.31	1.27	1.60	1.48
	30 June 2016	1.20	1.41	1.50	1.37
	14 July 2016	1.24	1.23	1.52	1.42
	11 August 2016	1.11	1.09	1.30	1.22
	25 August 2016	1.13	1.06	1.36	1.25
	22 September 2016	1.16	1.12	1.37	1.24
	6 October 2016	1.32	1.35	1.58	1.46
HV	16 June 2016	0.91	0.83	0.93	0.77
	30 June 2016	0.76	0.67	0.79	0.57
	14 July 2016	0.82	0.74	0.82	0.63
	11 August 2016	0.74	0.65	0.65	0.48
	25 August 2016	0.74	0.63	0.67	0.47
	22 September 2016	0.70	0.56	0.65	0.44
	6 October 2016	0.85	0.79	0.83	0.64
VV	16 June 2016	1.14	1.24	1.44	1.38
	30 June 2016	1.04	1.14	1.36	1.26
	14 July 2016	1.09	1.20	1.39	1.29
	11 August 2016	1.01	1.11	1.20	1.14
	25 August 2016	1.01	1.06	1.22	1.11
	22 September 2016	1.01	1.10	1.23	1.13
	6 October 2016	1.13	1.29	1.42	1.34

. From Table 2, it can be seen that the optimum n values of HV polarization is within the range of zero to one, while the HH and VV polarizations are greater than one.

4.2. Results of Terrain Correction

Terrain correction of all the data was performed by using the correction factors obtained in the previous section. In the AVE correction stage, the forest area was divided into four cover types through using the land-use data product, and then the final correction results were merged for analysis. In this section, only the data results of 14 July 2016 were presented to analyze the effects of each correction stage. Figure 6 presents the backscatter coefficients’ variation of the original data and each correction stage in different polarization channels. The horizontal axis and vertical axis are the longitude and latitude of the image, and the different colors represent the intensity of backscatter coefficients. Figure 6a1, b1, and c1 show the backscatter coefficients’ distribution of different polarization channels in the original data, and the results after POA compensation are shown in Figure 6a2, b2, and c2. According to a visual inspection, no evident differences can be seen in the corresponding polarization channels. This means that the contribution of POA compensation to terrain correction is limited. That is because the impact of the azimuth slope is relatively weak compared to the distance direction for the full polarimetric data. Figure 6a3, b3, and c3 show the results of ESA correction. Obviously, the topographic effects have been improved in all of the polarization channels. However, there are still some topographic effects in high elevation areas, such as the ridge where local incidence angles are usually relatively large, which requires further correction through the AVE correction stage. As shown in Figure 6a4, b4, and c3, in the three polarization channels, the topographic effects of ridges have effectively been removed.

To further illuminate the effects of terrain correction, we show the relationship between the local incidence angle and the backscatter coefficients of different polarization channels. The results are shown in Figure 7, where the red dashed line is the fitting curve of the backscatter coefficients and local incidence angle, and the different colors represent the density of points. We notice that there is a linear relationship between the backscatter coefficients and local incidence angle. The linear slope is relatively large in the POA compensation stage, but the linear slope becomes smaller following ESA correction and AVE correction. It indicates that POA compensation is limited in its ability to eliminate the impact of local complex terrain. In addition, the effect of the ESA correction stage is more considerable at low-local incidence angles than at high-local incidence angles, where it can effectively limit the overestimation of backscatter intensity. As shown in Figure 7a and Figure 7d, in the range of 5° to 15°, the distribution of backscatter coefficients is from −20 dB to −5 dB at the ESA correction stage, which is much lower than that of the POA correction stage (−10 dB to 5 dB). However, in the range of 70° to 80°, the distribution range of backscatter coefficients does not change in the two correction stages, staying between −10 dB and −20 dB. In contrast, after AVE correction, the distribution of backscatter coefficients is from −15 dB to 0 dB, while it remains unchanged at low-local incidence angle areas. It indicates that the AVE correction method is more effective at high-local incidence angles than at low-local incidence angles, and it can limit the underestimation of backscatter intensity at high-local incidence angle areas.

4.3. Backscatter Sensitivity to Forest GSV

In this section, we analyze the sensitivity between forest GSV of all 60 plots and the individual polarization channel backscatter in seven scenes of PALSAR-2 data. As an example, Figure 8 shows the scatterplots for different polarizations on 14 July 2016, which describe the relationship between the forest GSV and backscatter coefficents of the original and terrain-corrected data.

Compared with terrain correction, the dynamic range of these scatters is larger in the original data (Figure 8a–c) and shows low sensitivity to forest GSV. After POA compensation (Figure 8d–f), the correlation coefficient values are 0.129, 0.29, and 0.122 at HH, HV, and VV polarization, respectively, which indicate that the sensitivity has not been improved. After ESA correction (Figure 8g–i), the correlation coefficients are increased by 0.227, 0.27, and 0.24 at HH, HV, and VV polarization. The AVE correction stage also contributes to enhancing the sensitivity between forest GSV and backscatter where the correlation coefficients are 0.489, 0.643, and 0.473, respectively (Figure 8j–l). Clearly, the terrain correction can improve the sensibility between forest GSV and backscatter coefficients in this study area. However, we note that sample plots with too low GSV values (less than 37.06 m³/ha) still have negative effects on the sensitivity between forest GSV and backscatter at the HH and VV channels.

We also compare the correlations between forest GSV and backscatter coefficients of all the SAR data to select the most relevant data for GSV estimation. The results are summarized in

Table 3. Correlation coefficients between GSV and backscatter coefficent of all data.

Acquisition Time	HH	HV	VV
16 June 2016	0.418	0.495	0.370
30 June 2016	0374	0.561	0.216
14 July 2016	0.489	0.643	0.473
11 August 2016	0.435	0.564	0.377
25 August 2016	0.469	0.563	0.428
22 September 2016	0.182	0.545	0.123
6 October 2016	0.381	0.487	0.349

. We find that the correlation coefficient values are different for data acquired on different dates, and the acquired data on 14 July 2016 are much higher than other data. Moreover, the HV backscattering intensities of each scene of PALSAR-2 data show stronger correlations with the GSV than with HH or VV. Therefore, the PolSAR data obtained on 14 July 2016 is selected for further analysis, and the element HV backscatter is used as an individual measurement for single-variable regression models.

4.4. GSV Estimation and Mapping

Based on the results mentioned in Section 4.3, the PolSAR data obtained on 14 July 2016 was used to estimate forest GSV. The three elements of the diagonal of covariance matrix generated by this data are used as variables of regression models; among them, the element HV backscatter is used as an individual measurement for the single-variable regression models.

Table 4. Summary of regression model results.

Model	Regression Equation	R2
Direct linear	${\sigma }_{HV}^0=-17.525+0.013GSV$	0.529
Logarithmic	${\sigma }_{HV}^0=-26.608+2.296ln\left(GSV\right)$	0.601
Quadratic	${\sigma }_{HV}^0=-28.155+2.941ln\left(GSV\right)-0.066{\left(ln\left(\mathrm{GSV}\right)\right)}^2$	0.603
Exponential	${\sigma }_{HV}^0=-19.914+0.377sqrt\left(GSV\right)$	0.579
Water-Cloud analysis	${\sigma }_{HV}^0=-12.932-7.163\exp \left(-0.008GSV\right)$	0.612
Multi-variable	$ln\left(GSV\right)=-2.611+0.531{\sigma }_{HH}^0+0.031{\left({\sigma }_{HH}^0\right)}^2-1.693{\sigma }_{HV}^0-0.063{\left({\sigma }_{HV}^0\right)}^2$ $+0.255{\sigma }_{VV}^0+0.01{\left({\sigma }_{VV}^0\right)}^2$	0.674

provides the results of all the models’ fitting based on the training sample dataset. The fitting curves generated by the single-variable models are shown in Figure 9. The decisive coefficients’ R² values are 0.539 (Direct linear model, Figure 9a), 0.601 (Logarithmic model, Figure 9b), 0.603 (Quadratic model, Figure 9c), 0.579 (Exponential model, Figure 9d), and 0.612 (Water-Cloud analysis model, Figure 9e), respectively. The direct linear relationship (Figure 9a) between the backscatter coefficient and forest GSV is weak, but this phenomenon can be improved through the transformation of parameters, such as the natural logarithmic transformation of forest GSV (Figure 9b). In addition, the Water-Cloud analysis model is found to be the most reliable in the capacity of single-variable regression models, as it produces the highest coefficient of determination in the five models. Howeve, from Figure 9e, when the forest GSV is greater than 300 m³/ha, the change of the fitting curve tends to be gentle, which may limit its ability regarding estimation in higher GSV areas.

Compared with the Water-Cloud analysis model, the established multi-variable model may have greater potential to provide useful GSV estimation. We validate this possibility based on the test sample dataset, and the results are displayed in Figure 10. The validated plot of the Water-Cloud analysis model is characterized by a correlation coefficient R² value of 0.417, whose root mean square error (RMSE) is 91.075 m³/ha. For the multi-variable model, the correlation coefficient R² value is 0.630, whose RMSE is 70.965 m³/ha. Obviously, the accuracy of the multi-variable model inversion is higher than that of the Water-Cloud analysis model inversion. Therefore, the multi-variable model is used to estimate the forest GSV for the whole study region. The results are shown in Figure 11. Figure 11a is the schematic diagram for the spatial distribution of the forest GSV at the pixel scale, which shows that the range of the estimated forest GSV is 0 to 450 m³/ha. Figure 11b is the histogram of the GSV map. The mean and standard deviation of the GSV in the region are 135.759 m³/ha and 47.255 m³/ha. Furthermore, we also calculated the GSV of different land-cover types. The mean GSV of woodland, shrubbery, sparse woodland, and other forest were 137.701 m³/ha, 130.541 m³/ha, 125.991 m³/ha, and 113.759 m³/ha, respectively. The standard deviations were 45.906 m³/ha, 42.172 m³/ha, 56.274 m³/ha, and 62.051 m³/ha, respectively.

5. Discussion

PALSAR-2 L-band full polarimetric data has been widely used in the estimation of forest parameters. However, its use for the estimation of forest GSV presents a challenge in subtropical mountain areas, where the underlying topography is complex and diverse, seriously affecting the radiometric quality of SAR images [59]. We performed terrain correction through integrating three stages (POA, ESA, and AVE) to reduce the negative influence of topography on the full polarimetric data. In these three steps, the AVE correction step is based on a semi-empirical cosine model, which can be considered as a function of parameter n (Equation (6)). Therefore, the key to the effectiveness of the AVE correction is whether the value of n can be accurately obtained. A traditional way to obtain the value of n is to use an empirical value of one [60,61], which corrects the experimental data as a whole, and does not need to mask non-forest areas. However, it ignores the difference between polarization channels and the influence of forest features, which is often used for the terrain correction of single-polarimetric and dual-polarimetric data [39,62,63]. In this study, we calculated and evaluated the correlation results between the local incidence angle and the backscatter coefficient to generate the correction factor n of different polarization channels. It is adaptive and takes into account the difference between polarization channels. In addition, according to the results of Figure 5 and Table 2, we would like to stress the necessity of land-cover types for reducing the impact of microtopography in AVE correction, although it has many classification criteria.

In this paper, the SRTM DEM is used as an auxiliary data for SAR dataset geocoding and terrain correction. According to the results of Figure 7 and Figure 8, the effects of removing terrain and improving sensitivity are obvious. However, the DEM data are digital surface models (i.e., DSM), not digital terrain models (i.e., DTM). Thus, the experimental process is carried out under the assumption that the fluctuation of the forest canopy top is consistent with that of the underlying topography. For coarse resolution SAR data, this assumption is not a serious limitation, because it is difficult to reflect the information of individual canopy fluctuation in DSM data with coarse resolution [50]. Therefore, in the premise of multi-look processing of SAR data, SRTM DEM (30 m or 90 m) can be used to assist terrain correction [50,61]. Based on the resolution of DEM data being far lower than the original PolSAR data, we used the projection angle method instead of the area integration method [48] in the ESA correction stage. Although the use of globally shared DEM products can reduce the impact of terrain, it is still worth looking forward to obtaining high-precision and high-resolution DTM data through PolInSAR technology.

Forest GSV has different sensitivity to PolSAR data backscatter coefficients from different dates, even for the same polarization channel (Table 3). This phenomenon may be related to external environmental conditions (e.g., moisture and wind speed variations). This is because irregular variation of the environment affects the interaction between the electromagnetic waves and vegetation components. Although multi-temporal SAR data with significant climatic difference have been used to assess the relationships between the backscatter coefficients and forest parameters [28,42], the specific impact of the external environment on PolSAR data in subtropical regions is still unknown, and should be further studied. In addition, from Table 3, we find that the cross (HV)-polarized backscatter intensity of each scene of PALSAR-2 data is more sensitive to forest GSV than co (HH, VV)-polarization in subtropical mountain areas. The most likely reason is that the cross-polarized backscatter mainly occurs from multiple scattering within the tree canopy, and is less affected by the external environment [13].

We use five single-variable models to establish relationships between the GSV and the backscatter coefficients of the HV polarization (Table 4 and Figure 9). The results suggest that the direct linear relationship between the backscatter coefficient and forest GSV is weak. However, this phenomenon can be improved through the transformation of parameters, such as natural logarithmic transformation of forest GSV. Through the contrast analysis of the fitting performance of the five models, the Water-Cloud analysis model is found to be the most reliable. However, the multi-variable model has greater potential to provide a useful estimation of GSV than the Water-Cloud model. In fact, the correlation coefficient R² values of the two models only differ by 0.062, but the former has a higher accuracy of GSV estimation than the latter in the model test. This indicates that co-polarization can also make a certain contribution in GSV estimation. Therefore, our study recommends using the multi-variable model to map the GSV in the study area.

6. Conclusions

This research investigated the capability of full polarized L-band backscattering for the estimation of forest GSV in a subtropical mountain region of eastern Hunan, China. However, it was challenging to estimate GSV due to complex topography of the region. In this paper, we proposed a strategy for fine terrain correction through integrating three stages (POA, ESA, and AVE) and taking into account the impact of land-cover types. In the AVE correction stage, we calculated and evaluated the correlation results between the local incidence angle and the backscatter coefficient to generate the correction factor n of different polarization channels. We found that the optimum n values of HV polarization were within the range of zero to one, while the HH and VV polarizations were greater than one. The results of terrain correction demonstrated that the terrain correction strategy effectively reduced the negative influence of topography and improved the sensitivity between the forest GSV and backscatter coefficients. The results also showed that the land-cover types were necessary data for reducing the impact of microtopography in AVE correction. In the three primary diagonal elements of the PolSAR covariance matrix, the cross-polarized backscatter was more sensitive to forest GSV than co-polarization, and could be used as a single variable for GSV estimation. To estimate and map the GSV of the study area, five single-variable models and a multi-variable model were built using field measurements and corrected PolSAR data. The multi-variable model that was constructed by combining three diagonal elements had greater potential to provide a useful estimation of GSV than the single-variable models, whose correlation coefficient value was 0.630 and RMSE was 70.965 m³/ha. Therefore, our study recommended using the multi-variable model to map the GSV in the study area. The range of estimated forest GSV was 0 to 450 m³/ha. The mean value and stander deviation were 135.759 m³/ha and 47.255 m³/ha, respectively. The study expands the application potential of PolSAR data in complex topographic areas; thus, it is helpful and valuable for the large-scale (e.g., national or global scale) estimation of forest parameters.