2.2.1. Field Data

The field campaign was carried out from the end of August to the beginning of October in 2017. The stratified sampling design was adopted. Non-forest areas were masked out and the distribution of sampling plots was randomly generated in forest areas, while the plot sites which were impossible to access were replaced by the nearest sites. Following the national guidelines for forest resource survey [73], ten teams took part in the field campaign and collected measured data under the same protocol. A total of 1803 30 m by 30 m plots were located and sampled (Figure 1). At each sample site, tree species, diameter at breast height (DBH, the diameter at 1.3 m from the ground), tree height, soil depth and types, as well as the number of trees, were measured and recorded. Overhead photos of the canopy by fisheye lens were taken at the center of each sample site. The measured parameters were acquired based on field sample sites and processes shown in Table 1. Stand volume of each sample

was estimated by DBH and tree height according to the National Standard of China: Tree volume tables (LY/T 1353–1999) [74]. Canopy closure was estimated from fisheye photos after processing of circle clipping, binary conversion, and raster calculation [75]. Soil fertility was assessed by the depth weighted according to soil types as shown in Table 1. The six types of soils in the study area were assigned three different weights based on previous studies of their organic carbon density [76,77]. The descriptive statistics of field-based forest parameters were shown in Table 2. The 1803 sample sites were randomly split into training (n = 1202) and validation (n = 601) sets (Figure 1) for modeling and validating the spatial variation of forest parameters, respectively.

**Figure 1.** The figure illustrates the shape of the Changbai Mountain National Nature Reserve (CMNNR), and locations of field sample sites, and employed satellite remote sensing data derived from Advanced Land Observing Satellite (ALOS), ALOS-2, Sentinel-1, and Sentinel-2 series.


**Table 1.** Field measurements and processing to acquire measured forest parameters.

**Table 2.** Descriptive statistics of field-measured forest parameters.


#### 2.2.2. Remote Sensing Data

In this study, 18 predictor variables related to forest parameters were selected and extracted from multi-sensors imagery (Tables 3 and 4) [12,78–80]. ALOS-2 PALSAR-2 yearly mosaic image of 2017 was masked and converted to gamma naught values in decibel unit (dB) from 16-bit digital number (DN) (Table 4) using the following equation [81]:

$$
\gamma^0 = \ ^{10}\text{Log}\_{10}\text{[DN}^2\text{]} - 83\tag{1}
$$

where γ0 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.

