*Article* **Multi-Sensor Prediction of Stand Volume by a Hybrid Model of Support Vector Machine for Regression Kriging**

#### **Lin Chen 1,2, Chunying Ren 1,\*, Bai Zhang <sup>1</sup> and Zongming Wang 1,3**


Received: 23 January 2020; Accepted: 5 March 2020; Published: 6 March 2020

**Abstract:** Quantifying stand volume through open-access satellite remote sensing data supports proper management of forest stand. Because of limitations on single sensor and support vector machine for regression (SVR) as well as benefits from hybrid models, this study innovatively builds a hybrid model as support vector machine for regression kriging (SVRK) to map stand volume of the Changbai Mountains mixed forests covering 171,450 ha area based on a small training dataset (n = 928). This SVRK model integrated SVR and its residuals interpolated by ordinary kriging. To determine the importance of multi-sensor predictors from ALOS and Sentinel series, the increase in root mean square error (RMSE) of SVR was calculated by removing the variable after the standardization. The SVRK model achieved accuracy with mean error, RMSE and correlation coefficient in –2.67%, 25.30% and 0.76, respectively, based on an independent dataset (n = 464). The SVRK improved the accuracy of 9% than SVR based on RMSE values. Topographic indices from L band InSAR, backscatters of L band SAR, and texture features of VV channel from C band SAR, as well as vegetation indices of the optical sensor were contributive to explain spatial variations of stand volume. This study concluded that SVRK was a promising approach for mapping stand volume in the heterogeneous temperate forests with limited samples.

**Keywords:** ALOS-2 L band SAR; Sentinel-1 C band SAR; Sentinel-2 MSI; ALOS DSM; stand volume; support vector machine for regression; ordinary kriging

#### **1. Introduction**

Forest stand volume, as an ecosystem service, forms the basis for decision-making at diverse levels [1]. Spatial explicit information on forest stand volume is critical for indirect estimation of aboveground biomass for quantifying carbon sequestration and carbon dioxide exchange [2]. Field-based inventories of forest stand volume, the conventional approach, is costly and spatially limited [3]. Progress has been made in mapping forest volume by remote sensing modeling based on multisource satellite and inventory data for spatially continuous and temporally uniform predictions [4,5]. Those remote sensing algorithms were divided into two categories, i.e., physical and empirical models, and the latter included statistical regressions, machine learning techniques, and hybrid approaches [6–8]. Physically based models depend on numerous geometry and biochemistry factors, which may not be readily available [9,10]. Statistical regressions model stand volume by estimating equation parameters related to remote sensing variables [11,12]. These regressions have advantages on modeling explicit relationships and applications at large scales [13,14]. Machine learning algorithms have no assumption

on input variable distribution, type and number, which achieve robust and accurate predictions on complex relationships [15,16]. Among the various machine learning techniques, support vector machine is acclaimed for its capacity of dealing with small training datasets in remote sensing-based classification [17,18]. After the re-design to predict quantitative outputs and solve regression problems, this algorithm came to be the support vector machine for regression (SVR) and acquired wide successes in stand volume modeling [19,20]. Hybrid approaches involve either the statistical regression or machine learning model between the target variable and remote sensing predictors, interpolating residuals of predictions by kriging, and combining them [21–23]. Those two-step approaches both consider the spatial heterogeneity conveyed by remote sensing predictors and autocorrelation of neighboring observed data [24,25]. Those approaches, especially machine learning combined ordinary kriging of residuals such as artificial neural network kriging (ANNK) and random forest kriging (RFK), have yielded accurately spatial predictions [26,27]. However, support vector machine for regression kriging (SVRK) modeling for mapping forest volume has rarely been tested and reported.

