Mapping Growing Stem Volume of Chinese Fir Plantation Using a Saturation-based Multivariate Method and Quad-polarimetric SAR Images

For the planning and sustainable management of forest resources, well-managed plantations are of great significance to mitigate the decrease of forested areas. Monitoring these planted forests is essential for forest resource inventories. In this study, two ALOS PALSAR-2 quad-polarimetric synthetic aperture radar (SAR) images and ground measurements were employed to estimate growing stem volume (GSV) of Chinese fir plantations located in a hilly area of southern China. To investigate the relationships between forest GSV and polarization characteristics, single and fused variables were derived by the Yamaguchi decomposition and the saturation value of GSV was estimated using a semi-exponential empirical model as a base model. Based on the estimated saturation values and relative root mean square error (RRMSE), the single and fused characteristics and corresponding models were selected and integrated, which led to a novel saturation-based multivariate method used to improve the GSV estimation and mapping of Chinese fir plantations. The new findings included: (1) All the original polarimetric characteristics, statistically, were not significantly correlated with the forest GSV, and their logarithm and ratio transformation fused variables greatly improved the correlations, thus the estimation accuracy of the forest GSV. (2) The logarithm transformation of surface scattering resulted in the greatest saturation, value but the logarithm transformation of double-bounce scattering resulted in the smallest RRMSE of the GSV estimates. (3) Compared with the single transformations, the fused variables led to more reasonable saturation values and obviously reduced the values of RRMSE. (4) The saturation-based multivariate method led to more accurate estimates of the forest GSV than the univariate method, with the smallest value (29.64%) of RRMSE achieved using the set of six variables. This implied that the novel saturation-based multivariate method provided greater potential to improve the estimation and mapping of the forest GSV.


Introduction
Through carbon sequestration and biomass accumulation, planted forests play a major role in mitigation of global climate change, due to the reduction of natural forests and increase of plantations The previous studies of using SAR images to determine the saturation values and estimate GSV mainly focused on the use of semi-empirical models and boreal and natural forests (Table 1), and the obtained saturation values of GSV varied from 80 to 595 m 3 /ha. There have been very few reports that dealt with planted forests. The area of planted forests in China ranks the highest in the world. Chinese Remote Sens. 2019, 11, 1872 4 of 20 fir is the top tree species out of the plantations in China and distributed in more than 10 provinces located in the East, Southeast, South and central south, and southwestern parts of China. Thus, Chinese fir is the most important species for wood supply in China and its contribution to sustainable forest ecosystems and global carbon cycling is also very significant. There have been few studies to explore the saturation of GSV using SAR and optical images for estimation of Chinese fir forest biomass and GSV. The planted Chinese fir forests have complicated and different canopy structures from other species because of different climatic, topographic and soil properties. Determining the saturation values of GSV and accurately deriving its spatial distribution become very important.
The previous methods for mapping GSV of Chinese fir plantations have concentrated on the use of multi-source remote sensing data, but using quad-polarimetric SAR images has been rarely reported [41][42][43]. In this study, two ALOS PALSAR-2 quad-polarimetric SAR images acquired on June 30 and August 25 in 2016 were selected to estimate the GSV values of Chinese fir plantations located in a hilly area of southern China. We analyzed and compared the single and fused polarimetric characteristics to interpret forest structure parameters and estimate the saturation levels using the aforementioned semi-exponential empirical model as a base model. We then developed a novel saturation-based multivariate method in which the saturation values were considered as a model parameter and the saturation-based models with different characteristics and saturation levels were selected and integrated to map the forest GSV of the study area.
The paper was organized as follows. Section 2 introduced the study area and data from the measured plots and the acquired quad-polarimetric SAR images. The polarimetric SAR data pre-processing and the methods for estimating the saturation levels and mapping forest GSV were demonstrated in Section 3. In Section 4, the results of the saturation levels and the forest GSV estimated by the univariate and multivariate methods were compared. The factors affecting the saturation levels and the estimation accuracy of forest GSV were discussed in Section 5 and the conclusions were drawn in Section 6.

Study Area
This study was conducted in the Huangfengqiao state-owned forest farm in Youxian county (113 • 24 N, 27 • 15 E), Hunan province of China ( Figure 1a). The study area has a hilly landscape with the elevation varying from 115 m to 1270 m and a slope ranging from 20 to 35 degrees. Located in the subtropical monsoon climate zone, this area has an annual average temperature of 17.8 • C, an average precipitation of 1410.8 mm and an average frost-free period of 292 days. About 86.24% of the region (10,122.6 ha) is covered by Chinese fir (Cunninghamia Lanceolata), Pinus massoniana Lamb, bamboo, Liriodendron chinense and Cinnamomum camphora. The planted Chinese fir is dominant species in the hilly terrain (Figure 1b).