Monthly mosaic predictors from C band SAR were generated from seven Sentinel-1 Ground Range Detected images by masking and mosaic using Google Earth Engine (GEE). Those data were pre-processed by thermal noise removal, radiometric calibration, and terrain correction stored in dB via log scaling [82]. The cloud-free Sentinel-2B Level-1C images acquired on 25 September 2017 were downloaded from the Copernicus Sentinel Scientific Data Hub (https://scihub.copernicus.eu/) to extract vegetation and soil indices, as well as biophysical variables. Previous studies explored numerous Sentinel-2 spectral indices. They found that four vegetation indices and two soil indicators were useful in modelling forest age and soil fertility, respectively [78,80,83–85]. The MSI data had 13 spectral bands with 10 m (bands 2–4, 8), 20 m (band 5–7, 8a, 11–12), and 60 m (band 1, 9–10) spatial resolutions. Bottom-of-atmosphere-corrected reflectance (Level-2A) images were atmospherically corrected from the Level-1C data by SEN2COR atmospheric correction processor based on the radiative transfer model (version 2.5.5, European Space Agency, Paris, France). Then, eight predictors were acquired by resampling, band math, biophysical processor, and mosaic from Level-2A images based on SNAP software (version 6.0, European Space Agency).

The DSM data from ALOS as AW3D30 were download from Japan Aerospace Exploration Agency (https:ww.eroc.jaxa.jp/ALOS/en/aw3d30/data/index/htm) to extract topographic indices based on Spatial Analyst of ArcGIS software (version 10.0, ESRI, RedLands, CA, USA). All predictor data were re-projected into UTM Zone 52 WGS84, and then resampled to the 30 m pixel size by ArcGIS.


**Table 3.** The ALOS-2, Sentinel-1, Sentinel-2, and digital surface model (DSM) data used in this study.



#### *2.3. Assessment of Forest Conditions*

2.3.1. Spatial Modeling of Canopy Closure and Stand Density by Statistical Regressions

Stand density is a prominent component of forest structure, which governs elemental processing and retention, competition, and habitat suitability [86,87]. Canopy closure is the proportion of the sky hemisphere occupied by tree crowns when viewed from a single ground point [88,89]. It is closely associated with understory light and has wide-reaching e ffects for ecological processes in forests [90]. Whereas FVC is defined as the percentage of the forest area covered by the vertical projection of trees [91,92]. LAI is one half of the total green leaf area per unit ground surface area [93]. FVC and LAI are critical biodiversity variables as recognized by international organizations such as Global Climate Observation System and Global Terrestrial Observation System [35]. The generalized linear correlations were discovered between stand density and LAI, canopy closure and FVC, and LAI and FVC in previous studies [42,94,95]. Thus, it was assumed in this study that canopy closure and stand density could be modeled by generalized linear regressions based on FVC and LAI from Sentinel-2 images. In this study, five types of regression models (linear, quadratic polynomial, power, exponential, and logarithmic) were built. LAI or FVC, and canopy closure or stand density, were used as the input variables to derive the empirical parameters for the models. This study selected the model with the largest value of coe fficient of determination ( *R*2) to map the canopy closure and stand density.

2.3.2. Spatial Modeling of Stand Volume and Forest Age by Random Forests

Firstly, a semi-physical simple water cloud model (WCM) was used for the investigation of the relationship between stand volume and backscatters (HH and HV) derived from ALOS-2 data. The prerequisite assumption of WCM was that the dielectric constant of dry vegetation matter was much smaller than that of the water content of vegetation, and almost all volume backscatters were composed of air in the vegetation canopy [96]. Therefore, WCM was developed assuming that the canopy "cloud", called the water cloud, contained identical water droplets showed the random distribution within the canopy [55]. In this study, the WCM was adopted for the initial exploration [97], which was written as Equation (2):

$$\boldsymbol{\gamma}\_f^0 = \boldsymbol{\gamma}\_\mathcal{S}^0 \mathbf{e}^{-\delta SV} + \boldsymbol{\gamma}\_v^0 (\mathbf{1} - \mathbf{e}^{-\delta SV}) \tag{2}$$

where γ0 *f* is the backscatter from the forest, as the gamma naught value of HH or HV (dB); γ0 *g* is the direct backscatter from the forest floor through gaps in the canopy (dB); γ0 *v* is the volume backscatter from an opaque canopy without gaps (dB); *SV* is stand volume (m<sup>3</sup>/ha); and δ is the extinction coe fficient.

**Figure 2.** Steps of spatial modeling on stand volume and forest age by random forests and that on soil fertility by random forest kriging.

Then, RF was used to model the spatial distribution of stand volume with predictors of HH and HV. The RF was an ensemble of decision trees, which was created by a subset of training sample through replacement as a bagging approach [98]. Each decision tree was independently developed without any pruning and each node was divided using a user-defined number of features selected at random [99]. By producing the forest up to a user-defined number of trees, RF creates trees with large variance and small bias [98,99]. The abovementioned two user-defined parameters, i.e., numFeatures and numIterations, were selected by the smallest root mean square error (RMSE) in WEKA software (version 3.8, The University of Waikato, Hamilton, New Zealand), and the attribute importance was also estimated [100]. The new unlabeled data were input to evaluate and vote, and the finial prediction was the average of the membership (Figure 2a). Additionally, forest age, as indirect and complex retrieved parameters for remote sensing techniques, was also modeled by RF with multi-sensor predictors (Table 4 and Figure 3).

**Figure 3.** The flowchart for spatial modeling of parameters and application on forest condition assessment.

2.3.3. Spatial Modeling of Soil Fertility by Random Forest Kriging

Random forest kriging (RFK) was the extension of RF, which integrated RF prediction values and estimation of the residuals by ordinary kriging (OK) using Equation (3) [101]. It considered spatial parametric non-stationarity with the effects of environmental variables derived from the benefits of RF [102,103]. RFK also added the spatial dependence of the residuals interpolated through OK to the estimated trend, as part of the spatial autocorrelation. RFK has been conducted in soil attribute mapping and had a greater accuracy than RF [104,105]. Its implementation included two steps (Figure 2). First, RF was used to model the relationship between soil fertility and multi-sensor predictors. Second, the result of RFK was predicted as the sum of the RF result and its residuals interpolated by the OK approach using Equation (4).

$$SF\_{RF} = \ S F\_{RF} + R\_{OK} \tag{3}$$

$$\begin{cases} \begin{array}{l} \mathcal{R}\_{OK} = \sum\_{i=1}^{n} \omega\_{i} \mathcal{R}\_{i} \\ \sum\_{i=1}^{n} \omega\_{i} = 1 \end{array} \end{cases} \tag{4}$$

where *SFRFK*, *SFRF* are predication of soil fertility based on RFK and RF, respectively; *ROK* is the estimated residuals of soil fertility from RF models; *Ri* is the residuals of soil fertility from RF models at a measured sample *i*; *wi* is the weight estimated by the stationary OK system as an error variance model to minimize the error from the semivariogram modeling [106]; and *n* is the number of measured values within a neighborhood.

#### 2.3.4. Model Evaluation and Forest Condition Assessment

The validation set (Figure 1) was used to test performances of spatial modeling of forest parameters based on the root mean squared error (Equation (5)), mean absolute error (MAE, Equation (6)), mean error (ME, Equation (7)), and correlation coefficient between the measured and predicted parameters (*r*, Equation (8)). In order to better estimate accuracy, the mean measured value of the parameter (*y*) was applied to divide the RMSE, MAE, and ME (Equations (5)–(7)).

$$\text{RMSE} \ = \frac{1}{\overline{y}} \sqrt{\sum\_{1}^{n} \frac{\left(y\_i - \hat{y}\_i\right)^2}{n}} \times 100\% \tag{5}$$

$$\text{MAE} = \frac{1}{\overline{y}} \sum\_{1}^{n} \frac{\left| y\_i - \hat{y}\_i \right|}{n} \times 100\% \tag{6}$$

$$\text{ME} = \frac{1}{\overline{y}} \sum\_{1}^{n} \frac{y\_i - \hat{y}\_i}{n} \times 100\% \tag{7}$$

$$r = \frac{\sum\_{1}^{n} (y\_i - \overline{y}) \left(\hat{y}\_i - \overline{\hat{y}}\right)}{\sqrt{\sum\_{1}^{n} (y\_{i^-} - \overline{y})} \sqrt{\sum\_{1}^{n} (\hat{y}\_i - \overline{\hat{y}})}} \tag{8}$$

where *yi* is the measured parameter value; *y*ˆ*i* is the predicted parameter value; *y* and *y*ˆ are the average of measured and predicted values of the parameter, respectively; and *n* is 601 in this study. The RMSE and MAE should be as small as possible. The ME should be close to zero, while *r* should be larger.

After that, each map of forest parameters was transformed into the spatial distribution of the parameter score as Equation (9). All parameters were positive indicators for forest condition, except that stand density was considered as a complex indicator. Indeed, excessive or insufficient stand density was harmful to forest conditions [107]. In this study, optimum stand density for forest condition was assigned as the median of measured values (500 tree/ha in Table 2). In other words, stand density below 500 tree/ha was regarded as a positive parameter along with canopy closure, stand volume, forest age, and soil fertility. While stand density above 500 tree/ha was regarded as a negative parameter.

$$Score\_f = \begin{cases} \frac{P\_j - P\_{\min}}{P\_{\max} - P\_{\min}} \times 100, & P \text{ is}, \text{ CC, SV, FA, SF or SD} \le 500\\ \frac{P\_{\max} - P\_j}{P\_{\max} - P\_{\min}} \times 100, & P \text{ is SD} > 500 \end{cases} \tag{9}$$

where *Scorej* is the score of parameter *j*; *Pi*, *P*min, and *P*max are raw data, minimum, and maximum values of spatial modeled parameters, respectively; CC, SV, FA, SF, and SD are canopy closure, stand volume, forest age, soil fertility, and stand density, respectively.

To estimate quantitatively the weight of each parameters, the principal component analysis (PCA) is a common method to use [108,109]. PCA was performed under the factor analysis in SPSS (version 21.0, IBM, Armonk, NY, USA) using Equation (10) by three elements, i.e., coefficients of parameters in linear combinations of different principal components, variance contribution rate of principal components, and normalization of weights. Finally, the forest condition assessment map was generated by the score and weight of each parameter according to Equation (11).

$$\begin{cases} w\_{j} = \frac{\frac{\sum\_{k=1}^{q} \frac{\mathbb{C}\_{jk}}{\sum\_{k=1}^{q} \mathbb{V}\_{k}}}{\sum\_{k=1}^{q} \frac{\mathbb{C}\_{jk}}{\sum\_{k=1}^{q} \frac{\mathbb{C}\_{jk}}{\sum\_{k=1}^{q} \mathbb{V}\_{k}}}} \\ \qquad \frac{\sum\_{j=1}^{5} w\_{j} = 1}{\sum\_{j=1}^{5} \mathbb{V}\_{k} \ge 8} \end{cases} \tag{10}$$

where *wj* is the weight of parameter *j*; *Cjk* is the component matrix value of parameter *j* in component *k*; *Ek* is eigenvalue of component *k*; *Vk* is variation contribution rate of component *k*; and *q* is principal component number.

$$\text{Condition}\_{\text{score}} = \sum \mathbf{\bar{s}}\_{j=1} \mathbf{Score}\_j \times \mathbf{w}\_j \tag{11}$$

where *Conditionscore* is the score of forest condition s; *Scorej* is the score of parameter *j*; and *wj* is the weight of parameter *j*.