Stand volume modeling with open-access satellite data has been comparable, repeatable, and has long-term monitoring [28–30]. With the global coverage, Sentinel-1 C band synthetic aperture radar (SAR) and Sentinel-2 multispectral instrument (MSI) images provide capabilities for stand volume modeling using both active and passive remote sensing techniques [31,32]. The Advanced Land Observing Satellite (ALOS/ALOS-2) Phased Array type L band SAR (PALSAR/PALSAR-2) from L band SAR have penetrability, which contain comprehensive information on the orientation and structure of tree canopy and stems within the pixel [33,34]. It makes the ALOS/ALOS-2 images with global observations particularly useful for stand volume mapping [35,36]. The ALOS digital surface model (DSM) from L band interferometric SAR (InSAR) with accurate values of elevation and can provide useful topographic indices to estimate stand volume [37–39]. Reported studies have explored the potential of multi-sensor data using the SVR in volume mapping [31,40,41]. However, how volume predictions would be affected by using SAR and MSI predictors based on the SVR deserves further exploration.

The Changbai Mountains Mixed forests, as the richest eco-region in temperate forests of northeastern China, play a key role in carbon cycles and ecosystem services both at regional and global scales [42–45]. Hence, in this study, we innovatively developed a SVRK model based on limited samples and open-access satellite predictors, and adopted it to map stand volume of the Changbai Mountains Mixed forests, a vital eco-region of temperate ecosystems. The specific objectives were to: (1) determine and compare the relationships of forest volume with multi-sensor variables from ALOS-2, Sentinel-1, Sentinel-2 and ALOS DSM; (2) map stand volume by the SVRK modeling; and (3) analysis spatial variations of stand volume and provide managerial suggestions for forest farms in the study area.

#### **2. Materials and Methods**

#### *2.1. Study Area and Field-Measured Stand Volume*

The study area covers 171,450 ha and 12 forest farms belonging to Forestry Bureau of Dunhua County (Figure 1). The site is located within the western mountainous area of Yanbian Korean Autonomous Prefecture of Jilin Province, northeast China. The climate is four-season, monsoon-influenced and humid continental, with an annual average temperature and precipitation of 3.28 ◦C and 632 mm, respectively [25]. Characterized by the dense cover of the Changbai Mountains mixed forests, the major forest types include deciduous broadleaved forest and mixed broadleaf-conifer forest with natural vegetation [43]. Dominant tree species include *Tilia amurensis* (Rupr.), *Juglans mandshurica* (Maxim.), *Fraxinus mandschurica* (Rupr.), *Mongolian oak* (Quercus spp.) and *Betula platyphylla* (Suk.). Typical soils are dark-brown earths, meadow, bog, chernozem, and peat soil.

**Figure 1.** The outline of the study area, sampling sites, and employed open-access satellite remote sensing data derived from Advanced Land Observing Satellite (ALOS), ALOS-2, Sentinel-1, and Sentinel-2 series. Date include (**a**) Sentinel-2 Level-1C, (**b**) ALOS-2 yearly mosaic, (**c**) Sentinel-1 Level-1 GRD, (**d**) Sentinel-2 Level-2A, and (**e**) ALOS Digital Surface Model products.

The field campaign was carried out in September 2017. Stratified sampling design was used by masking non-forest areas and randomly generating the distribution of sampling plots in forest areas, while the plots that were impossible to access were replaced by the nearest sites. Following the national guidelines for forest resource survey [46], eight teams took part in collecting measured data under the same protocol. A total of 1392 squared 30 m by 30 m samples were established (Figure 1a). At each sample site, tree species, diameter at breast height (DBH, the diameter at 1.3 m from the ground), and tree height were measured and recorded. Age classes from young to over-mature were acquired from the forest manager's archives at the local forestry bureau for further analysis (Figure 1a). Stand volume was estimated by DBH and tree height according to the National Standard of China: Tree volume tables (LY/T 1353–1999) [47]. The field-measured stand volume was from 1 to 499.8 m3/ha, and was divided into six levels with the same frequency, with the median and standard deviation (SD) value of 146.3 and 56.2 m3/ha, respectively (Figure 2a). The values of measured volume were mainly below 200 m3/ha with 85.57 % (Figure 2b). The 1,392 sampling sites were randomly divided into training (n = 928) and validation (n = 464) sets (Figure 1a) for training and assessing the models.