Ground Data Collection and Processing
Most Chinese fir plantations were distributed in the northern and eastern parts of the study area. By the approach of random stratification sampling based on the spatial distribution of the Chinese fir plantations and age groups, a total of 50 plots were measured from 2016 to 2017 ( Figure  2a) and as examples, the photos in Figure 2b showed three ground measured plots with young, immature and mature forests. The random stratification sampling led to the sample plots that had a consistent spatial distribution with that of the Chinese fir plantations in the study area (Figure 1b versus Figure 2a). All the sample plots were pure Chinese fir plantation forests. There were only 6 sample plots in which Chinese fir was mixed with broad-leaved tree species and Masson pine, but the percentages of broad-leaved tree species and Masson pine were less than 8%. The ground measured plot GSV values were used to investigate the relationship between polarimetric characteristics and the Chinese fir forest GSV. The plots had a size of 30 m × 30 m or 20 m × 20 m depending on the topographic features. The locations of the selected plots, including the corners and central points, were accurately obtained using the real-time dynamic measurement (RTK) global positioning system (GPS). In each plot, the height and the diameter at breast height (DBH) of each tree were measured. Trees with the DBH equal to or greater than 5 cm were considered in the inventory.
In each plot, the stem volume of the i-th tree was estimated by: in which i V is the stem volume, i g is the cross-section area related to DBH, i H is the height of the ith tree and  f is the trunk taper coefficient of the planted Chinese fir associated with height and DBH [44][45][46]. The GSV of each plot was the sum of all the tree stem volumes within the plot. In the study, the maximum DBH was 29.48 cm and the maximum height was 20.5 m. The GSV values of the forests at different age groups varied between 0 m 3 /ha and 480 m 3 /ha and statistical parameters of ground measured plots are listed in Table 2.

Ground Data Collection and Processing
Most Chinese fir plantations were distributed in the northern and eastern parts of the study area. By the approach of random stratification sampling based on the spatial distribution of the Chinese fir plantations and age groups, a total of 50 plots were measured from 2016 to 2017 ( Figure 2a) and as examples, the photos in Figure 2b showed three ground measured plots with young, immature and mature forests. The random stratification sampling led to the sample plots that had a consistent spatial distribution with that of the Chinese fir plantations in the study area (Figure 1b versus Figure 2a). All the sample plots were pure Chinese fir plantation forests. There were only 6 sample plots in which Chinese fir was mixed with broad-leaved tree species and Masson pine, but the percentages of broad-leaved tree species and Masson pine were less than 8%. The ground measured plot GSV values were used to investigate the relationship between polarimetric characteristics and the Chinese fir forest GSV. The plots had a size of 30 m × 30 m or 20 m × 20 m depending on the topographic features. The locations of the selected plots, including the corners and central points, were accurately obtained using the real-time dynamic measurement (RTK) global positioning system (GPS). In each plot, the height and the diameter at breast height (DBH) of each tree were measured. Trees with the DBH equal to or greater than 5 cm were considered in the inventory.
In each plot, the stem volume of the i-th tree was estimated by: in which V i is the stem volume, g i is the cross-section area related to DBH, H i is the height of the ith tree and f ε is the trunk taper coefficient of the planted Chinese fir associated with height and DBH [44][45][46]. The GSV of each plot was the sum of all the tree stem volumes within the plot. In the study, the maximum DBH was 29.48 cm and the maximum height was 20.5 m. The GSV values of the forests at different age groups varied between 0 m 3 /ha and 480 m 3 /ha and statistical parameters of ground measured plots are listed in Table 2.

