Mapping Forest Canopy Height in Mountainous Areas Using ZiYuan-3 Stereo Images and Landsat Data

Forest canopy height is an important parameter for studying biodiversity and the carbon cycle. A variety of techniques for mapping forest height using remote sensing data have been successfully developed in recent years. However, the demands for forest height mapping in practical applications are often not met, due to the lack of corresponding remote sensing data. In such cases, it would be useful to exploit the latest, cheaper datasets and combine them with free datasets for the mapping of forest canopy height. In this study, we proposed a method that combined ZiYuan-3 (ZY-3) stereo images, Shuttle Radar Topography Mission global 1 arc second data (SRTMGL1), and Landsat 8 Operational Land Imager (OLI) surface reflectance data. The method consisted of three procedures: First, we extracted a digital surface model (DSM) from the ZY-3, using photogrammetry methods and subtracted the SRTMGL1 to obtain a crude canopy height model (CHM). Second, we refined the crude CHM and correlated it with the topographically corrected Landsat 8 surface reflectance data, the vegetation indices, and the forest types through a Random Forest model. Third, we extrapolated the model to the entire study area covered by the Landsat data, and obtained a wall-to-wall forest canopy height product with 30 m × 30 m spatial resolution. The performance of the model was evaluated by the Random Forest’s out-of-bag estimation, which yielded a coefficient of determination (R2) of 0.53 and a root mean square error (RMSE) of 3.28 m. We validated the predicted forest canopy height using the mean forest height measured in the field survey plots. The validation result showed an R2 of 0.62 and a RMSE of 2.64 m.


Introduction
Forest canopy height refers to the mean height of the top surface of a forest's canopy.It is useful in studying biodiversity and the carbon cycle.Many studies have mapped forest canopy height using remote sensing data [1][2][3][4][5][6][7][8][9].The majority of these are based on the light detection and ranging (LiDAR), synthetic aperture radar (SAR) or digital photogrammetry (DP) techniques [5,8,10].
Researchers have demonstrated the feasibility of airborne laser scanning (ALS) in mapping forest canopy height under various conditions [10][11][12].The Geoscience Laser Altimeter System (GLAS) onboard NASA's Ice, Cloud, and land Elevation Satellite (ICESat) produced the most commonly used spaceborne LiDAR data, which can also be used to map forest height [13].Researchers also mapped forest height by combining the LiDAR and auxiliary data, such as the meteorological data and the