**Figure 2.** The values of measured stand volume. (**a**) Field sample profiles of volume in the study site from Plot 1 to 1,392; (**b**) Components of volume.

#### *2.2. Satellite Data Pre-Processing and Derived Variables*

The adopted multi-sensor satellite data are listed in Table 1. The 25-m ALOS-2 L band SAR yearly mosaic images of 2017 were downloaded from the ALOS Research and Application Project of EORC, the Japan Aerospace Exploration Agency to acquire the normalized backscatter coefficients (gamma naught values) (Table 2), which was sensitive to stand volume [36,48,49]. Images were converted to gamma naught values in decibel unit (dB) from 16-bit digital number (DN) (Figure 1b) using the following Equation (1) [50]:

$$
\gamma^0 = 10 \log\_{10} \left( DN^2 \right) - 83 \tag{1}
$$

where γ<sup>0</sup> is gamma naught backscatter coefficient of horizontal transmit-horizontal channel (HH) or horizontal transmit-vertical channel (HV); *DN* is the polarization data in HH or HV.


**Table 1.** The adopted ALOS-2, Sentinel-1, Sentinel-2, and digital surface model (DSM) data.


**Table 2.** Remote sensing indices from the ALOS and Sentinel series data for volume mapping.

Sentinel series images were downloaded from the Copernicus Sentinel Scientific Data Hub. The data included one Sentinel-1 C-band SAR and three Sentinel-2 MSI images. The SAR image was at a high-resolution (HR) Level-1 ground range detected (GRD) processing level with a pixel size of 10 m [51]. Promising results demonstrated that the normalized backscatter coefficients and texture features from Sentinel-1 images could improve forest parameter estimation [32,52]. Sentinel-1 Toolbox in SNAP software (version 6.0, European Space Agency, Paris, France) was used to acquire Sentinel-1 variables with a map projection (Figure 1c) by image calibration, speckle reduction using the Refined Lee Filter, terrain correction by the Range-Doppler, and grey level co-occurrence matrix analysis with 3 × 3-pixel window [52–54]. The Sentinel-2 Level 1C data were top-of-atmosphere reflectance, which were processed by orthorectification and registration [55]. The MSI data had 13 spectral bands, including four in 10 m (bands 2–4, 8), six in 20 m (band 5–7, 8a, 11–12), and three in 60 m (band 1, 9–10) spatial resolutions, respectively [55]. The 10-m Sentinel-2 Level 2A data (Figure 1d) were atmospherically corrected from the Level 1C data by the radiative transfer model-based SEN2COR atmospheric correction processor (version 2.5.5, European Space Agency, Paris, France), and were resampled by Sentinel-2 Toolbox in SNAP. Spectral indices were strongly related to reflectance, and were useful in volume mapping, especially some with red edge bands (band 5, 6, 7, and 8A) [56,57]. Totally, 26 variables from Sentinel-2 were selected and extracted based on previous findings (Table 2) [58,59].

The ALOS Global Digital Surface Model (AW3D30) used in this study was a global dataset generated from L band SAR images collected using the ALOS from 2006 to 2011 (Figure 1e). The data were download from the Japan Aerospace Exploration Agency to extract topographic indices from previous researches by Spatial Analyst of ArcGIS software (version 10.0, ESRI, RedLands, CA, USA) [60,61]. All remote sensing variables were re-projected into UTM Zone 52 WGS84, and then resampled to the 30 m pixel size by ArcGIS.

#### *2.3. Support Vector Machine for Regression Kriging (SVRK) and Modeling Evaluation*

The pairwise Pearson's product-moment correlation analysis was operated to determine predictor variables from multi-sensor indices. It consisted of two steps: the selection of variables which were significantly related to field-measured volume (*p* < 0.05) as candidates; the disposal of candidates that were collinear (*r* ≥ 0.8), except the one that had the largest correlation coefficient with volume [62]. Those analyses were performed in SPSS software (version 21.0, IBM, Armonk, NY, USA).