Quad-Polarimetric SAR Data and Digital Elevation Model (DEM)
The two L-band quad-polarimetric SAR images over the study area were downloaded from the Japanese Aerospace Exploration Agency (http://global.jaxa.jp/), which were acquired on 30 June 2016 and 25 August 2016, on descending orbit with an incidence angle of about 38.99 degrees. The resolutions in the azimuth direction and the slant range direction were 2.83 m and 2.86 m, respectively. Based on the information from the local weather forecast, the weather when the SAR images were acquired on 30 June and 25 August was cloudy and showering, respectively. Moreover, it was rainy for three days before the image was acquired on 25 August. Figure 3a shows a Pauli

Quad-Polarimetric SAR Data and Digital Elevation Model (DEM)
The two L-band quad-polarimetric SAR images over the study area were downloaded from the Japanese Aerospace Exploration Agency (http://global.jaxa.jp/), which were acquired on 30 June 2016 and 25 August 2016, on descending orbit with an incidence angle of about 38.99 degrees. The resolutions in the azimuth direction and the slant range direction were 2.83 m and 2.86 m, respectively. Based on the information from the local weather forecast, the weather when the SAR images were acquired on 30 June and 25 August was cloudy and showering, respectively. Moreover, it was rainy for three days before the image was acquired on 25 August.

SAR Data Pre-processing
Two single-look complex (SLC) PALSAR data sets were selected to investigate the saturation level and map the GSV of Chinese fir plantations. The polarimetrical calibration was performed to reduce the impact of Faraday rotation. Then, the Papathanassiou algorithm was employed to calibrate the L-band ALOS PALSAR-2 quad-polarimetric SAR images [47][48][49][50][51][52]. The speckle noise was reduced using the Lee filter (7  7). After that, the coherency matrix was generated from the calibrated and filtered scatter matrix. The terrain slope is also an important factor for SAR image calibration. The polarization orientation angle (POA) between the assumed and the local horizontal polarization vectors [53][54][55] should be compensated. The POA was estimated by the circular polarization as follows [56][57][58][59][60][61]: where, Arg is the phase of complex data， RR S and LL S are the circular polarization components [7,48]. The facet method related to DEM was applied to deal with the terrain radiometric correction (TRC). The relationship between the backscatter coefficient of object 0  and radar brightness 0  is [50][51][52]:

SAR Data Pre-processing
Two single-look complex (SLC) PALSAR data sets were selected to investigate the saturation level and map the GSV of Chinese fir plantations. The polarimetrical calibration was performed to reduce the impact of Faraday rotation. Then, the Papathanassiou algorithm was employed to calibrate the L-band ALOS PALSAR-2 quad-polarimetric SAR images [47][48][49][50][51][52]. The speckle noise was reduced using the Lee filter (7 × 7). After that, the coherency matrix was generated from the calibrated and filtered scatter matrix. The terrain slope is also an important factor for SAR image calibration. The polarization orientation angle (POA) between the assumed and the local horizontal polarization vectors [53][54][55] should be compensated. The POA was estimated by the circular polarization as follows [56][57][58][59][60][61]: where, Arg is the phase of complex data, S RR and S LL are the circular polarization components [7,48]. The facet method related to DEM was applied to deal with the terrain radiometric correction (TRC). The relationship between the backscatter coefficient of object σ 0 and radar brightness β 0 is [50][51][52]: where σ 0 is the cross area, δ r is the range resolution, δ a is the azimuth resolution and dσ is related to the local terrain and geometric parameters of SAR images. The software Polsarpro was used to process the SAR data and obtain the polarimetric characteristics, and the software SARscape was utilized to conduct the geocoding.

Retrival of Polarimetric Characteristics
The four-component decomposition approach proposed by Yamaguchi can lead to polarimetric characteristics from the coherency matrix without using any ground measurements, and was used to extract the scattering features, including the surface scattering (Odd), double-bounce scattering (Dbl), volume scattering (Vol) and helix scattering (Hlx) [24,25,28,29]. The contributions of the four scatterings are related to the wavelength, incidence angle, forest structure properties, canopy shape and terrain. The power of these scattering features was used to estimate the span (P t ) as follows: where P Odd , P Dbl , P Vol and P Hlx are the powers of surface scattering, double-bounce scattering, volume scattering and helix scattering, respectively. A series of polarimetric characteristics after logarithm and ratio transformation were also used to estimate the GSV. The logarithm transformation used to extract the polarimetric characteristics is: where P i−dB is the polarimetric characteristic total after the logarithm transformation, and P i is the original polarimetric characteristic total. The ratio transformation was used to get the relative changes between different scattering mechanisms as follows in which, P i−RA is the polarimetric characteristics after the ratio transformation and P t is the span of each pixel. In addition, polarimetric characteristic combinations formed by multiplication or division of the single polarimetric characteristics were called the fused variables such as P Dbl/Odd , P Vol/Odd , P Dbl×Vol and P Dbl×Vol/ Odd , and employed to estimate GSV. The single and fused polarimetric characteristics were compared in terms of the sensitivity to forest GSV.

Forest GSV Estimation
In this study, an empirical model with an exponential form proposed by Wagner et al. [39,40] was employed as a base model to estimate the saturation level, since the model obeys the scattering mechanisms, and can simply account for the relationships between the polarimetric characteristics and forest GSV [12][13][14][15][16][17][18][19][20]39,40]: where σ 0 GSV is one of the selected polarimetric characteristics and GSV is the measured GSV (m 3 /ha). β n refers to the polarimetric characteristic of non-vegetated area, β s refers to the characteristics of the forests with the highest GSV, and k is the saturation level of the forest's GSV. β n , β s and k are unknown parameters, whose initial values were determined by the range of polarimetric characteristics. The non-linear algorithm was employed to estimate the unknown parameters. The forest GSV (m 3 /ha) was retrieved using the obtained univariate exponential model [12,40]: The forest GSV reaches the saturation level when σ 0 GSV is larger than β s . Both the range of GSV and the accuracy of GSV estimates are dependent on the saturation level. If the values of GSV exceed the saturation level, the GSV estimation errors of the univariate method would obviously increase. The model was used as a base model to explore the sensitivity of the single and fused characteristics from the images to the estimation of the saturation levels and forest GSV.
The polarimetric characteristics with different saturation levels are sensitive to the forest GSV differently. Selecting and integrating these polarimetric characteristics may provide the great potential to improve the forest GSV estimation accuracy. In this study, a saturation-based multivariate approach for estimating the forest GSV was proposed by selecting the single and fused characteristics that accurately captured the forest canopy structures and tree trunk features, and then combining the corresponding models that could lead to reasonable saturation levels. The selection was conducted based on the smallest relative root mean square error (RRMSE). In each pixel, the saturation levels were considered as weight coefficients. The saturation-based multivariate model could be expressed as follows: where GSV f (i, j) is the GSV value of pixel (i, j) obtained by the saturation-based multivariate method. GSV(i, j, n) is the GSV value of pixel (i, j) estimated by the univariate model using the nth independent variable, and k n is the saturation level of the nth independent variable estimated by Equation (7); that is, the nth univariate model. The N is the number of the selected univariate models. The r denotes an adjustment coefficient that is used to help indirectly quantify the reasonability of a GSV estimate. In this study, r was determined to range from 0.8 to 2, implying that a GSV estimate would be reasonable if it fell within the range of a 0.8 timing saturation value to a 2 timing saturation value. The adjustment coefficient (r) changes from 0.8 to 2 by an interval of 0.05, which could lead to a total of 24 r values.
With each of the r values, a product of r with k n could be created and utilized to help measure the reasonability of GSV estimates. Given an r value, the reasonable estimates were then compared with the observed GSV values at the plot level and an RRMSE value was yielded. The optimal value of r could be finally determined using the smallest RRMSE value.
In order to obtain a forest GSV map using Equation (9), first, the unknown parameters of each semi-exponential empirical model Equation (7) with the selected polarimetric characteristics were solved by a non-linear algorithm. After that, the GSV in Equation (8) was estimated by the selected single polarimetric characteristics on the basis of pixel by pixel. Given a set of multiple polarimetric characteristics, the same set of forest GSV maps were generated. This meant that at each pixel, the same set of GSV estimates were obtained. The plot estimates were used to determine the optimal value of the adjustment coefficient r based on the smallest RRMSE. The optimal value of r was utilized to select and integrate the reasonable estimates in Equation (9), which would lead to the final forest product of GSV. Estimating GSV was conducted using the programs made by the first author based on Matlab, and mapping forest GSV was finished using ArcGIS.

Polarimetric Characteristics
The polarimetric characteristics derived by the Yamaguchi decomposition were employed to interpret the interaction between forest GSV and the independent variables at the plot level ( Figure 4). The Pearson correlation coefficient (γ) was adopted to select the characteristics for the GSV estimation based on a statistical test of correlation: Whether or not the correlation coefficients statistically were significantly different from zero at the significance level of 0.01. To match the size of the measured plots, the power value of scatterings in each plot was obtained by a window spatial averaging. Since the power of helix scattering (−30 dB to 10 dB) was too low (Figure 4d) compared with other scatterings (Figure 4a-c), it was not considered in the subsequent analysis. Moreover, the power of volume scattering (less than −8 dB) (Figure 4c) was also lower than those of the surface scattering (−5 dB to −2 dB) ( Figure 4a) and double-bounce scattering (−7 dB to −3 dB) (Figure 4b). Before the transformation, the power of double-bounce ( Figure 4b) and volume scattering (Figure 4c) increased slightly with the increasing GSV. After the logarithm and ratio transformations, the positive correlation became stronger (Figure 4f,g,j,k). All the correlation coefficients of the original scattering components were significantly smaller than the critical value of 0.345 at the risk level of 0.01 and these components were thus omitted in the subsequent analysis.
interpret the interaction between forest GSV and the independent variables at the plot level ( Figure   4). The Pearson correlation coefficient (  ) was adopted to select the characteristics for the GSV estimation based on a statistical test of correlation: Whether or not the correlation coefficients statistically were significantly different from zero at the significance level of 0.01. To match the size of the measured plots, the power value of scatterings in each plot was obtained by a window spatial averaging. Since the power of helix scattering (-30 dB to 10 dB) was too low (Figure 4d) compared with other scatterings (Figure 4a,b,c), it was not considered in the subsequent analysis. Moreover, the power of volume scattering (less than -8 dB) (Figure c) was also lower than those of the surface scattering (-5 dB to -2 dB) ( Figure 4a) and double-bounce scattering (-7 dB to -3 dB) (Figure 4b). Before the transformation, the power of double-bounce ( Figure 4b) and volume scattering ( Figure  4c) increased slightly with the increasing GSV. After the logarithm and ratio transformations, the positive correlation became stronger (Figure 4f,g,j,k). All the correlation coefficients of the original scattering components were significantly smaller than the critical value of 0.345 at the risk level of 0.01 and these components were thus omitted in the subsequent analysis.  Table 3, after the logarithm and ratio transformation, the Pearson's correlation coefficients were greatly improved and become statistically different from zero, with the maximum of 0.70 achieved by the logarithm transformation of double-bounce scattering. Compared with the original polarimetric characteristics, the combined or fused characteristics also significantly  Table 3, after the logarithm and ratio transformation, the Pearson's correlation coefficients were greatly improved and become statistically different from zero, with the maximum of 0.70 achieved by the logarithm transformation of double-bounce scattering. Compared with the original polarimetric characteristics, the combined or fused characteristics also significantly increased the correlation coefficients with the forest plot GSV. After the transformations, the single and fused polarimetric characteristics were significantly correlated with the forest GSV and considered as potential independent variables to estimate the forest GSV.

Saturation Level of Planted Forest
Before the estimation, the plots in the shadow regions were discarded. The values of three unknown parameters were obtained by solving the base empirical model using the non-linear algorithm and the proposed initial values (Table 4). For the image acquired on 30 June 2016, the estimated saturation values ranged from 140.05 m 3 /ha to 349.84 m 3 /ha, and the greatest value was obtained by the power of surface scattering. Moreover, the saturation values varied with the polarimetric characteristics. The saturation value obtained using the power of the double-bounce scattering after the ratio transformation, 260.88 m 3 /ha, was slightly larger than that obtained using the double-bounce scattering after the logarithm transformation, 188.13 m 3 /ha. Additionally, the saturation values derived from the image acquired on 25 August 2016 were much smaller than those derived from the images acquired on 30 June 2016. Some of the estimated saturation values were smaller than 100 m 3 /ha. The saturation values derived from the fused characteristics ranged from 140.05 m 3 /ha to 206.42 m 3 /ha for the image acquired on 30 June and from 89.80 m 3 /ha to 193.71 m 3 /ha for the image acquired on 25 August, which were slightly smaller than those by single variables ( Table 4). As shown in Figure 5, the power of the surface scattering (Figure 5a,c) and the double-bounce scattering (Figure 5b,d) had the highest saturation level for the image dated on 30 June. On the contrary, it was obvious that the ranges of the saturation levels obtained by the fused variables were more stable (Figure 5e-h).
Remote Sens. 2019, 11 12 of 4 that the ranges of the saturation levels obtained by the fused variables were more stable (Figure 5eh).

Forest GSV Estimated by the Univariate Method
With the estimated parameters, the GSV values of the planted forests were retrieved by the univariate models ( Table 5). The leave-one-out cross-validation (LOOCV) was utilized for the accuracy assessment [60] based on the root mean square error (RMSE) and the coefficient of determination ( 2 R ) between the estimated and observed GSV values. The relative RMSE (RRMSE = RMSE × 100 / sample mean) was also employed. We used 120 m 3 /ha as a deviation threshold value between the measured and estimated GSV values and counted the number of the plots with the errors exceeding the threshold. The reason for using the 120 m 3 /ha as the deviation threshold value was because this value was 50% of the sample mean 240 m 3 /ha and when the error of an estimate was larger than 50% of the sample mean, the estimate could be considered to be highly problematic.
In Table 5, scattering mechanisms had great impacts on the forest GSV estimation accuracy. As mentioned previously, the power of volume scattering was insufficient to describe the relationship between the forest GSV and the polarimetric characteristics. The power of surface scattering could successfully estimate the saturation value but more than one third of the selected plots were excluded, since their errors exceeded the given threshold. The minimum RMSE (69.19 m 3 /ha) and RRMSE (33.94%) between the estimated and measured GSV were obtained by the logarithm transformation of Dbl. The Dbl was thus more sensitive to the change of the forest GSV than the others.
The fused polarimetric characteristics were also used to estimate the forest GSV (Table 5). For the image acquired on 30 June, the fused characteristics led to the determination coefficients (R 2 ) ranging from 0.53 to 0.66 and the maximum coefficient (0.66) was obtained by Dbl × Vol. The RMSE ranged from 71.04 m 3 /ha to 75.62 m 3 /ha. Therefore, using the fused characteristics could improve the accuracy of forest GSV estimation. The correlation coefficients and saturation levels indicated that it

Forest GSV Estimated by the Univariate Method
With the estimated parameters, the GSV values of the planted forests were retrieved by the univariate models ( Table 5). The leave-one-out cross-validation (LOOCV) was utilized for the accuracy assessment [60] based on the root mean square error (RMSE) and the coefficient of determination (R 2 ) between the estimated and observed GSV values. The relative RMSE (RRMSE = RMSE × 100/sample mean) was also employed. We used 120 m 3 /ha as a deviation threshold value between the measured and estimated GSV values and counted the number of the plots with the errors exceeding the threshold. The reason for using the 120 m 3 /ha as the deviation threshold value was because this value was 50% of the sample mean 240 m 3 /ha and when the error of an estimate was larger than 50% of the sample mean, the estimate could be considered to be highly problematic. Table 5. The GSV estimation accuracy of the selected variables (Note: * and ** indicate the percentages of the plots with errors exceeding the threshold are larger than 10% and 20%, respectively).

Variable
Image of Table 5, scattering mechanisms had great impacts on the forest GSV estimation accuracy. As mentioned previously, the power of volume scattering was insufficient to describe the relationship Remote Sens. 2019, 11, 1872 13 of 20 between the forest GSV and the polarimetric characteristics. The power of surface scattering could successfully estimate the saturation value but more than one third of the selected plots were excluded, since their errors exceeded the given threshold. The minimum RMSE (69.19 m 3 /ha) and RRMSE (33.94%) between the estimated and measured GSV were obtained by the logarithm transformation of Dbl. The Dbl was thus more sensitive to the change of the forest GSV than the others.
The fused polarimetric characteristics were also used to estimate the forest GSV (Table 5). For the image acquired on 30 June, the fused characteristics led to the determination coefficients (R 2 ) ranging from 0.53 to 0.66 and the maximum coefficient (0.66) was obtained by Dbl × Vol. The RMSE ranged from 71.04 m 3 /ha to 75.62 m 3 /ha. Therefore, using the fused characteristics could improve the accuracy of forest GSV estimation. The correlation coefficients and saturation levels indicated that it was very hard to accurately estimate forest GSV using the single variables from the image acquired on 25 August. Figure 6 compared the measured and estimated GSV values by the univariate method. The estimated GSV values were highly affected by the saturation value. The power of Dbl with high saturation levels (Figure 6a,b) had fewer plots with errors exceeding the threshold than that with low saturation levels (Figure 6c was very hard to accurately estimate forest GSV using the single variables from the image acquired on 25 August. Table 5. The GSV estimation accuracy of the selected variables (Note: * and ** indicate the percentages of the plots with errors exceeding the threshold are larger than 10% and 20%, respectively).  Figure 6 compared the measured and estimated GSV values by the univariate method. The estimated GSV values were highly affected by the saturation value. The power of Dbl with high saturation levels (Figure 6a,b) had fewer plots with errors exceeding the threshold than that with low saturation levels ( Fig. 6c-f). The selected polarimetric characteristics became insensitive to the change of forest GSV if the GSV values exceeded the saturation value (Figure 6e,f). The independent variables with different saturation levels had different sensitivities to the forest GSV. The range of the forest GSV estimates depended largely on both the independent variables and the obtained saturation levels. The selected polarimetric characteristics with low saturation levels could not be used to generate high values of forest GSV. The base empirical model is a univariate approach, and in this approach, just one polarimetric characteristic was selected and others were discarded because of low saturation levels. In fact, the discarded polarimetric characteristics might also contain useful information and should be introduced into the multivariate method to estimate forest GSV.

Forest GSV Estimated by the Saturation-Based Multivariate Method
To overcome the disadvantage of the base empirical model, we developed the saturation-based multivariate method, Equation (9), to improve the estimation of the forest GSV. According to the property of scattering mechanisms and the estimated saturation levels, three multivariate sets consisting of eight, six and four variables, respectively, were selected from each image to investigate the sensitivities of polarimetric characteristics. The set of the four variables consisted of four fused characteristics, including Dbl/Odd, Vol/Odd, Dbl ×Vol and Dbl × Vol/Odd. The set of the six variables consisted of the four fused variables, and the logarithm and ratio transformations of Dbl. The eight variables included the four fused variables, and the logarithm and ratio transformations of Odd and Dbl. The results of the saturation-based multivariate method are listed in Table 6. Compared with the univariate method, the saturation-based multivariate method led to higher estimation accuracy. The RRMSE ranged from 29.64% to 35.09% for the image acquired on June 30 and from 36.23% to 44.70% for the image acquired on 25 August. Moreover, the smallest RMSE and RRMSE values were 65.79 m 3 /ha and 29.64%, obtained by the set of the six variables that were involved in Equation (9). Additionally, the scatter graphs between the observed and estimated GSV values (Figure 6g,i) showed that compared with those from the univariate method, overall, the overestimations and underestimation were greatly reduced by the saturation-based multivariate method. Specially, for the plot that had the forest GSV values larger than the saturation points, the estimation accuracy was significantly improved by the saturation-based multivariate method. In Figure 7

Polarimetric Characteristics Related to Tree Species
In this study, the measured GSV of Chinese fir plantations showed strong positive correlations with the power of double-bounce scattering, and strong negative correlations with the power of surface scattering (Table 3), but weak correlations with the power of volume scattering. Several factors may influence the correlations between the scattering mechanisms and forest GSV, and tree species might be the dominant factor. In the mixed forests of Siberia, the power of the volume scattering derived by the Yamaguchi decomposition showed high correlations with the measured GSV using L-band ALOS PALSAR data [20,21,29,31,36]. In the fast-growing forests in Indonesia, the surface scattering showed a negative correlation with the forest GSV (−0.70) and the volume scattering had a positive correlation with the forest GSV (0.45), but the double-bounce scattering was

Polarimetric Characteristics Related to Tree Species
In this study, the measured GSV of Chinese fir plantations showed strong positive correlations with the power of double-bounce scattering, and strong negative correlations with the power of surface scattering (Table 3), but weak correlations with the power of volume scattering. Several factors may influence the correlations between the scattering mechanisms and forest GSV, and tree species might be the dominant factor. In the mixed forests of Siberia, the power of the volume scattering derived by the Yamaguchi decomposition showed high correlations with the measured GSV using L-band ALOS PALSAR data [20,21,29,31,36]. In the fast-growing forests in Indonesia, the surface scattering showed a negative correlation with the forest GSV (−0.70) and the volume scattering had a positive correlation with the forest GSV (0.45), but the double-bounce scattering was weakly correlated with the forest GSV [29]. In tropical forests, the power of volume scattering derived from the airborne L-band SAR data had positive correlations with the forest GSV [36]. The reasons for the different results might be mainly due to different tree species, thus different forest canopy structures. In Siberia, the percentage of broad-leaf species ranged from 29% to 44% [20,21,31], and in Indonesia and tropical forests [29,36], broad-leaf species were considered as dominated tree species and the power of volume scattering became strong by increasing the density of leaves. The strong correlations between the power of volume scattering and forest GSV were thus observed in the studies. In this study, however, the planted Chinese fir forests had narrow and short needles that led to many canopy gaps, thus the power of volume scattering could not account for the interaction between the forest canopy structures and the polarimetric characteristics. Moreover, tree trunks and big branches were the dominant scattering objects under the canopies, thus the double-bounce scattering was much stronger than that of volume scattering. This was especially true in the mature and over mature stands.

Saturation Level of Forest GSV
For a single SAR image, the saturation value is determined by two factors, including the estimation model and the polarimetric characteristics used. The empirical models have been widely employed to map GSV, in which the saturation level is estimated as a parameter. In this study, using the base empirical model with different independent variables, led to the different saturation values of GSV for the Chinese fir forests. The saturation values ranged from 140.05 m 3 /ha to 349.84 m 3 /ha (Table 4) and the greatest value was obtained using the power of surface scattering. The saturation values estimated by the double bounce scattering acquired on 25 August were about 63 m 3 /ha, seeming unreasonable. Some of variables, such as volume scattering, failed to result in the saturation values of the forest GSV, due to their low correlation to the plot GSV.
Theoretically, the power of surface scattering seems to be more sensitive to the forest GSV but the results of this study showed different ( Table 5). The major reason was the existence of many forest canopy gaps in the Chinese fir plantations, which made the power of surface scattering related to the ground under the forest canopies, and led to a weak capacity to capture the forest GSV. This was especially true in the mature forests. The double bounce scattering was directly related to tree trunks and big branches and was very sensitive to the changes of the forest GSV. Furthermore, the fused variables formed by multiplication or division of the polarimetric characteristics greatly reduced the overestimations and underestimations of the saturation values, thus greatly improved the estimation accuracy. The saturation levels from the models that used the fused polarimetric characteristics were more reasonable and stable.
The great differences of the saturation levels derived from two SAR images were observed in Table 4. This might be mainly due to different weather conditions when the images were acquired. When the first image dated 30 June was collected, the weather was cloudy and the clouds had less impact on the quality of the image. However, there were showers on 25 August, on which date the second image was acquired. Moreover, it had been rainy for three days before 25 August. The showers and moisture might have great impact on the quality of the second image, and thus resulted in a low Pearson's correlation between the GSV and the polarimetric characteristics from the image acquired from 25 August. The saturation levels derived from the SAR image acquired on 25 August were obviously lower than those from the first image acquired on 30 June.

Estimated GSV of Chinese Fir Plantation
The base empirical model has been frequently applied to express the interaction between polarimetric characteristics and forest GSV [36][37][38][39][40]. When one independent variable was employed, the estimated GSV values were closely related to the saturation level obtained (Table 4). For the power of Dbl after the ratio transformation, the RRMSE of the estimated GSV values with a large saturation value (260.88 m 3 /ha) was much smaller than that with a small saturation value (82.84 m 3 /ha). The error between the measured and estimated GSV values was significantly large when the measured GSV was larger than the saturation value (Figure 6a,e,f). The accuracy of the estimated GSV values also depended on the sensitivity of the independent variables. The different polarimetric characteristics led to different saturation values of the forest GSV, and thus different estimation accuracies. Overall, the fused characteristics resulted in more reasonable saturation values, thus much higher estimation accuracies of the forest GSV (Table 5).
In this study, the univariate method led to the smallest RRMSEs of 31.13% and 36.89%, for the images acquired on 30 June and 25 August 2016, respectively. Using the saturation-based multivariate method, the corresponding smallest RRMSE values were 29.64% and 36.23%, for the images acquired on 30 June and 25 August 2016, respectively. Compared with the univariate method, the saturation-based multivariate method obviously decreased the errors of the estimated GSV values. Additionally, the accuracy of the estimated GSV values was not related to the number of the used variables. The most accurate results were obtained by the saturation-based multivariate method with six variables from the image acquired on 30 June.

Conclusions
In this study, the powers of scatterings derived by the Yamaguchi decomposition were used to estimate the saturation values and forest GSV of Chinese fir plantations located in the hilly area of southern China using two quad-polarimetric PALSAR-2 images. The saturation value of the forest GSV was estimated by the univariate empirical model based on the single and fused characteristics. A novel and saturation-based multivariate method was then proposed by selecting and integrating the results from the univariate empirical models to improve the estimation accuracy of mapping the forest GSV in the study area. It was newly found that all the original polarimetric characteristics had weak correlations with the forest GSV and the power of double-bounce scattering was more sensitive to the forest GSV than the other polarimetric characteristics. Overall, the logarithm and ratio transformations of the scatterings greatly improved the correlations with the forest GSV, and further improvement of the correlation was achieved by the fused polarimetric characteristics. Moreover, the greatest saturation value was obtained using the logarithm transformation of surface scattering, but among the transformations of the polarimetric characteristics, the logarithm transformation of double-bounce scattering resulted in the smallest RRMSEs of the GSV estimates. Compared with the single transformations, the fused variables led to more reasonable saturation values and obviously reduced the values of RRMSE. More importantly, compared with the univariate method, the saturation-based multivariate method led to more accurate estimates of the forest GSV by reducing the overestimations and underestimations, with the smallest value (29.64%) of RRMSE achieved using the set of six variables. In addition, the SAR image dated 30 June provided better performance than the image dated 25 August. In summary, the novel saturation-based multivariate method provided the potential to improve the estimation accuracy of the Chinese fir forest GSV.