*2.1. Study Area*

The study was conducted in Huzhong National Natural Reserve of Great Xing'an Mountains (Figure 1), where primary forests are well protected since all logging has been prohibited since 1958. In this area, post-fire forests follow a natural recovery process and successional trajectory without artificial seeding or plantation. A strict ban on logging of natural forests has been enforced across the Great Xing'an Mountains since April 2014, thus historical fires in the nature reserve will provide valuable insight for understanding how resilient this forest is in responding to wildfire and for refining post-fire managemen<sup>t</sup> practices. The particular fire (Lat/Lon: 122◦50E, 51◦53N) we studied was ignited by lightning on 17 June 2000 and lasted for seven days [45]. It burned an area approximately 7735 ha, most of which was moderate-high severity with high tree mortality rates [26,39].

**Figure 1.** Spatial location of the study area (**<sup>a</sup>**, red dot) in Great Xing'an Mountains (**<sup>a</sup>**, in blue) of northeastern China (**<sup>a</sup>**, in green). The false color Landsat image (**b**, R: TM7, G: TM4, B: TM3) acquired on 24 May 2002 exhibits an overview of two-year post-fire landscape, while the WorldView-2 image (**c**) shows the post-fire landscape after 14 years recovery. The red and blue dots are field plots explored in 2012 and 2013, respectively.

The pre-fire forests were dominated by Dahurian larch (*Larix gmelinii*), a deciduous conifer, and a few deciduous broadleaf species of birch (*Betula platyphylla*) and aspen (*Populus davidiana* and *Populus suaveolens*) [46]. Some evergreen conifers including Scotch pine (*Pinus sylvestris var. mongolica*) and Korean spruce (*Picea koraiensis*) were also sparsely distributed. The understory layer was the primary fuel bed and source of ladder fuels [45], which consisted of evergreen shrubs (e.g., *Pinus pumila*, *Ledum linnaeus* and *Vaccinium vitisidaea*), deciduous shrubs (e.g., *Betulafruticosa* and *Rhododendron dauricum linnaeus*) and some herbaceous plants (e.g., *Chamaenerion angustifolium*, *Carex appendiculata* and *Rubus linnaeus*) that varied with edaphic and topographic conditions [47]. This area has a relatively short growing seasons (150 days) and frost-free period (~80–100 days), with a mean annual temperature of approximately −4.7 ◦C and mean annual precipitation of about 460 mm [26]. The topography is mountainous with elevation ranges from 760 m to 1300 m.

#### *2.2. Field Data Collection*

For field data collection, we focused our attention on two structural attributes of post-fire vegetation, TSA and LAI, which can be easily evaluated in the field. In summers of 2012 and 2013, our field crews collected 70 plots of sapling inventory data for both TSA and LAI (Figure 1), with plot size set to 30 m × 30 m in order to match the spatial resolution of Landsat imagery and to minimize potential scale issues. This will improve the accuracy of linkage between in situ measurement and Landsat spectral properties. All sites were at least 150 m from roads to eliminate edge effect and were selected according to severity and topographic positions. Within each plot, we set three 5 m × 5 m quadrants along the diagonal (~42 m) to survey sapling stems with a height greater than 1.5 m (Figure 2). We did not distinguish between tree species (although most were white birch and larch) as they were not differentiable in VHR imagery due to very similar spectral signature and highly mixed in situ. All saplings were counted based on the number of stems (>1.5 m) because asexual resprouting of white birch is common. TSA was calculated as the mean number of saplings and was normalized to saplings per hectare.

**Figure 2.** Layout of sampling plot. Three 5 m × 5 m quadrates are used for tree sapling abundance survey. Nine fish-eye photos were used for estimation of leaf area index. Three upward digital hemispherical photography (DHP, top-left) pictures show examples of in situ canopy measurement, while tree processed pictures (low-right) show corresponding classification results (see marks) in CAN-EYE software. The black areas in mid-right picture (blue triangle) represent masked dead standing stems in corresponding DHP picture (top-middle).

LAI was measured using a digital hemispherical photography (DHP) system using a Canon EOS 60D Digital Single Lens Reflex camera and a Sigma 8-mm F3.5 EX DG Fisheye lens. We followed the user manual of CAN-EYE software (Version 6.3.3, National Institute of Agronomical Research, Toulouse, France) to calibrate our DHP system and derive the parameters for establishing the

projection functions [48]. Using a 1920 × 1280 pixels parameter for digital photograph, our calibration process generated accurate projection functions that are comparable with results as shown in user manual (see Appendix A). The DHP system was installed to a tripod with a fixed 0.70 m height. Within each sampling plot, we took nine upward digital photographs distributed as shown in Figure 2. Following the recommendation of the user manual, the limit of the circle of interest was set as 60◦ in CAN-EYE software to avoid mixed pixels effects in the blended photograph edges. Other input parameters were set as either default values or derived from photographs automatically. Given this procedure, dead standing stems with large diameter had the potential to influence the LAI estimation, so we manually masked them during the image process in CAN-EYE software to reduce the overestimation. The canopy of surviving trees was retained, however, because their greenness was involved in Landsat spectral signals and not possible to isolate.