The SVRK model built in this study is the extension of SVR, which integrated SVR prediction and estimation of the residuals by ordinary kriging using Equation (2). SVRK considers spatial parametric non-stationarity with the effects of multi-sensor predictors derived from the benefits of SVR. It also added the spatial dependence of the residuals interpolated through ordinary kriging to the estimated trend, as part of the spatial autocorrelation:

$$V\_{SVRK} = V\_{SVR} + R\_{OK} \tag{2}$$

where: *VSVRK*, *VSVR* are predication of stand volume based on SVRK and SVR, respectively; *ROK* is the estimated residuals of volume from the SVR prediction.

The implementation of SVRK includes two steps, as shown in Figure 3. SVR is firstly used to model the relationship between stand volume and multi-sensor predictors, as a non-linear machine learning method. It uses kernel functions to project the training data onto a new hyperspace where complex non-linear patterns can be simply illustrated (Figure 3a) [63,64]. The optimal hyperspace, constructed by SVR, fits training data and predicts with minimal empirical risk [65]. The SMO (sequential minimal optimization) algorithm is used to solve the quadratic programming optimization problem step-by-step. It updates the SVR function, as shown in Equation 3, to reflect the new values until the Lagrange multipliers converged [66]:

$$f(\mathbf{x}) = \sum\_{k=1}^{n} (\alpha\_k - \alpha\_k^\*) \mathbf{K} \{ \mathbf{x}\_{k\prime} | \mathbf{x}\_j \} + b$$

where *x* is a vector of the input predictors, *f*(*x*) . is an optimal function developed by SVR, *b* is a constant threshold, *K(xk, xj)* is the radial basis function (RBF) kernel with the best bandwidth parameter σ, and α*<sup>k</sup>* are α<sup>∗</sup> *<sup>k</sup>* the weights (Lagrange multipliers) with the constraints given in Equation 4.