Study Area
Our study area is located in the Yangtze River Basin, upstream from the Three Gorges Dam in Hubei Province, Central China (30 • 13'40"-31 • 34'32" N, 110 • 4'54"-111 • 39'11" E) (Figure 1).This is a transitional zone between the mountains and plains with a total area of 11,339.70km 2 .The region is ecologically fragile and dominated by forest ecosystems.It is an important area for ecological protection and biodiversity conservation.The region's ecological significance is the main reason we chose it as our study area.Due to the completion of the Three Gorges Reservoir, the changes in the ecological environment in this region have received considerable attention in the last decade.We expect that the forest canopy height product obtained in this study will serve as a reference for related future studies.The study area we selected is mountainous and the terrain is complex.The altitude ranges from 3 m to 2995 m (Figure 2).The ground slope is between 0° and 77°, with an average of 23.40°.The dominant communities include Pinus massoniana Lamb.forest, Cunninghamia lanceolate (Lamb.)Hook.forest, Quercus variabilis Blume forest, Cyclobalanopsis glauca (Thunb.)Oerst.forest and Pinus armandii Franch.forest.The main vegetation types are coniferous forest, evergreen broad-leaved forest, evergreen and deciduous broad-leaved mixed forest, deciduous broad-leaved forest, coniferous and broad-leaved mixed forest, and shrub [33].

ZY-3 Data
The ZY-3 satellite, launched on 9 January 2012, is China's first civilian stereo mapping satellite.It carries one multi-spectral sensor, and three panchromatic sensors pointing forward, backward, and nadir directions, respectively.All three panchromatic sensors have spectral ranges of 0.50 μm-0.80 μm.The forward and backward sensors have resolutions of 3.5 m × 3.5 m, while the nadir sensor has a resolution of 2.1 m × 2.1 m.The angle between the nadir sensor and the other two sensors is 23.5°.The study area we selected is mountainous and the terrain is complex.The altitude ranges from 3 m to 2995 m (Figure 2).The ground slope is between 0 • and 77 • , with an average of 23.40 • .The dominant communities include Pinus massoniana Lamb.forest, Cunninghamia lanceolate (Lamb.)Hook.forest, Quercus variabilis Blume forest, Cyclobalanopsis glauca (Thunb.)Oerst.forest and Pinus armandii Franch.forest.The main vegetation types are coniferous forest, evergreen broad-leaved forest, evergreen and deciduous broad-leaved mixed forest, deciduous broad-leaved forest, coniferous and broad-leaved mixed forest, and shrub [33].The study area we selected is mountainous and the terrain is complex.The altitude ranges from 3 m to 2995 m (Figure 2).The ground slope is between 0° and 77°, with an average of 23.40°.The dominant communities include Pinus massoniana Lamb.forest, Cunninghamia lanceolate (Lamb.)Hook.forest, Quercus variabilis Blume forest, Cyclobalanopsis glauca (Thunb.)Oerst.forest and Pinus armandii Franch.forest.The main vegetation types are coniferous forest, evergreen broad-leaved forest, evergreen and deciduous broad-leaved mixed forest, deciduous broad-leaved forest, coniferous and broad-leaved mixed forest, and shrub [33].

ZY-3 Data
The ZY-3 satellite, launched on 9 January 2012, is China's first civilian stereo mapping satellite.It carries one multi-spectral sensor, and three panchromatic sensors pointing forward, backward, and nadir directions, respectively.All three panchromatic sensors have spectral ranges of 0.50 μm-0.80 μm.The forward and backward sensors have resolutions of 3.5 m × 3.5 m, while the nadir sensor has a resolution of 2.1 m × 2.1 m.The angle between the nadir sensor and the other two sensors is 23.5°.

ZY-3 Data
The ZY-3 satellite, launched on 9 January 2012, is China's first civilian stereo mapping satellite.It carries one multi-spectral sensor, and three panchromatic sensors pointing forward, backward, and nadir directions, respectively.All three panchromatic sensors have spectral ranges of 0.50 µm-0.80 µm.The forward and backward sensors have resolutions of 3.5 m × 3.5 m, while the nadir sensor has a resolution of 2.1 m × 2.1 m.The angle between the nadir sensor and the other two Forests 2019, 10, 105 4 of 21 sensors is 23.5 • .The multi-spectral sensor consists of four channels: Blue (0.45 µm-0.52 µm), green (0.52 µm-0.59µm), red (0.63 µm-0.69 µm), and near infrared (0.77 µm-0.89µm), with resolutions of 5.8 m × 5.8 m.The positioning error of the ZY-3 is less than 4 m and the height measurement error is less than 3 m if the ground control points (GCPs) are available [34][35][36].The ground swath is about 50 km × 50 km and the digitalizing bit is 10 bits [37].The ZY-3 data used in this study were acquired on 8 October 2014, during the leaf-on season.We ordered the data from the website of the China Centre for Resources Satellite Data and Application (CRESDA, http://www.cresda.com/CN/).As the forward view image contained more mountain shadows, we only used the backward and nadir view panchromatic images to extract the DSM.We used the multi-spectral image to classify the land-cover and extract forest areas.

SRTMGL1 Data
The Shuttle Radar Topography Mission (SRTM) data were collected from 11 February 2000 to 22 February 2000, with a radar interferometry system [38][39][40][41][42].The SRTMGL1 data for China have been available since July 2015 [43].The spatial resolution of this data is about 30 m × 30 m in the study area.Elevation values are stored as 16-bit signed integers in the HGT files.The vertical accuracy of SRTMGL1 in Hubei Province ranges from 1.9 m to 18.6 m, with higher accuracy in flatter areas.In this study, we relied on the SRTMGL1 to determine the elevation of the GCPs and used it as the DTM to calculate the slope and aspect.We downloaded SRTMGL1 Version 3 product from NASA's website (https://earthdata.nasa.gov/).This is a void-filled version, with fill value information included in the ancillary NUM files.
We also used the SRTMGL1 as the DTM to extract the crude CHM in this study.The penetration capability of the SRTM C-radar was one of the main reasons for this choice.There are also several other near-global terrain data available [43].The ASTER Global Digital Elevation Model (GDEM) and ALOS World 3D-30 m DEM (AW3D30) are based on optical stereo images from multiple periods and cannot penetrate the forest, in principle.The terrain data produced by the TerraSAR-X add-on for Digital Elevation Measurements (TanDEM-X) mission are also DSMs because the X-band radar beam can hardly penetrate the canopy [44,45].It was reported that the phase center of the SRTM C-radar is located above the ground in many cases, but specific analysis is required.Forest height, structure, and density all influence penetration depth.Studies have shown that the SRTM radar phase center can still reach the land surface if the forest height is relatively low (e.g., less than 30 m, as was the case in our study area) [39,46].In addition, the SRTM radar beam should penetrate deeper in leaf-off seasons (i.e., February 11 to February 22, 2000), because during these times, the bare land surface is visible even in optical images.If the phase center is close to the land surface, or the penetration depth increases as the height of the forest increases, using the SRTMGL1 as a DTM remains worthwhile.

Landsat 8 OLI Surface Reflectance Data
The Landsat 8 Operational Land Imager (OLI) surface reflectance product is derived from the Landsat 8 Level-1 precision and terrain corrected product (L1TP), using the Landsat 8 Surface Reflectance Code (LaSRC) algorithm [47].The acquisition date of the Landsat 8 data used in this study is 15 September 2013.We used them for the land-cover classification, principal component analysis (PCA), image segmentation, and to calculate the vegetation indices as well as to predict the forest canopy height.We also collected the Landsat 8 data acquired in the non-growing season to help with the classification of forest types.The acquisition date of the non-growing season data is 21 January 2014.We obtained these data from NASA's website (https://earthdata.nasa.gov/).The images acquired in 2013 were contaminated by a few clouds and haze on the southwestern and eastern edges of the study area (Figure 14).The effects of haze were mitigated by the LaSRC algorithm.The clouds were masked using the Level-2 Pixel Quality Assurance band in the Landsat 8 surface reflectance product.The cloud flags in the Level-2 Pixel Quality Assurance band were derived from the Function of Mask (Fmask) algorithm, with the high-confidence cloud pixels dilated [47,48].A few remaining cloud pixels were eliminated after the land-cover classification.The images acquired in 2014 were cloudless but contained more mountain shadows.The shadows in the images were also masked after the land-cover classification.

Field Data
We obtained field survey data from the forestry administration.The survey plots were set near the coordinate grid nodes on topographic maps and were located in intact forests.The plots are square or rectangular and cover an area of 0.667 hectares.Field measurements were undertaken in the leaf-off season of 2013.The diameter at breast height (DBH) of all the trees with a DBH greater than 5 cm in each plot was recorded.Then, three to five trees, with a DBH close to the average, were selected and their heights were derived from angle and distance measurements.Angles were measured using a mechanical clinometer.Distances were measured using a laser range finder.The average height of the selected trees was recorded as the mean height of the forest.The measured mean forest height was used to validate the predicted forest canopy height of this study.

Methods
Figure 3 displays a flow chart of the methods we used in this study.These methods can be divided into four sections.First, we derived a DSM from the ZY-3 stereo images and subtracted the SRTMGL1 to obtain a crude CHM.We then applied a series of refining measures to the crude CHM, including the masking of non-forest areas using land-cover classification results obtained from the ZY-3 multi-spectral image.Second, we correlated the refined CHM with the topographically corrected Landsat 8 surface reflectance data, the vegetation indices, and the forest types using a Random Forest regression model.Third, we extrapolated the model to the entire study area covered by the Landsat 8 data.Finally, we compared our predicted forest canopy height with the mean forest height measured in the field survey plots.

Extraction and Refinement of Crude CHM
We extracted the DSM from the ZY-3 stereo images using the OrthoEngine module of the PCI Geomatica software (PCI Geomatics Enterprises, Inc., Canada).The backward view image was automatically re-sampled to the same resolution as the nadir view image.The module uses the polynomial coefficients, GCPs, and tie points to compute a math model that relates the rows and columns of the matching pixels with ground coordinates and elevations.The polynomial coefficients are provided in the rational polynomial coefficients (RPC) files distributed with the ZY-3 data.We set forty-nine GCPs manually in the vegetation-free areas referencing the ZY-3 multi-spectral image, ESRI's online World Imagery, and Google Earth (Figure 4).The elevation of the GCPs was derived from the SRTMGL1 to minimize the elevation differences between the ZY-3 DSM and SRTMGL1.We collected a total of 208 tie points interactively in the nadir and backward view images.We set the sampling distance of the output DSM to 2 m × 2 m, to avoid the loss of accuracy.The extracted ZY-3 DSM can only represent the elevation of the upper surface of the forest canopy.We used bi-linear interpolated 2 m × 2 m resolution SRTMGL1 data as the DTM and subtracted it from the 2 m × 2 m resolution ZY-3 DSM to obtain a crude CHM (Figure 5).We had removed the filled values from the SRTMGL1 in advance according to the auxiliary NUM files.In the crude CHM, forests with canopy heights larger than 15 m are mostly concentrated in the north and mid-west.This corresponds to the locations of protected areas and scenic spots.Subgraphs (b) and (f), (c) and (g), and (d) and (h) in Figure 6 show the vertical information of different surface objects, such as forest patches and new houses.
Rising water levels after the 2009 completion of the Three Gorges Dam are responsible for the extracted CHM above the Yangtze River.We also found numerous erroneous values in rugged places, especially on the left and right edges of the ZY-3 coverage, which could come from shadows and steep terrain.As Figure 6e shows, both the CHM values above 45 m (the red areas) and below - The extracted ZY-3 DSM can only represent the elevation of the upper surface of the forest canopy.We used bi-linear interpolated 2 m × 2 m resolution SRTMGL1 data as the DTM and subtracted it from the 2 m × 2 m resolution ZY-3 DSM to obtain a crude CHM (Figure 5).We had removed the filled values from the SRTMGL1 in advance according to the auxiliary NUM files.The extracted ZY-3 DSM can only represent the elevation of the upper surface of the forest canopy.We used bi-linear interpolated 2 m × 2 m resolution SRTMGL1 data as the DTM and subtracted it from the 2 m × 2 m resolution ZY-3 DSM to obtain a crude CHM (Figure 5).We had removed the filled values from the SRTMGL1 in advance according to the auxiliary NUM files.In the crude CHM, forests with canopy heights larger than 15 m are mostly concentrated in the north and mid-west.This corresponds to the locations of protected areas and scenic spots.Subgraphs (b) and (f), (c) and (g), and (d) and (h) in Figure 6 show the vertical information of different surface objects, such as forest patches and new houses.
Rising water levels after the 2009 completion of the Three Gorges Dam are responsible for the extracted CHM above the Yangtze River.We also found numerous erroneous values in rugged places, especially on the left and right edges of the ZY-3 coverage, which could come from shadows and steep terrain.As Figure 6e shows, both the CHM values above 45 m (the red areas) and below - In the crude CHM, forests with canopy heights larger than 15 m are mostly concentrated in the north and mid-west.This corresponds to the locations of protected areas and scenic spots.Subgraphs (b) and (f), (c) and (g), and (d) and (h) in Figure 6 show the vertical information of different surface objects, such as forest patches and new houses.
Rising water levels after the 2009 completion of the Three Gorges Dam are responsible for the extracted CHM above the Yangtze River.We also found numerous erroneous values in rugged places, especially on the left and right edges of the ZY-3 coverage, which could come from shadows and steep terrain.As Figure 6e shows, both the CHM values above 45 m (the red areas) and below −15 m (the black areas), were erroneous.They correspond to the areas with shadows and steep terrain in Figure 6a.In addition, resolution differences could also cause errors, which would be mentioned later.The erroneous values in the crude CHM highlighted the need for data refining.
15 m (the black areas), were erroneous.They correspond to the areas with shadows and steep terrain in Figure 6a.In addition, resolution differences could also cause errors, which would be mentioned later.The erroneous values in the crude CHM highlighted the need for data refining.We refined the crude CHM through a series of measures.First, the areas with poor correlation between stereo images were excluded.We achieved this by masking pixels with scores equal to zero.The score layer was generated during the DSM extraction process.Pixels with scores equal to zero are places where the matching between the stereo images failed, and their values were filled based on the surrounding values.Since they were not calculated from the geometric model, they needed to be excluded.
Second, we ortho-rectified the multi-spectral image, acquired simultaneously with the stereo images using the ZY-3 DSM.We then radiometrically corrected it, based on the calibration parameters downloaded from the CRESDA (http://www.cresda.com/CN/).Next, we used the maximum likelihood supervised classification method to classify the multi-spectral image into five categories.Training samples were selected in Google Earth.The number of samples for forest, water body, shadow, farmland, and bare surface were 80, 20, 50, 95, and 90, respectively.After the classification, only the CHM pixels in forest areas were preserved.Third, we averaged the 2 m × 2 m resolution CHM pixels into 30 m × 30 m to align with the pixels of the Landsat.In the process of averaging, we removed the areas with pixel counts less than 5, or standard deviations larger than 1/3 of the mean values in each 30 m × 30 m cell, which correspond to the regions with poor uniformity.These regions are susceptible to the averaging operation.A total of 61.85% of the averaged CHM pixels were excluded in this way.Fourth, we eliminated 2.80% of CHM pixels, with values less than 0 m or larger than 30 m.This threshold was determined based on the information provided by the forestry department.Few forests in this region have mean canopy heights greater than 30 m.Therefore, canopy heights greater than 30 m are more likely to represent abnormal values.
Fifth, the ZY-3 DSM has a finer resolution and contains more topographical details such as hills and depressions than the SRTMGL1.This means that when we subtracted the SRTMGL1 from the ZY-3 DSM, these topographical details would be transferred to the crude CHM and mix with the forest canopy heights.To solve this problem, we used the slope difference layer between the ZY-3 DSM and the SRTMGL1 as a mask.In order to rule out the undulations on the canopy surface, we averaged the ZY-3 DSM to the same resolution as the SRTMGL1 before calculating the slopes and their differences.As Figure 7 shows, we only preserved the CHM pixels with slope differences less than 4.5°.We relied on visual inspections referencing the stereoscopic display of the terrain data to select the slope difference threshold.We found that the areas being masked increased sharply as the threshold decreased.The 4.5° threshold was a compromise option to minimize the interferences from terrain, while retaining enough samples.The masked parts correspond to the rugged areas where the CHM is abnormal.A total of 17.85% of the CHM pixels were excluded in this step.We refined the crude CHM through a series of measures.First, the areas with poor correlation between stereo images were excluded.We achieved this by masking pixels with scores equal to zero.The score layer was generated during the DSM extraction process.Pixels with scores equal to zero are places where the matching between the stereo images failed, and their values were filled based on the surrounding values.Since they were not calculated from the geometric model, they needed to be excluded.
Second, we ortho-rectified the multi-spectral image, acquired simultaneously with the stereo images using the ZY-3 DSM.We then radiometrically corrected it, based on the calibration parameters downloaded from the CRESDA (http://www.cresda.com/CN/).Next, we used the maximum likelihood supervised classification method to classify the multi-spectral image into five categories.Training samples were selected in Google Earth.The number of samples for forest, water body, shadow, farmland, and bare surface were 80, 20, 50, 95, and 90, respectively.After the classification, only the CHM pixels in forest areas were preserved.Third, we averaged the 2 m × 2 m resolution CHM pixels into 30 m × 30 m to align with the pixels of the Landsat.In the process of averaging, we removed the areas with pixel counts less than 5, or standard deviations larger than 1/3 of the mean values in each 30 m × 30 m cell, which correspond to the regions with poor uniformity.These regions are susceptible to the averaging operation.A total of 61.85% of the averaged CHM pixels were excluded in this way.Fourth, we eliminated 2.80% of CHM pixels, with values less than 0 m or larger than 30 m.This threshold was determined based on the information provided by the forestry department.Few forests in this region have mean canopy heights greater than 30 m.Therefore, canopy heights greater than 30 m are more likely to represent abnormal values.
Fifth, the ZY-3 DSM has a finer resolution and contains more topographical details such as hills and depressions than the SRTMGL1.This means that when we subtracted the SRTMGL1 from the ZY-3 DSM, these topographical details would be transferred to the crude CHM and mix with the forest canopy heights.To solve this problem, we used the slope difference layer between the ZY-3 DSM and the SRTMGL1 as a mask.In order to rule out the undulations on the canopy surface, we averaged the ZY-3 DSM to the same resolution as the SRTMGL1 before calculating the slopes and their differences.As Figure 7 shows, we only preserved the CHM pixels with slope differences less than 4.5 • .We relied on visual inspections referencing the stereoscopic display of the terrain data to select the slope difference threshold.We found that the areas being masked increased sharply as the threshold decreased.The 4.5 • threshold was a compromise option to minimize the interferences from terrain, while retaining enough samples.The masked parts correspond to the rugged areas where the CHM is abnormal.A total of 17.85% of the CHM pixels were excluded in this step.Sixth, to further reduce the impact of the terrain, the CHM pixels with slopes larger than 7.5° were removed, accounting for 16.96% of all CHM pixels.This threshold was also a compromise option to suppress terrain interference while maintaining sufficient CHM pixels.In addition, we referenced the filtering methods of the GLAS data in mountainous areas.In those cases, the slope threshold was set to 5° to 7° [3,49].

Preparation of Landsat Data
We used the Level-2 Pixel Quality Assurance band in the Landsat 8 surface reflectance product, in order to assess the per-pixel quality.In each scene, only the pixels flagged as clear land pixels were preserved.We then classified the clear land pixels using the maximum likelihood supervised classification method.Training samples were selected in Google Earth.The number of samples for forest, water body, shadow, cloud, farmland, and bare surface were 200, 45, 120, 30, 52, and 150, respectively.We only used the pixels classified as forest in the subsequent processing and modeling.
We applied the C-correction [50] to correct the effects of illumination and terrain on surface reflectance.It can be calculated as follows: cos(i) = cos(sz) × cos(sl) + sin(sz) × sin(sl) × cos(az − as) (1) where cos(i) is the cosine of the solar incidence angle, sz is solar zenith angle, sl is slope, az is solar azimuth angle, as is aspect, ρλ,t is the spectral reflectance influenced by topography, bλ and mλ are intercept and slope of the linear regression model between ρλ,t and cos(i).cλ is the empirical parameter calculated separately for each spectral band λ. ρλ,n is the topographically corrected spectral reflectance.The slope and aspect data used in Equation (1) were derived from the SRTMGL1.The solar azimuth and zenith data used in Equation (1) and Equation (2) were provided in the Landsat surface reflectance product.Next, we calculated the following four vegetation indices: The normalized difference vegetation index (NDVI), the enhanced vegetation index (EVI), the soil adjusted vegetation index (SAVI), and the modified soil adjusted vegetation index (MSAVI).
We classified the pixels of the forest area into different forest types.The classification of forest types was based on four vegetation index layers derived from the topographically corrected Landsat 8 non-growing season scenes (acquired on January 21, 2014).Again, in this case, we used the maximum likelihood supervised classification method.A total of 100 training samples of deciduous forest, and 100 training samples of evergreen forest, were selected by referencing multi-temporal high spatial resolution images from Google Earth.We also performed a forest type classification based on Sixth, to further reduce the impact of the terrain, the CHM pixels with slopes larger than 7.5 • were removed, accounting for 16.96% of all CHM pixels.This threshold was also a compromise option to suppress terrain interference while maintaining sufficient CHM pixels.In addition, we referenced the filtering methods of the GLAS data in mountainous areas.In those cases, the slope threshold was set to 5 • to 7 • [3,49].

Preparation of Landsat Data
We used the Level-2 Pixel Quality Assurance band in the Landsat 8 surface reflectance product, in order to assess the per-pixel quality.In each scene, only the pixels flagged as clear land pixels were preserved.We then classified the clear land pixels using the maximum likelihood supervised classification method.Training samples were selected in Google Earth.The number of samples for forest, water body, shadow, cloud, farmland, and bare surface were 200, 45, 120, 30, 52, and 150, respectively.We only used the pixels classified as forest in the subsequent processing and modeling.
We applied the C-correction [50] to correct the effects of illumination and terrain on surface reflectance.It can be calculated as follows: cos(i) = cos(sz) × cos(sl) + sin(sz) × sin(sl) × cos(az − as) (1) where cos(i) is the cosine of the solar incidence angle, sz is solar zenith angle, sl is slope, az is solar azimuth angle, as is aspect, ρ λ,t is the spectral reflectance influenced by topography, b λ and m λ are intercept and slope of the linear regression model between ρ λ,t and cos(i).c λ is the empirical parameter calculated separately for each spectral band λ. ρ λ,n is the topographically corrected spectral reflectance.The slope and aspect data used in Equation (1) were derived from the SRTMGL1.The solar azimuth and zenith data used in Equation (1) and Equation (2) were provided in the Landsat surface reflectance product.Next, we calculated the following four vegetation indices: The normalized difference vegetation index (NDVI), the enhanced vegetation index (EVI), the soil adjusted vegetation index (SAVI), and the modified soil adjusted vegetation index (MSAVI).
We classified the pixels of the forest area into different forest types.The classification of forest types was based on four vegetation index layers derived from the topographically corrected Landsat 8 non-growing season scenes (acquired on January 21, 2014).Again, in this case, we used the maximum likelihood supervised classification method.A total of 100 training samples of deciduous forest, and 100 training samples of evergreen forest, were selected by referencing multi-temporal high spatial resolution images from Google Earth.We also performed a forest type classification based on the reflectance bands and NDVI of the topographically corrected Landsat 8 growing season scenes (acquired on September 15, 2013), using the same set of training samples.This was used to fill the void areas caused by shadows in the non-growing season classification result.
We performed PCA on seven reflectance bands and four vegetation indices of the Landsat data.The first and second principal components contributed 93.22% of the accumulative eigenvalue.They were stacked with the forest type layer to serve as the input data of the multi-resolution segmentation.We used the eCognition Developer software (Trimble Germany GmbH, Germany) to implement the multi-resolution segmentation.Then, we calculated the average slope in each segment based on the slope layer derived from the SRTMGL1.The segments with average slopes of no more than 7.5 • were selected to clip the CHM layer.This threshold was also determined by experiments to suppress terrain interferences, while retaining sufficient CHM pixels.A total of 0.5% CHM pixels were excluded in this way.Finally, after excluding a total of 99.96% of crude CHM pixels, we obtained the refined CHM samples for use in subsequent modeling.

Random Forest Modeling and Extrapolation
Statistical modeling and extrapolation were used in this study for two reasons.First, because we discarded most of the crude CHM pixels during the refining process, the refined CHM samples might be very sparse.Thus, we needed to use the statistical model to obtain a wall-to-wall forest canopy height map.Second, the coverage of the study area is much larger than that of the ZY-3.The use of a statistical model is more economically efficient.It will reduce the dependence on the abundance of the ZY-3 data and increase the update frequency of the forest canopy height product.
Random Forest [51] is an ensemble of unpruned decision trees, induced from bootstrap samples.Each tree is grown with randomized subsets of variables and functions by recursively partitioning the dependent variable into less varying subsets.A prediction is made by averaging the predictions of the ensemble.It is robust to noise and has good generalization ability.We employed the Random Forest algorithm to correlate the refined CHM with the topographically corrected Landsat 8 surface reflectance data, the vegetation indices, and the forest types.The SAVI and MSAVI were not used in the regression model as they showed strong correlations with the NDVI and EVI.All the 10 predictor layers (seven reflectance layers, two vegetation index layers, and a forest type layer) were clipped to the ZY-3 coverage.The Random Forest's built-in out-of-bag estimation was used to evaluate the performance of the model.There are two important parameters can be used to tune the model.The mtry is the number of variables tried at each split.The ntree is the number of decision trees.We determined these parameters by experiments.The values that produced the best out-of-bag estimation result were selected.The mtry was set to 4, and the ntree was set to 200 after the experiments.Variable importance was evaluated by the percentage increase in the mean squared error (%IncMSE) and the increase in node purity (IncNodePurity).The %IncMSE is the increase of mean squared error after randomly permuting the values of the variable, averaged over all trees and normalized by the standard deviation of the differences.The IncNodePurity is calculated by adding up the decreases of residual sum of squares when the variable splits nodes, calculated over all trees.
We used the random Forest package in R environment (R Development Core Team, http://www.R-project.org) to implement the Random Forest algorithm.We extrapolated the regression model to the entire study area to obtain a wall-to-wall forest canopy height product.The predicted forest canopy height was independently validated using the mean forest height measured in the field survey plots.

Refined CHM Samples
Figure 8 shows the refined CHM samples.Most of the samples were located in the relatively flat areas between the mountains.Some were at the tops of the mountains.The number of samples was 1,020.Their mean was 9.55 m, close to the reference mean forest height of 8.59 m provided by the forestry department for this area.The spatial resolution of the refined CHM samples was 30 m × 30 m.

Forest Type Classification
Figure 9 shows the forest type classification result.To evaluate the classification accuracy, we randomly set 500 points in the study area and used manual interpretation in Google Earth to determine the reference categories for each point.
The number of points used for assessing the accuracy of classification was determined based on a balance between what was statistically sound and what was practically attainable.Congalton proposed an empirical minimum number of points of 50-100 for each category.He also pointed out that this number should be adjusted according to the relative importance and the inherent variability of each category [52].In our case, a total of 321 points fell into the deciduous forest areas, 143 points fell into the evergreen forest areas, and 36 points fell into the non-forest areas.Non-forest areas (accounting for 7.27% of the total area) got relatively less points.But this was acceptable, since they were mainly composed of waters with small variabilities and were less important relative to forests.At the same time, the workload of manually interpreting 500 points remained attainable.
The confusion matrix showed producer accuracy values of 0.62, 0.94, and 0.78 for non-forest area, deciduous forest, and evergreen forest, respectively.The user accuracy values were 0.89, 0.84, and 0.87, respectively.The total accuracy was 0.86 and the kappa coefficient was 0.73.

Forest Type Classification
Figure 9 shows the forest type classification result.To evaluate the classification accuracy, we randomly set 500 points in the study area and used manual interpretation in Google Earth to determine the reference categories for each point.
The number of points used for assessing the accuracy of classification was determined based on a balance between what was statistically sound and what was practically attainable.Congalton proposed an empirical minimum number of points of 50-100 for each category.He also pointed out that this number should be adjusted according to the relative importance and the inherent variability of each category [52].In our case, a total of 321 points fell into the deciduous forest areas, 143 points fell into the evergreen forest areas, and 36 points fell into the non-forest areas.Non-forest areas (accounting for 7.27% of the total area) got relatively less points.But this was acceptable, since they were mainly composed of waters with small variabilities and were less important relative to forests.At the same time, the workload of manually interpreting 500 points remained attainable.
The confusion matrix showed producer accuracy values of 0.62, 0.94, and 0.78 for non-forest area, deciduous forest, and evergreen forest, respectively.The user accuracy values were 0.89, 0.84, and 0.87, respectively.The total accuracy was 0.86 and the kappa coefficient was 0.73.
It was found that the deciduous forests were widely distributed throughout the study area.They covered a total area of 7,275.02km 2 and had an average distribution-altitude of 923 m.The evergreen forests showed a clustered distribution.They covered a total area of 3,240.54km 2 and had an average distribution-altitude of 917 m.It was found that the deciduous forests were widely distributed throughout the study area.They covered a total area of 7,275.02km 2 and had an average distribution-altitude of 923 m.The evergreen forests showed a clustered distribution.They covered a total area of 3,240.54km 2 and had an average distribution-altitude of 917 m.It was found that the deciduous forests were widely distributed throughout the study area.They covered a total area of 7,275.02km 2 and had an average distribution-altitude of 923 m.The evergreen forests showed a clustered distribution.They covered a total area of 3,240.54km 2 and had an average distribution-altitude of 917 m.We performed multi-resolution segmentation based on the composited layer.Segmentation parameters were set interactively to ensure that different forest types were separated and the segment size was about 1 km.Finally, we obtained 357,628 segments with a mean area of 3.17 hectares.

30 m × 30 m Resolution Forest Canopy Height Mapping
As Figure 11 shows, band5 (Near infrared band) and band3 (Green band) were the most important variables in the Random Forest model.Band5 had the highest %IncMSE (18.84%) and band3 had the Forests 2019, 10, 105 13 of 21 most significant impact on the IncNodePurity.Band1 (Coastal/Aerosol band) and band4 (Red band) also had relatively high levels of %IncMSE, and IncNodePurity, respectively.
The importance of bands 3 (0.53 µm-0.59µm) and 5 (0.85 µm-0.88 µm) have a bio-physical basis.Band3 is related to chlorophyll.A more intense reflectance in this band usually corresponds to younger forests and therefore, corresponds to a lower canopy height.Band5 is associated with biomass, and biomass is highly correlated with forest canopy height.Therefore, we suggest that they deserve closer attention in future research.If studies use multi-spectral data, other than Landsat, they should make the existence of these two bands a prerequisite.
We performed multi-resolution segmentation based on the composited layer.Segmentation parameters were set interactively to ensure that different forest types were separated and the segment size was about 1 km.Finally, we obtained 357,628 segments with a mean area of 3.17 hectares.

30 m × 30 m Resolution Forest Canopy Height Mapping
As Figure 11 shows, band5 (Near infrared band) and band3 (Green band) were the most important variables in the Random Forest model.Band5 had the highest %IncMSE (18.84%) and band3 had the most significant impact on the IncNodePurity.Band1 (Coastal/Aerosol band) and band4 (Red band) also had relatively high levels of %IncMSE, and IncNodePurity, respectively.
The importance of bands 3 (0.53 μm-0.59μm) and 5 (0.85 μm-0.88 μm) have a bio-physical basis.Band3 is related to chlorophyll.A more intense reflectance in this band usually corresponds to younger forests and therefore, corresponds to a lower canopy height.Band5 is associated with biomass, and biomass is highly correlated with forest canopy height.Therefore, we suggest that they deserve closer attention in future research.If studies use multi-spectral data, other than Landsat, they should make the existence of these two bands a prerequisite.We performed a piecewise averaging of all samples to reveal the relationships between predictor variables and refined CHM (Appendix A).It was found that, with the increase of the refined CHM, the reflectance of all bands showed a downward trend, and the NDVI and EVI showed upward trends.Evergreen forests usually accounted for larger proportions in the groups with higher refined CHMs.We used out-of-bag estimation to evaluate the Random Forest regression model.The canopy height of each sample was predicted, based on about 1/3 of the regression trees.The estimation showed a coefficient of determination (R 2 ) of 0.53 and a root mean square error (RMSE) of 3.28 m (Figure 12).We performed a piecewise averaging of all samples to reveal the relationships between predictor variables and refined CHM (Appendix A).It was found that, with the increase of the refined CHM, the reflectance of all bands showed a downward trend, and the NDVI and EVI showed upward trends.Evergreen forests usually accounted for larger proportions in the groups with higher refined CHMs.We used out-of-bag estimation to evaluate the Random Forest regression model.The canopy height of each sample was predicted, based on about 1/3 of the regression trees.The estimation showed a coefficient of determination (R 2 ) of 0.53 and a root mean square error (RMSE) of 3.28 m (Figure 12).size was about 1 km.Finally, we obtained 357,628 segments with a mean area of 3.17 hectares.

30 m × 30 m Resolution Forest Canopy Height Mapping
As Figure 11 shows, band5 (Near infrared band) and band3 (Green band) were the most important variables in the Random Forest model.Band5 had the highest %IncMSE (18.84%) and band3 had the most significant impact on the IncNodePurity.Band1 (Coastal/Aerosol band) and band4 (Red band) also had relatively high levels of %IncMSE, and IncNodePurity, respectively.
The importance of bands 3 (0.53 μm-0.59μm) and 5 (0.85 μm-0.88 μm) have a bio-physical basis.Band3 is related to chlorophyll.A more intense reflectance in this band usually corresponds to younger forests and therefore, corresponds to a lower canopy height.Band5 is associated with biomass, and biomass is highly correlated with forest canopy height.Therefore, we suggest that they deserve closer attention in future research.If studies use multi-spectral data, other than Landsat, they should make the existence of these two bands a prerequisite.We performed a piecewise averaging of all samples to reveal the relationships between predictor variables and refined CHM (Appendix A).It was found that, with the increase of the refined CHM, the reflectance of all bands showed a downward trend, and the NDVI and EVI showed upward trends.Evergreen forests usually accounted for larger proportions in the groups with higher refined CHMs.We used out-of-bag estimation to evaluate the Random Forest regression model.The canopy height of each sample was predicted, based on about 1/3 of the regression trees.The estimation showed a coefficient of determination (R 2 ) of 0.53 and a root mean square error (RMSE) of 3.28 m (Figure 12). Figure 13 shows the predicted forest canopy height, which was related to the forest type to some extent.For example, two areas, one northeast of the Three Gorges Dam and the other east of the study area, had relatively large canopy heights and both of them were evergreen coniferous forest areas.Figure 13 shows the predicted forest canopy height, which was related to the forest type to some extent.For example, two areas, one northeast of the Three Gorges Dam and the other east of the study area, had relatively large canopy heights and both of them were evergreen coniferous forest areas.The predicted forest canopy height also responded to altitude.In high mountain areas, forest canopy heights gradually decreased as the altitude rose, forming sunken patches centered on the mountain peaks.Some alpine regions in the northwest corner (red areas in Figure 2) of the study area are covered with snow for months every year and the forest canopy heights in these areas were also low.
The study area has long been a region with frequent human activity.Except for some planted forests on the hillsides, orchards and tea trees dominate the areas below the altitude of 800 m.The predicted forest canopy height confirmed this.As Figure 13 shows, except for some planted forests in the east, the forests located below the altitude of 800 m had lower canopy heights than other forests.The spatial distributions of these forests were also fragmented.In contrast, forests in protected areas and scenic spots had higher canopy heights and more continuous distributions.
Our methods produced three kinds of canopy heights: (1) The crude CHM had a resolution of 2 m × 2 m and represented the heights of different positions on the canopy surface.(2) During the refining process, we averaged the crude CHM in 30 m × 30 m cells.Thus, the refined CHM had a resolution of 30 m × 30 m and could be regarded as a mean canopy height.(3) We used the refined CHM as the target value in the Random Forest model, thus the predicted canopy height also represented the mean canopy height.The predicted forest canopy height also responded to altitude.In high mountain areas, forest canopy heights gradually decreased as the altitude rose, forming sunken patches centered on the mountain peaks.Some alpine regions in the northwest corner (red areas in Figure 2) of the study area are covered with snow for months every year and the forest canopy heights in these areas were also low.

Independent Validation of the Predicted Forest Canopy Height
The study area has long been a region with frequent human activity.Except for some planted forests on the hillsides, orchards and tea trees dominate the areas below the altitude of 800 m.The predicted forest canopy height confirmed this.As Figure 13 shows, except for some planted forests in the east, the forests located below the altitude of 800 m had lower canopy heights than other forests.The spatial distributions of these forests were also fragmented.In contrast, forests in protected areas and scenic spots had higher canopy heights and more continuous distributions.
Our methods produced three kinds of canopy heights: (1) The crude CHM had a resolution of 2 m × 2 m and represented the heights of different positions on the canopy surface.(2) During the refining process, we averaged the crude CHM in 30 m × 30 m cells.Thus, the refined CHM had a resolution of 30 m × 30 m and could be regarded as a mean canopy height.(3) We used the refined CHM as the target value in the Random Forest model, thus the predicted canopy height also represented the mean canopy height.

Independent Validation of the Predicted Forest Canopy Height
Figure 14 shows the spatial distribution of the 20 field survey plots on the Landsat image.They are distributed between 390 m and 1520 m above sea level, and their slopes range from 8 • to 42 • .The age of the forests surveyed is between 8 years and 48 years.The measured mean forest heights range from 5.2 m to 15.2 m, with an average of 8.79 m.All the survey plots used for validation are occupied by coniferous and broad-leaved mixed forests.Considering that the forest type classification in this study was rough and the forest type did not play a significant role in the model (Figure 11), the survey plots occupied by other forest types were not used.
Before validation, we averaged the predicted forest canopy height in a 4 × 4 pixel window around the GPS position of each survey plot to reduce the effects of GPS positioning errors and the residuals of topographic correction.The validation result (Figure 15) showed an R 2 of 0.62 and a RMSE of 4.66 m.An offset of 6.5 m was observed.After compensating the offset, the RMSE dropped to 2.64 m.
by coniferous and broad-leaved mixed forests.Considering that the forest type classification in this study was rough and the forest type did not play a significant role in the model (Figure 11), the survey plots occupied by other forest types were not used.
Before validation, we averaged the predicted forest canopy height in a 4 × 4 pixel window around the GPS position of each survey plot to reduce the effects of GPS positioning errors and the residuals of topographic correction.The validation result (Figure 15) showed an R 2 of 0.62 and a RMSE of 4.66 m.An offset of 6.5 m was observed.After compensating the offset, the RMSE dropped to 2.64 m.Two factors might be responsible for the offset in the validation.First, we used 49 GCPs set in vegetation-free areas when extracting the ZY-3 DSM.Most of the GCPs were located on artificial surfaces, such as roads and farmlands in valleys.Since the SRTMGL1 has a spatial resolution of 30 m × 30 m, the extracted elevation of GCPs might be affected by the surrounding hillsides and be higher than the actual elevation.This might shift the ZY-3 DSM upwards, causing an over-estimation of the CHM and, in turn, an over-estimation of the predicted forest canopy height.The second factor, perhaps a more influential factor, might be the topographic correction.An insufficient topographic correction would result in lower reflectance in shady slopes.This would also lead to an overestimation of the predicted forest canopy height, since the reflectance had negative correlation with canopy height, as shown in Figure A1.Further studies may require more data, with higher spatial by coniferous and broad-leaved mixed forests.Considering that the forest type classification in this study was rough and the forest type did not play a significant role in the model (Figure 11), the survey plots occupied by other forest types were not used.
Before validation, we averaged the predicted forest canopy height in a 4 × 4 pixel window around the GPS position of each survey plot to reduce the effects of GPS positioning errors and the residuals of topographic correction.The validation result (Figure 15) showed an R 2 of 0.62 and a RMSE of 4.66 m.An offset of 6.5 m was observed.After compensating the offset, the RMSE dropped to 2.64 m.Two factors might be responsible for the offset in the validation.First, we used 49 GCPs set in vegetation-free areas when extracting the ZY-3 DSM.Most of the GCPs were located on artificial surfaces, such as roads and farmlands in valleys.Since the SRTMGL1 has a spatial resolution of 30 m × 30 m, the extracted elevation of GCPs might be affected by the surrounding hillsides and be higher than the actual elevation.This might shift the ZY-3 DSM upwards, causing an over-estimation of the CHM and, in turn, an over-estimation of the predicted forest canopy height.The second factor, perhaps a more influential factor, might be the topographic correction.An insufficient topographic correction would result in lower reflectance in shady slopes.This would also lead to an overestimation of the predicted forest canopy height, since the reflectance had negative correlation with canopy height, as shown in Figure A1.Further studies may require more data, with higher spatial factors might be responsible for the offset in the validation.First, we used 49 GCPs set in vegetation-free areas when extracting the ZY-3 DSM.Most of the GCPs were located on artificial surfaces, such as roads and farmlands in valleys.Since the SRTMGL1 has a spatial resolution of 30 m × 30 m, the extracted elevation of GCPs might be affected by the surrounding hillsides and be higher than the actual elevation.This might shift the ZY-3 DSM upwards, causing an over-estimation of the CHM and, in turn, an over-estimation of the predicted forest canopy height.The second factor, perhaps a more influential factor, might be the topographic correction.An insufficient topographic correction would result in lower reflectance in shady slopes.This would also lead to an over-estimation of the predicted forest canopy height, since the reflectance had negative correlation with canopy height, as shown in Figure A1.Further studies may require more data, with higher spatial resolution and vertical accuracy.In the future, we plan to apply this method to areas covered by ALS data for a more comprehensive assessment.The ALS data can also be used to calibrate the possible offset of the predicted forest canopy height.

Limitations and Future Outlook
In the refining of the crude CHM, we excluded pixels from areas of moderate or high relief, which is a limitation of this study.However, it was necessary for exploring the correlation between the CHM and predictor variables, since the forest height information was overwhelmed by erroneous values in many places.We believe that to remove this limitation, a more complicated and ingenious method is needed to refine the crude CHM, which is the direction of our future research.In addition, two assumptions guided the extrapolation of the model.First, we assumed that the topographic correction effectively eliminated the effects of terrain and illumination on the reflectance.This meant that the same object should have the same Landsat reflectance, regardless of the slope, aspect, and light source direction.Second, we assumed that the forests with the same type and Landsat reflectance would have the same canopy heights, no matter where they were.However, we had observed many small-scale residuals from the topographic correction in the rugged areas.We also greatly simplified the classification of forest types.These are issues that need to be addressed future research.
The resolution difference between the ZY-3 DSM and SRTMGL1 had a significant impact on the CHM.The ZY-3 DSM had a finer resolution and contained more topographical details than the SRTMGL1.Although we interpolated the SRTMGL1 to 2 m × 2 m resolution before subtracting it from the ZY-3 DSM, small-scale elevation changes on the ZY-3 DSM would still be transferred to the crude CHM and mix with the forest canopy heights.The resolution of the DTM also limited the use of other multi-spectral images with higher resolutions.It was used for the topographic correction of multi-spectral image in mountainous areas, which had a direct impact on the model-predicted forest canopy height.These problems indicated the importance of the DTM.
In the validation of the predicted forest canopy height, only 20 survey plots were used.This was determined by the rough forest type classification and the lack of other validation data, such as LiDAR in the study area.Our next work will attempt to conduct a more detailed forest type classification by using multi-temporal or hyperspectral data.We also plan to apply this method to areas covered by LiDAR data for a more comprehensive assessment.Researchers familiar with LiDAR, especially airborne LiDAR, may not be satisfied with the accuracy of the proposed method.However, we would like to clarify that our goal was to develop a complementary approach to existing techniques, and the complex terrain in the study area was also a big challenge.The method is expected to achieve better accuracy in areas where the composition of forests is simpler and the terrain is gentler.The results in this study could be considered as a conservative assessment of the proposed method.
The use of the ZY-3 and SRTMGL1 data could facilitate the forest canopy height mapping over large regions.The processing steps mentioned in this study were important for solving the problems caused by the resolution differences between data, and for suppressing the interferences from rugged terrain.In this study, the slope difference layer was very important.Multi-resolution segmentation also played a significant role in refining the crude CHM.If we had not used the average slope of segments as a criterion, the R 2 of the out-of-bag estimation would not have exceeded 0.36, and the RMSE would have been larger than 4 m.

Conclusions
This study explored the feasibility of mapping forest canopy height in a mountainous area of Hubei Province, Central China using the China's stereo mapping satellite ZiYuan-3 (ZY-3) data, the Shuttle Radar Topography Mission global 1 arc second data (SRTMGL1), and the Landsat 8 Operational Land Imager (OLI) surface reflectance data.We found that the crude CHM, derived from the ZY-3 and SRTMGL1, contained information on the heights of surface objects, but was contaminated by many erroneous values.The slope difference and the average slope of segments criteria played important roles in refining the crude CHM.The near infrared band and green band of the Landsat 8 data were the most important variables in the Random Forest model.The evaluation of the model using out-of-bag estimation showed an R 2 of 0.53 and a RMSE of 3.28 m.The predicted forest canopy height responded to forest types, altitude, and human activity.We validated the predicted forest canopy height using the mean forest height measured in the field survey plots.The validation showed an R 2 of 0.62 and a RMSE of 2.64 m.It will be easy to combine the proposed method with other techniques.For example, a DSM derived from the X-band SAR images through radargrammetry can be used as a substitute for the ZY-3 DSM in this study.Combining the ZY-3 DSM, with a ALS-derived DTM may improve the quality of the CHM, simplify the refining process and increase the number of training samples.Forest heights derived from the spaceborne LiDAR data, that do not provide full coverage, such as those from the ICESat-2 and GEDI missions, can be used as a substitute for the refined CHM samples in this study to establish and extrapolate the regression model in combination with the Landsat 8 data.

Figure 1 .
Figure 1.Location of the study area.

Figure 1 .
Figure 1.Location of the study area.

Forests 2019, 10 , 21 Figure 3 .
Figure 3. Flow chart of methods.The blue box indicates that the data and operations were based on the ZY-3 coverage.The red box indicates that the data and operations were based on the entire study area.In the flow chart, RPC stands for Rational Polynomial Coefficients files, GCP stands for Ground Control Points, and PCA stands for Principal Component Analysis.

Figure 3 .
Figure 3. Flow chart of methods.The blue box indicates that the data and operations were based on the ZY-3 coverage.The red box indicates that the data and operations were based on the entire study area.In the flow chart, RPC stands for Rational Polynomial Coefficients files, GCP stands for Ground Control Points, and PCA stands for Principal Component Analysis.

Figure 4 .
Figure 4.The GCPs and ZY-3 multi-spectral image.We set 49 GCPs manually in the vegetation-free areas and derived the elevation of the GCPs from the SRTMGL1.

Figure 4 .
Figure 4.The GCPs and ZY-3 multi-spectral image.We set 49 GCPs manually in the vegetation-free areas and derived the elevation of the GCPs from the SRTMGL1.

Figure 4 .
Figure 4.The GCPs and ZY-3 multi-spectral image.We set 49 GCPs manually in the vegetation-free areas and derived the elevation of the GCPs from the SRTMGL1.

Figure 6 .
Figure 6.The ZY-3 nadir view images and the corresponding crude CHMs.The crude CHM's legend is the same as Figure 5's.Subgraphs (a) and (e) show the impact of shadows and steep terrain, with the positive error labelled in red and the negative error labelled in black.(b) and (f), (c) and (g), and (d) and (h) show the vertical information of different surface objects.

Figure 6 .
Figure 6.The ZY-3 nadir view images and the corresponding crude CHMs.The crude CHM's legend is the same as Figure 5's.Subgraphs (a) and (e) show the impact of shadows and steep terrain, with the positive error labelled in red and the negative error labelled in black.(b) and (f), (c) and (g), and (d) and (h) show the vertical information of different surface objects.

ForestsFigure 7 .
Figure 7.The slope difference mask.Subgraph (a) shows the ZY-3 multi-spectral image of the demonstration area.Subgraph (b) shows the crude CHM that was contaminated by terrain fluctuations (the legend is the same as Figure5's).Subgraph (c) shows the slope difference mask.We only preserved the areas with slope differences less than 4.5° (displayed in white).The 4.5° threshold was a compromise option to minimize the interferences from terrain, while retaining enough samples.

Figure 7 .
Figure 7.The slope difference mask.Subgraph (a) shows the ZY-3 multi-spectral image of the demonstration area.Subgraph (b) shows the crude CHM that was contaminated by terrain fluctuations (the legend is the same as Figure 5's).Subgraph (c) shows the slope difference mask.We only preserved the areas with slope differences less than 4.5 • (displayed in white).The 4.5 • threshold was a compromise option to minimize the interferences from terrain, while retaining enough samples.

Forests 2019 ,
10 FOR PEER REVIEW 11 of 21 forestry department for this area.The spatial resolution of the refined CHM samples was 30 m × 30 m.

Figure 8 .
Figure 8. Refined CHM samples on the elevation layer.

Figure 8 .
Figure 8. Refined CHM samples on the elevation layer.

Figure 9 .
Figure 9. Classification result of forest types.

Figure 10
Figure 10 shows a piece of the composited layer of the first principal component (PC1), the second principal component (PC2), and the forest type.The red channel represents PC1, the green channel represents PC2, and the blue channel represents forest type.

Figure 10 .
Figure 10.A piece of the composited layer of the first principal component (PC1), the second principal component (PC2), and the forest type.We performed multi-resolution segmentation, based on the composited layer.

Figure 9 .
Figure 9. Classification result of forest types.

Figure 10
Figure 10 shows a piece of the composited layer of the first principal component (PC1), the second principal component (PC2), and the forest type.The red channel represents PC1, the green channel represents PC2, and the blue channel represents forest type.

Figure 9 .
Figure 9. Classification result of forest types.

Figure 10
Figure 10 shows a piece of the composited layer of the first principal component (PC1), the second principal component (PC2), and the forest type.The red channel represents PC1, the green channel represents PC2, and the blue channel represents forest type.

Figure 10 .
Figure 10.A piece of the composited layer of the first principal component (PC1), the second principal component (PC2), and the forest type.We performed multi-resolution segmentation, based on the composited layer.

Figure 10 .
Figure 10.A piece of the composited layer of the first principal component (PC1), the second principal component (PC2), and the forest type.We performed multi-resolution segmentation, based on the composited layer.

Figure 11 .
Figure 11.Variable importance of the Random Forest model.

Figure 11 .
Figure 11.Variable importance of the Random Forest model.

Figure 11 .
Figure 11.Variable importance of the Random Forest model.

Figure 12 .
Figure 12.Out-of-bag estimation of the Random Forest model.

Figure 12 .
Figure 12.Out-of-bag estimation of the Random Forest model.

Figure 13 .
Figure 13.Forest canopy height predicted by the Random Forest model.

Figure 14
Figure14shows the spatial distribution of the 20 field survey plots on the Landsat image.They are distributed between 390 m and 1520 m above sea level, and their slopes range from 8° to 42°.The age of the forests surveyed is between 8 years and 48 years.The measured mean forest heights range from 5.2 m to 15.2 m, with an average of 8.79 m.All the survey plots used for validation are occupied

Figure 13 .
Figure 13.Forest canopy height predicted by the Random Forest model.

Figure 14 .
Figure 14.Spatial distribution of the survey plots.

Figure 15 .
Figure 15.Independent validation of the predicted forest canopy height.

Figure 14 .
Figure 14.Spatial distribution of the survey plots.

Figure 14 .
Figure 14.Spatial distribution of the survey plots.

Figure 15 .
Figure 15.Independent validation of the predicted forest canopy height.

Figure 15 .
Figure 15.Independent validation of the predicted forest canopy height.