#### *2.3. Remote Sensing for Estimating Vegetation Coverage*

#### 2.3.1. Remote Sensing Data Process

The cloud-free WorldView-2 (DigitalGlobe. Inc., Herndon, VA, USA, Figure 1) image was acquired on 1 June 2014 at 2 m spatial resolution including four multispectral bands (Blue 450–510 nm; Green 510–580 nm, Red 705–745 nm, and near infrared (NIR) 860–1040 nm) and 0.5 m spatial resolution for panchromatic band (450–800 nm). We used the Pan Sharpening toolbox in ENVI 5.3 software (Harris Geospatial Solution. Inc., Washington, DC, USA) to perform image fusion, and improved the spatial resolution of multispectral bands to 0.5 m. The WorldView-2 imagery had been geo-spatially projected using the Universal Transverse Mercator (UTM) coordinate system based on WGS-84 ellipsoid with a nominal positioning accuracy of 3.5 m. It is accurately overlapped with Landsat imageries, which have 30-m spatial resolution.

To produce time-synchronous comparisons between Landsat and field survey, the spectral indices were extracted according to the year when the TSA and LAI data were collected in the field. In addition, we used the same field data for the Landsat-derived indices of 2014 in order to form an unbiased comparison with WorldView-2 data. We obtained 8 Landsat-5 TM, Landsat-7 ETM+, and Landsat-8 OLI surface reflectance images (path-row: 122/24 and 121/24, time: 2010, 2012–2014) from the EROS Science Processing Architecture (ESPA) system of the USGS (Table 1). Atmospheric and topographic correction had been carried out for these before distribution. Due to the lack of cloud-free days over the study area during the time periods assessed, we used "cloud masking" to remove cloud and shadow pixels [49], and then carried out image mosaic to fill data-gap for each year. Although the intra-annual Landsat imageries have very close observation dates, we still applied the iteratively re-weighted multivariate alteration detection (IR-MAD) approach to minimize spectral inconsistency among Landsat sensors [50,51]. The Landsat-5 imagery of 2011 was used as the reference imagery for band-by-band radiometric normalization. The mean correlation coefficients (for 6 bands) of four target imageries were each higher than 0.95, indicating good performance of IR-MAD for radiometric normalization.

**Table 1.** Detailed information of Landsat scenes. Invalid pixels are cloud, shadow and null value pixels.


#### 2.3.2. Landsat-Derived Spectral Indices

Vegetation indices developed from remote-sensed Red and NIR bands, such as the Normalized Differenced Vegetation Index (NDVI), are often used as indicators of vegetation recovery in many forest ecosystems [31]. But saturation issues with vegetation indices are widely reported [34,52], especially in areas covered with dense vegetation, which may limit the strength of vegetation indices for detecting structural changes related to forest recovery [30]. Spectral indices derived from the shortwave infrared (SWIR) bands, which are sensitive to water content in vegetation foliage, are proving to be good indicators of post-fire recovery in many forest ecosystems [34,53,54]. In addition, components of Tasseled Cap (TC) transformation (i.e., brightness, greenness and wetness) and their modified spectral indices are also commonly used for monitoring forest disturbance and forest recovery [55]. Different indices may provide different perspectives of forest recovery. Therefore, it is critical to elucidate which Landsat-derived spectral indices have the feasibility to depict post-fire tree sapling attributes in our ecosystem. Such efforts can be helpful to provide an in-depth understanding of what remote sensing indices mean in the field. To comprehensively investigate the performance of Landsat-derived spectral indices on estimating TSA and LAI, we calculated nine commonly-used spectral indices (Table 2) to analyze their explanatory strength.

**Table 2.** Formulas of nine Landsat-derived spectral indices. Spectral Index Abbreviation: NDVI—Normalized Difference of Vegetation Index; EVI—Enhanced Vegetation Index [56]; SAVI—Soil Adjusted Vegetation Index [52]; MSAVI—Modified Soil Adjusted Vegetation Index [52]; NBR—Normalized Burn Ratio [54]; NDMI—Normalized Difference Moisture Index [57]; TCA—Tasseled Cap Angle [55]; TCW—Tasseled Cap Wetness. Note: spectral bands are different between Landsat5/7 and Landsat 8.


#### 2.3.3. Image Textures of WorldView-2 Imagery