$$\begin{cases} \begin{array}{l} \Sigma\_{k=1}^{n} \left( \alpha\_{k} + \alpha\_{k}^{\*} \right) = 0 \\ 0 \le \alpha\_{k'} \; \alpha\_{k}^{\*} \le \mathbb{C} \end{array} \end{cases} \tag{4}$$

where *C* is the regularization parameter for balancing between the training error and model complexity.

**Figure 3.** Illustration of support vector machine for regression kriging (SVRK). It includes (**a**) the first step as the support vector machine for regression (SVR) algorithm and (**b**) the second step as the sum of SVR and residuals interpolated by ordinary kriging. I, Imean, and Isd are raw, mean and standard deviation values, respectively.

In this study, SVR was conducted in WEKA software (version 3.8, The University of Waikato, Hamilton, New Zealand). Parameters of SVR, *C,* and σ, were selected by the smallest root mean square error (RMSE) based on field-measured volume in training dataset (Figure 1a). In order to determine the importance of multi-sensor predictors on volume mapping, the training data was standardized (Figure 3a), and then the increases in RMSEs were calculated as the predictors were excluded one by one from the SVR model.

At the second step, the residuals resulting from SVR are estimated using the ordinary kriging approach (*ROK*) (Figure 3b). Ordinary kriging, a widely used geostatistical technique, generates an optimal unbiased estimation by the semivariogram [67]. The semivariogram can be modeled by spherical, exponential, and Gaussian functions with three parameters—nugget, range and sill [68]. The nugget is an observation error, and sill is the magnitude of spatial autocorrelation [69]. Thus, the stronger spatial autocorrelation is denoted by the larger value of sill relative to nugget [69]. The range parameter shows on which distance the spatial autocorrelation does not influence any more [69]. The Kolmogorov-Smirnov test (K-S) was used to examine the distribution of residuals based on the stationarity assumption of ordinary kriging. The interpolation of residuals by ordinary kriging was conducted in ArcGIS with the smallest RMSE. Finally, volume prediction by SVRK (*VSVRK*) were acquired as the sum of *VSVR* and *ROK*.

The validation set (Figure 1a) was used to test the performance of volume mapping by SVRK based on the mean error (ME), RMSE, and correlation coefficient between the measured and predicted parameters (*r*) [70]. In order to better estimate accuracy, the mean measured value of stand volume (146.1 m3/ha in Figure 2a) was applied to divide the ME and RMSE. The relative improvement (RI) based on RMSE of SVRK over SVR was used as another index for accuracy evaluation [25].

#### **3. Results**

#### *3.1. Relationship between Field-Measured Volume and Remote Sensing Variables*

In total, 31 variables were significantly related to stand volume (*p* < 0.05) (Table 3), including seven from SAR, 20 from MSI, and four from DSM. The backscatters from different wavelengths all had strongly positive correlations with volume, while variables from ALOS-2 had the closer relationship than that from Sentinel-1. The backscatter from HH was more sensitive to volume than that from HH. Among 10 kinds of texture features from Sentinel-1, only the GLCM mean and variation were actively related to volume. In other words, the increasing texture regularity and variety of VV and VH backscatters indicated the growth of stand volume. It was shown that backscatter texture from VV was more relevant to volume than that from VH.

**Table 3.** Related variables and predictors derived from multi-sensor satellite data for stand volume mapping. \* denotes significance with a *p*-value of the t-test being below 0.05; \*\* denotes strong significance with a *p*-value below 0.01.


As for Sentinel-2 variables, the reflectance of B2–B5, B11, and B12 as well as TCW were negatively related to volume, while the other 13 variables represented the positive correlation. All Sentinel-2 volume-related variables displayed the strong correlation (*p* < 0.01), excluding NLI5 (*p* < 0.05). The vegetation indices that were calculated by characteristic red-edge bands of Sentinel-2 closely connected with volume. Variables from Sentinel-2 had similar performances with that from ALOS-2, which showed the greater sensitivity to volume than Sentinel-1 indices.

All four topographic indicators from ALOS DSM showed the strongly positive influence on the increase of volume. It was indicated that variables from DSM was distinguished in the correlation analysis with volume than that from MSI and SAR. Above all, elevation, ALOS-2 backscatters, the texture features of VV channel of Sentinel-1, and the vegetation indices from Sentinel-2 were comparatively vital for stand volume prediction.

#### *3.2. Modeling Forest Volume by SVRK*

#### 3.2.1. Support Vector Machine for Regression (SVR) Modeling for Volume Mapping

To degrade the redundancy, 15 variables that had *r* values of the correlation analysis among predictor candidates above 0.8 were disposed [62]. The predictors involved in modeling were the following 16 list in Table 3. After standardization of training data, the optimal SVR model was built by *C* and σ setting as 1000 and 0.01, respectively, with the minimum RMSE being 40.58 m3/ha. Based on the magnitude of increase in RMSEs (Figure 4), the SVR model showed topographic indicators as the most important predictor for explaining the spatial variations of stand volume, followed by ALOS-2 backscatters, Sentinel-2 indices and texture features of VV channel from Sentinel-1. The VV backscatter from Sentinel-1 was marginal in volume prediction by SVR.

**Figure 4.** Variable importance shown by increase in the root mean square errors (RMSEs) of SVR models after excluding a predictor.

By the optimal SVR model, the predicted values of stand volume in the study area ranged from 5.37 to 523.84 m3/ha, with the mean and SD of 150.26 and 26.04 m3/ha, respectively (Figure 5a). Predicted values were divided into six levels by intervals of field-measured volume values in Figure 2a. The map depicted that the high-altitude region (Figure 1e) was the large forest volume area, with values ranging from 195.51 to 523.94 m3/ha. Zones with small values of volume (5.37 to 94.40 m3/ha) were located close to the non-forest area. Among six levels of stand volume, the smallest and largest occupy the minority of the study area. It was revealed that the SVR model overestimated the small values, and underestimated the large volume, compared to field-measured data.

**Figure 5.** Stand volume mapping by SVR (**a**), and its residuals interpolated by ordinary kriging (%, divided by the mean measured value of stand volume) (**b**) as well as final prediction by SVRK (**c**).

3.2.2. Integration of SVR Prediction and its Residuals by Ordinary Kriging

The residuals were calculated by the field-measured volume and SVR predicted values based on training data. The result of K-S showed that volume residuals from the SVR model possessed a normal distribution (*p* < 0.05), which could be used to calculate experimental semivariograms for ordinary kriging interpolation (Figure 6). Nugget values of spherical, Gaussian, and exponential models were 1478.7, 1509.1, and 1500.5, respectively. Range values were 16.97, 11.44, and 3.16 km, respectively. Sill values were 1689.04, 1673.94, and 1550.52, respectively. The strongest spatial autocorrelation is shown in the spherical model with the largest values of sill relative to nugget. While, the exponential model of ordinary kriging in Figure 6c was chosen to interpolate residuals from SVR with smallest RMSE 39.18 m3/ha. Based on Equation (2), the SVRK model was built.

By the optimal ordinary kriging model (Figure 6c), the distribution of volume residuals from the SVR model were obtained (Figure 5b). The interpolated values of volume residuals ranged from –29.29 to 45.85% (–42.79 to 66.99 m3/ha), with the mean and SD of 0.74 and 14.56 m3/ha, respectively. It was demonstrated that the overestimation of small volume values was located in the western and southern parts of the study area with residuals ranging from –29.29% to –20%. While the SVR model underestimated large volume values in the northern part of the study area, and residuals were from 40.01% to 45.85%.

**Figure 6.** Experimental variograms and fitted models of residuals from SVR by (**a**) spherical, (**b**) Gaussian, and (**c**) exponential models.

#### *3.3. Models Assessment and Volume Mapping*

Table 4 presented the accuracy of the SVR and SVRK models for estimating volume of the validation set (n = 464). The comparison of SVR and SVRK models demonstrated that additional prediction of residuals by ordinary kriging as the spatial autocorrelation, was more accurate than only considering influences of predictor variables from multi-sensor satellite data (Figure 7). It was indicated by ME values that both two models overestimated stand volume. SVRK remarkably improved accuracy of volume prediction over SVR by 9% (3.77 m3/ha) based on RMSE values.

**Table 4.** Accuracy assessment of stand volume modeling based on independent validation data.

**Figure 7.** Scatter plots of predicted versus observed volume from validation data based on SVR (**a**) and SVRK (**b**) models.

The distribution of forest volume based on the SVRK model was acquired by combing the Figure 5a,b as the Figure 5c. By the optimal SVRK model, the predicted values of stand volume in the study area ranged from 0.18 to 532.83 m3/ha, with the mean and SD of 150.99 and 30.83 m3/ha, respectively (Figure 5c). Based on SVRK mapping, the northern part of the study area with high altitude had the largest volume values ranging from 195.51 to 532.83 m3/ha. In the south with low altitude and nearby the non-forest area, the smallest volume values ranged from 0.18 to 94.40 m3/ha. The map showed different distribution of stand volume with the SVR result, while values remained similar (Figure 5a,c). The six levels of stand volume in the SVRK map covered relatively equal areas than that in the SVR result, especially the largest (<sup>≥</sup> 195.91 m3/ha). It was illustrated that the error, which was caused by the SVR model with the overestimation of small values and underestimation of large volume, was reduced. Forest volume of the SVRK map showed the greater spatial variation than that of the SVR. Namely, combining interpolation values of residuals, the spatial distribution of forest volume was much closer to the measured data (SD = 56.2 m3/ha).

#### **4. Discussion**

#### *4.1. Multi-Sensor Satellite Predictors of Forest Volume Mapping*

The role of multi-sensor variables on volume mapping was revealed by correlation coefficients (Table 3) and importance (Figure 4). SAR was able to penetrate forest canopy to a certain depth, and related to roughness and water content of vegetation [71], so that its variables were valuable for volume prediction. The elevation as a proxy of InSAR height, was dominant in volume prediction of this study. It was support from previous findings that InSAR height and its slope parameter were directly proportional to volume [72,73]. It was found that HV was more contributive to stand volume prediction than HH and VV channel; yet, VH backscatter was not significantly related. It was owing to the stronger sensitivity of HV backscatter to the forest growth stage than the HH and VV polarizations [74,75]. All backscatters showed positive relationships with volume. This was likely due to an increase in the volume of scattering with the growth of trees [76].

The penetrability was weaker with shorter wavelength. It resulted in the weaker capability of C band SAR for volume prediction than that of L band SAR according to our findings. It also revealed saturation problems of backscatters. The measured forest volume values in the study area were partially above 200 Mg/ha (Figure 2a), which were larger than the common saturation value of C band and smaller than that of L band SAR backscatters [48,77]. The results indicated that texture features of SAR data were much more helpful than original backscatters to forest volume prediction, which was consistent with existing researches [52,78]. However, textural indices from Sentinel-1 was marginal in this study for volume mapping compared to the previous finding [79]. It was resulted from the decrease in the heterogeneity by texture analysis and large variations of stand volume in the study area.

As optical sensor data, MSI variables, i.e., reflectance and spectral indices, were powerful for the retrieval of horizontal forest structures such as vegetation types, canopy cover and DBH [80,81]. Results revealed that the short-wave infrared (SWIR) band was the highly ranked variable for predicting stand volume. It was explained by the closer relationship between SWIR spectral band and vegetation properties, i.e., canopy biomass and water content, compared to other electromagnetic spectrum regions [82]. In line with existing studies, reflectance of optical bands and spectral indices were quite helpful in volume prediction [83,84], whereas, the role of red-edge bands and their vegetation indices in this study was minor than that in previous researches [58,85]. This may be caused by the diversity of tree species in the study area with different responses to various red-edge bands, the average relationships of which were weaker. In a word, topographic indices from L band InSAR, backscatters of L band SAR, texture features of VV channel from C band SAR and vegetation indices of MSI were recommended for stand volume mapping based on open-access satellite data in the heterogeneous temperate forests.

#### *4.2. SVR versus SVRK*

It was a pioneering study that built the SVRK hybrid model and utilized it to map stand volume. The optimal RBF-kernel SVR model trained in this study as the first step achieved higher accuracy than multiple linear regressions and SVR models with various kernels based on similar multi-sensor satellite data [31,32], while the SVR model in this study was less accurate than that built by ALOS optical and SAR variables [20]. It was attributed to coarser spatial resolution of L band SAR data and the complex composition of tree species in the study area. Moreover, the density of training dataset of this study (928 samples/171450 ha) was quite smaller than that of the reported research (77 samples/83.71 ha).

Results demonstrated that SVRK improved the mapping accuracy by incorporating interpolation values of residuals to SVR models (Figure 5 and Table 4). The value of the accuracy improvement of SVRK was smaller than that of RFK and ANNK for soil carbon prediction as reported [86–88]. It was resulted from the weaker autocorrelation of residuals from SVR compared with that of soil attributes, as well as the smaller sampling density. It was denoted that stand volume was influenced more by multi-sensor variables, and values of nearby sites affected less. The autocorrelation of volume residuals from SVR was weaker than that of biomass errors from RF, while the accuracy improvement of SVRK (RI = 0.09) for volume mapping was much more than that of RFK (RI = 0.07) for biomass prediction [25]. This is due to the higher spatial heterogeneity and the smaller training dataset of this study. The study concluded that SVRK was a promising approach for mapping stand volume with a small training dataset in heterogeneous temperate forests.

#### *4.3. Spatial Variations of Stand Volume and Forest Management*

The spatial distribution of forest volume derived by SVRK with more equal area of each level than the result of the commonly used SVR model was much closer to the measured data (Figures 2 and 5a,c). Whereas, the stand volume map (SD = 30.83 m3/ha) displayed smaller spatial variations than measured data (SD = 56.2 m3/ha). The large variations of measured values of stand volume was resulted from the positively skew distribution with the majority below 300 m3/ha (Figure 2b). The maximum measured volume of 499.84 m3/ ha belonged to a mature *Pinus koraiensis* (Sieb. et Zucc.) dominant natural forest site in the northern part of the study area. The smaller variations of SVRK-derived volume mainly resulted from the coarse mapping resolution as 30 m in this highly heterogeneous forest landscape. The smallest and largest levels of stand volume (< 94.90 and > 195.51 m3/ha) still occupied smaller areas than other four levels. It was illustrated that 30-m multi-sensor data from mosaic L band SAR and InSAR, C band SAR and MSI displayed a saturation problem in detecting small and large values of stand volume.

The volume values of different forest ages were summarized as Table 5 from the SVRK prediction by multi-sensor satellite data. Based on spatial and age variations of stand volume, certain measures can be taken for the sustainable forest management. In young forests, minimum, maximum and mean were all the smallest among five classes, while the variation was largest. The larger values of stand volume in young forests was mainly attributed to the high stand density. With tree growth, more space and resource competition occur among individual trees. Thus, young forests with volume above 195.51 m3/ha should be the critical focus areas, which need thinning (Figure 5c). However, the young forests with volume below 94.90 m3/ha should be enclosed for cultivation. The volume variation of middle-age forests was second-largest. Middle-age forests with volume above 195.51 m3/ha also require thinning, while the increment felling should be conducted in smaller volume areas. Near-mature forests obtained the smallest maximum of stand volume. Management measures in these forests can include the artificial promotion of natural regeneration and beforehand regeneration. Forest manager can selectively cut weak, pest-infested, and diseased trees in mature and over-mature forests.


**Table 5.** Stand volume of forests with different ages in the study area.

#### **5. Conclusions**

Machine learning modeling with remote sensing data combined sample plot data has become a well adopted method to generate spatially explicit estimates of forest parameters. Among that, SVR has achieved wide success in application and has been praised for its ability to deal with small training datasets. A major shortcoming of machine learning is that it ignores the spatial autocorrelation of neighboring observed data. The main objective of the study was to build a hybrid model, i.e., SVRK, which integrated SVR and its residuals by ordinary kriging, based on a small training dataset. Then SVRK was used to map stand volume, the most common forest parameter needed for sustainable forest management at all scales. This study also determined the potential of open-access satellite predictors from multi-frequency SAR data in predicting volume in the heterogeneous temperate forests. As the first exploration of the SVRK modeling, this study provides an informative foundation for decision makers and other researchers on stand volume mapping with limited samples in northeastern China.

Based on the results of this study, the following was concluded:


**Author Contributions:** L.C., C.R., and B.Z. designed this research. L.C. and C.R. conducted field sampling. L.C. performed the experiments, conducted the analysis, and drafted the manuscript. L.C., C.R., B.Z., and Z.W. revised and finalized the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Key Research and Development Project of China (No. 2016YFC0500300), the Jilin Scientific and Technological Development Program (No. 20170301001NY), the funding from Youth Innovation Promotion Association of Chinese Academy of Sciences (No. 2017277, 2012178), and National Earth System Science Data Center of China.

**Acknowledgments:** The authors are grateful to the support from colleagues and local forestry bureau who participated in the field surveys and data collection. We thank the National Earth System Science Data Center (http://www.geodata.cn) for providing geographic information data. This study is supported by the National Key Research and Development Project of China (No. 2016YFC0500300), the Jilin Scientific and Technological Development Program (No. 20170301001NY), funding from Youth Innovation Promotion Association of Chinese Academy of Sciences (No. 2017277, 2012178), and National Earth System Science Data Center of China.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