The WorldView-2 has limited spectral bands that cannot sufficiently support to calculate much spectral indices, but it offers image textures to improve the land cover mapping [58,59] and estimation of forest structures [60,61], ye<sup>t</sup> its performance in estimating attributes of sapling tree was not verified. Using the panchromatic band of WorldView-2 imagery, we calculated 13 image texture variables, which consist of five occurrence measures and eight co-occurrence measures. The occurrence measures directly use the number of occurrences of each gray level within a given moving window for calculating summaries of statistics (i.e., Range, Mean, Variance, Entropy and Skewness). The co-occurrence measures are second-order statistical calculations involving the spatial relationships among neighboring pixels. The grey level co-occurrence matrix, which is a function of both the angular relationship and distance between a central pixel and its neighboring pixels, was used to calculate eight statistical measures including mean, variance, homogeneity, contrast, dissimilarity, entropy, angular second moment (ASM), and correlation.

Calculation of image texture is influenced by the window size as it determines the number of surrounding pixels. Wood, et al. [62] suggested that a small window size for image texture analysis may best match the spatial scale at which the forest structure varies. We tested three kinds of window size, including 5 × 5, 15 × 15 and 25 × 25, but did not find notable differences for subsequent analysis. We chose a 5 × 5 window size for all texture variable calculations because it matches the spatial size of a single tree sapling and since it was the least costly in computing time. We did not calculate image texture for Landsat imagery because our field plots had identical spatial size with Landsat pixel, which means the measured TSA and LAI can only reflect within-pixel (30 m) variations. All image texture variables were computed using ENVI 5.3 software (Harris Geospatial Solution. Inc., Herndon, VA, USA).

#### 2.3.4. Land Cover Mapping of WorldView-2 Imagery

We used the per-pixel classification method to obtain the land cover information from the 0.5 m WorldView-2 imagery. According to our field knowledge and visual interpretation of the imagery, we used a classification scheme consisting of 9 popular objects in the field during the preliminary analysis. We used the support vector machine (SVM) approach to carry out the land cover classification. The SVM is a supervised classifier that has been widely applied for land cover and land use mapping through the local to global spatial scales [63,64]. Training pixels were selected with very high confidence based on our knowledge. We chose over 3000 training pixels for each land cover type in order to provide sufficient model training and avoid overfitting problem [65]. We calculated the Jeffries–Matusita distance (JMD) to evaluate the separability of training pixels, defined as:

$$JMD\_{xy} \;=\; 2\left(1 - e^{-B}\right),$$

where *B* is the Bhattacharyya distance, which is usually used for indicating dissimilarity between two distributions:

$$B\_{\pm} = \frac{1}{8}(\mathbf{x} - \mathbf{y})^{\mathbf{t}} \left(\frac{\sum \mathbf{x} + \sum \mathbf{y}}{2}\right)^{-1} (\mathbf{x} - \mathbf{y}) + \frac{1}{2} \ln \left(\frac{\left|\frac{\sum \mathbf{x} + \sum \mathbf{y}}{2}\right|}{\sqrt{|\sum \mathbf{x}| \times |\sum \mathbf{y}|}}\right)$$

In this equation, *x* and *y* are two vectors of spectral and textural signatures and ∑ *x* and ∑ *y* are covariance matrix of *x* and *y* respectively [66]. The JMD will be close to 2 when spectral and textural signatures of two classes are completely different (high separability), and close to 0 when spectral and textural signatures are identical (low separability).

Classification accuracy was evaluated based on additional pixels (~1000 pixels for each class) that were independently selected from the training pixels. Using the confusion matrix approach [67], we calculated the Cohen's Kappa coefficient to evaluate overall classification accuracy, and used the omission error and commission error to evaluate accuracy of individual class. The ability to distinguish tree saplings from the surviving trees was desired to refine our research. However, we found low separability between these two classes (Table 3). We also found low separability between grasslands and shrublands and therefore combined those classes for the further analysis (but reported accuracies evaluated based on the nine-class scheme).

## *2.4. Statistical Analysis*

#### 2.4.1. Compare Performance of Landsat and WorldView-2 on Predicting TSA and LAI

To match the spatial size of field data, we conducted spatial aggregation to obtain statistical measures for information derived from WorldView-2 data. For the land cover map, we used a moving window (60 × 60 pixels), whose spatial size is consistent with spatial resolution of Landsat imagery, to compute the percentage of pixels that were classified as the canopy trees (i.e., tree saplings and mature trees), denoted as PPCT. For image texture, we conducted spatial aggregation to obtain the mean value for each texture variable. Given high collinearity among image texture variables, pairwise Pearson correlation coefficients (r) calculated from the "Hmisc" package in R 3.4.1 (R Development Core Team 2017, Boston, MA, USA) were used as the criteria based on randomly sampled 300 pixels.

We used a forward stepwise method combined with variance inflation factors (VIF) functions in the "car" package to drop variables with high collinearity by applying a stringent threshold of 3, following recommendation of Zuur, et al. [68]. Landsat spectral indices often show collinearity as their computation may use the same spectral bands. Because of this we only retained two co-occurrence texture variables (i.e., mean and ASM) and one occurrence texture variable (i.e., Entropy) as they represented low correlations (|r| < 0.50). We retained all spectral indices for the further analysis to allow us to investigate the usage of a spectral index in monitoring post-fire vegetation dynamics.

The presence of spatial autocorrelation in dependent variables may violate the assumption that all observations are independent, which will inflate the significance and affect the coefficient estimates. Using functions in the "ape" package [69], we calculated the Moran's I index to examine spatial autocorrelations in our measured TSA and LAI respectively. Moran's I is similar to a correlation coefficient ranging between −1 and 1. Higher positive values indicate greater similarity, while lower negative values indicate stronger dissimilarity. Autocorrelation was found significant at α < 0.05 level, but very weak for both TSA (Moran's I = 0.18) and LAI (Moran's I = 0.20). It suggested our field observations on TSA and LAI have low spatial dependence, and the subsequent analysis will not be strongly influenced by spatial autocorrelation.

The coefficient of determination (R2) was used to evaluate the explanatory power of remotely sensed variables for variations of field-measured TSA and LAI. Because the LAI is partially determined by tree density, it is not surprising that our LAI data represented very high correlation (r = 0.783, *p* < 0.01) with TSA data. We retained both variables for the remote sensing analysis as they reflect different perspectives of forest structure and may be sensitive to different spectral bands. Before parameterizing the regression model, we applied the Shapiro-Wilk test and histogram approach to detect whether the normality assumption is violated [70]. If the null hypothesis (i.e., samples came from a normally distributed population) was rejected at α < 0.05 level, we applied logarithm transformation to mitigate skew. Homogeneity of variance was evaluated using the Fligner–Killeen test [71]. The test results indicated that the null hypothesis (i.e., all populations variances are equal) was not violated in any regression model. Outliers (Cook's distance greater than four times of the mean) were detected and removed. The Landsat indices and WorldView-2 variables with the highest R2 were used to predict TSA and LAI through the post-fire landscape.

#### 2.4.2. Evaluation of Relative Importance Using Random Forest Model

We used the random forest (RF) model in the "randomForest" package to evaluate the relative importance of spatial controls on determining TSA and LAI. The RF model is a machine learning algorithm with advantages for dealing with nonlinear relationships, multi-collinearity, and complex interactions without imposing assumptions on data distribution [72]. Based on a bootstrap subsampling (bagging) scheme, the RF model can generate multiple regression trees with low variance that can be combined for an accurate prediction. The best split of node is chosen from the random subsets of predictor variables to ensure that all predictors are tested. For each individual tree, the remaining (i.e., out-of-bag, OOB) data that not drawn into training subset is used for unbiased model validation [73]. We used the variance explained (R2) to evaluate the goodness of fit of RF model for training data. It is calculated based on the internal OOB error rate in terms of mean of squared residuals (MSE):

$$MSE = n^{-1} \sum\_{1}^{n} \left( y\_i - \hat{y}\_i^{OOB} \right)^2$$

where *y*<sup>ˆ</sup>*OOB i*is the average of the OOB predictions for observation *i*. The R<sup>2</sup> is calculated as:

$$\mathbb{R}^2 = \frac{MSE}{\delta\_y^2}$$

The most important variable will cause highest degradation on model fit when it is omitted.

Predictors were spatial variables representing gradients of burn severity, topography and understory vegetation abundance (Table 3). Burn severity was evaluated based on a quadratic correlation model (R<sup>2</sup> = 0.856), which was developed from the difference between the normalized burned ratio (dNBR) and field-sampled severity measures [26]. We applied a 30 m resolution digital elevation model (DEM) to generate topographic variables reflecting elevation, topographic relief, topographic wetness, and the total solar radiation of growing season (April to October). Understory vegetation abundance was derived from WorldView-2 land cover map. Similar to PPCT, we calculated the coverage of understory vegetation (i.e., shrubland and grassland) using a 30 m × 30 m moving window aggregation approach. We also found that the shadow pixels largely reflected mature trees in our WorldView-2 imagery. We used the coverage of tree shadow as a surrogate to represent the number of surviving mature trees post-fire. This is similar to Berner, et al. [74], who applied the tree shadow coverage as a surrogate of tree biomass. To understand the edge effects of unburned areas on post-fire recovery, we also calculated the nearest distance to the edge of unburned area as a predictor. The remote sensed TSA and LAI were used as the response variables for two RF models respectively. We generated random samplings based on a 150 m sampling space to balance the spatial autocorrelation and the maximum sampling number (~500 points). To reduce the risk of stochastic errors and create stable model outputs, we carried out 50 RF modeling trials independently and used the average value as the final result.


**Table 3.** Description of predictors used in random forest models.
