Next Article in Journal
Atypical Pattern of Soil Carbon Stocks along the Slope Position in a Seasonally Dry Tropical Forest in Thailand
Next Article in Special Issue
Do High-Voltage Power Transmission Lines Affect Forest Landscape and Vegetation Growth: Evidence from a Case for Southeastern of China
Previous Article in Journal
Improving Forest Aboveground Biomass (AGB) Estimation by Incorporating Crown Density and Using Landsat 8 OLI Images of a Subtropical Forest in Western Hunan in Central China
Previous Article in Special Issue
Application of Terrestrial Laser Scanner to Evaluate the Influence of Root Collar Geometry on Stump Height after Mechanized Forest Operations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100101, China
2
University of Chinese Academy of Sciences, Beijing 100049, China
3
Academy of Forest Inventory and Planning, State Forestry Administration, Beijing 100714, China
*
Author to whom correspondence should be addressed.
Forests 2019, 10(2), 105; https://doi.org/10.3390/f10020105
Submission received: 11 December 2018 / Revised: 24 January 2019 / Accepted: 28 January 2019 / Published: 29 January 2019

Abstract

:
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.

1. 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 Moderate Resolution Imaging Spectroradiometer (MODIS) data [2,3,14,15,16]. The SAR is well known for its high temporal resolution, because it is less affected by weather and illumination conditions [17]. The SAR-based techniques for forest canopy height mapping include radargrammetry, interferometry SAR (InSAR), and polarimetric interferometric SAR (PolInSAR). The radargrammetry is based on SAR stereo images [5,18,19]. The InSAR is based on the phase differences between two complex SAR images [6,20]. The PolInSAR includes phase difference methods and the model-based methods, which rely on the complex coherent coefficient model [7,21]. In these methods, the TerraSAR-X (X-band), Radarsat-2 (C-band), and ALOS-2 (L-band) are commonly used SAR data sources. It is expected that the techniques for mapping forest canopy height using the LiDAR and SAR will be further developed when data from the ICESat-2 [22], the GEDI [23], the Biomass (P-band) (due for launch in 2021) [24], and the Tandem-L (L-band) (due for launch in 2022) [25] missions become available. The airborne DP systems have a large field of view, high cruising altitude, fast data acquisition speed, and easy flight planning [26]. The spaceborne DP systems can perform continuous and repetitive observations. Several spaceborne systems, such as the WorldView series, provide sub-meter resolution stereo images from which high-quality digital surface models (DSM) can be extracted [27]. Since the DSM represents the surface of the earth, including vegetation, buildings and other man-made features, extracting canopy height requires a digital terrain model (DTM) that provides the bare earth elevation (another term that needs to be distinguished here is the digital elevation model or DEM, which is a generic term that can refer to either a DSM or DTM). In recent years, some studies have mapped forest canopy height using airborne or spaceborne ultra-high spatial resolution (finer than 1 m × 1 m) stereo images and the ALS-derived DTM data [9,28,29,30,31,32].
Existing techniques provide a variety of means for mapping forest height. The ALS is probably the most promising and mature technique. However, the demands for forest canopy height mapping in practical applications are often not met. Because the data needed in the existing techniques (e.g., LiDAR or SAR data) may not be available when it comes to a given study area. The lack of data demonstrates the necessity of making full use of the free datasets, exploiting the latest, cheaper datasets, and combining them with models for forest height mapping. In this study, we used the ZiYuan-3 (ZY-3) stereo images to derive the DSM and the Shuttle Radar Topography Mission global 1 arc second data (SRTMGL1) as the DTM. We correlated the extracted canopy height model (CHM) with the Landsat data through a regression model and obtained a wall-to-wall forest canopy height product after extrapolation.

2. 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.70 km2. 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].

3. Data

3.1. 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 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.

3.2. 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.

3.3. 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.

3.4. 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.

4. 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.

4.1. 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 −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.
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].

4.2. 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)
ρλ,t = bλ + mλ × cos(i)
cλ = bλ / mλ
ρλ,n = ρλ,t × (cos(sz) + cλ) / (cos(i) + cλ)
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.

4.3. 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.

5. Results and Discussion

5.1. 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.

5.2. 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.02 km2 and had an average distribution-altitude of 923 m. The evergreen forests showed a clustered distribution. They covered a total area of 3,240.54 km2 and had an average distribution-altitude of 917 m.

5.3. PCA and Multi-resolution Segmentation

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.
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.

5.4. 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 (R2) 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.
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.

5.5. 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 R2 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 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.

5.6. 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 in 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 R2 of the out-of-bag estimation would not have exceeded 0.36, and the RMSE would have been larger than 4 m.

6. 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 R2 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 R2 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.

Author Contributions

Data curation, M.L., Y.D., and X.N.; Funding acquisition, C.C.; Investigation, M.L.; Methodology, M.L.; Project administration, C.C., Y.D., and X.N.; Resources, C.C. and Y.D.; Software, Y.D.; Supervision, C.C.; Writing-original draft, M.L.; Writing-review & editing, C.C., and X.N.

Funding

This research was funded by the National Key Research and Development Program of China grant number 2016YFB0501505.

Acknowledgments

Landsat surface reflectance product courtesy of the U.S. Geological Survey. SRTMGL1 Version 3 product was retrieved from the online Data Pool, courtesy of the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

We sorted all of the samples, based on their refined CHM values and segmented them at 0.5 m intervals. We then averaged the refined CHM and predictor variables in each interval (i.e., piecewise averaging) and created scatter plots to demonstrate the relationship between them. For the forest type, we counted the number of samples instead of averaging. Figure A1 shows the scatter plots of the piecewise averaged predictor variables and refined CHM. Data points with refined CHM less than 3 m or larger than 23 m were discrete, because these intervals contained few samples.
Figure A1. Scatter plots of the piecewise averaged predictor variables and refined CHM.
Figure A1. Scatter plots of the piecewise averaged predictor variables and refined CHM.
Forests 10 00105 g0a1

References

  1. Lefsky, M.A.; Keller, M.; Pang, Y.; de Camargo, P.B.; Hunter, M.O. Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveforms. J. Appl. Remote Sens. 2007, 1, 1–18. [Google Scholar]
  2. Lefsky, M.A. A global forest canopy height map from the Moderate Resolution Imaging Spectroradiometer and the Geoscience Laser Altimeter System. Geophys. Res. Lett. 2010, 37, 1–5. [Google Scholar] [CrossRef]
  3. Simard, M.; Pinto, N.; Fisher, J.B.; Baccini, A. Mapping forest canopy height globally with spaceborne lidar. J. Geophys. Res. Biogeosci. 2011, 116, 1–12. [Google Scholar] [CrossRef]
  4. Roussel, J.R.; Caspersen, J.; Béland, M.; Thomas, S.; Achim, A. Removing bias from LiDAR-based estimates of canopy height: Accounting for the effects of pulse density and footprint size. Remote Sens. Environ. 2017, 198, 1–16. [Google Scholar] [CrossRef]
  5. Vastaranta, M.; Holopainen, M.; Karjalainen, M.; Kankare, V.; Hyyppa, J.; Kaasalainen, S.; Hyyppa, H. SAR radargrammetry and scanning LiDAR in predicting forest canopy height. IEEE Int. Geosci. Remote Sens. Symp. 2012, 53, 6515–6518. [Google Scholar]
  6. Balzter, H.; Rowland, C.S.; Saich, P. Forest canopy height and carbon estimation at Monks Wood National Nature Reserve, UK, using dual-wavelength SAR interferometry. Remote Sens. Environ. 2007, 108, 224–239. [Google Scholar] [CrossRef]
  7. Ghasemi, N.; Tolpekin, V.; Stein, A. A modified model for estimating tree height from PolInSAR with compensation for temporal decorrelation. Int. J. Appl. Earth Obs. 2018, 73, 313–322. [Google Scholar] [CrossRef]
  8. Balenovic, I.; Seletkovic, A.; Pernar, R.; Jazbec, A. Estimation of the mean tree height of forest stands by photogrammetric measurement using digital aerial images of high spatial resolution. Ann. For. Res. 2015, 58, 125–143. [Google Scholar] [CrossRef]
  9. Herrero, H.M.; Felipe, G.B.; Belmar, L.S.; Hernandez, L.D.; Rodriguez, G.P.; Gonzalez, A.D. Dense Canopy Height Model from a low-cost photogrammetric platform and LiDAR data. Trees-Struct. Funct. 2016, 30, 1287–1301. [Google Scholar] [CrossRef]
  10. Zhou, T.; Popescu, S.C.; Krause, K.; Sheridan, R.D.; Putman, E. Gold—A novel deconvolution algorithm with optimization for waveform LiDAR processing. ISPRS-J. Photogramm. Remote Sens. 2017, 129, 131–150. [Google Scholar] [CrossRef]
  11. Anderson, J.; Martin, M.E.; Smith, M.L.; Dubayah, R.O.; Hofton, M.A.; Hyde, P.; Peterson, B.E.; Blair, J.B.; Knox, R.G. The use of waveform lidar to measure northern temperate mixed conifer and deciduous forest structure in New Hampshire. Remote Sens. Environ. 2006, 105, 248–261. [Google Scholar] [CrossRef]
  12. Boudreau, J.; Nelson, R.F.; Margolis, H.A.; Beaudoin, A.; Guindon, L.; Kimes, D.S. Regional aboveground forest biomass using airborne and spaceborne LiDAR in Québec. Remote Sens. Environ. 2008, 112, 3876–3890. [Google Scholar] [CrossRef]
  13. Huang, H.; Liu, C.; Wang, X.; Biging, G.S.; Chen, Y.; Yang, J.; Gong, P. Mapping vegetation heights in China using slope correction ICESat data, SRTM, MODIS-derived and climate data. ISPRS-J. Photogramm. Remote Sens. 2017, 129, 189–199. [Google Scholar] [CrossRef]
  14. Bolton, D.K.; Coops, N.C.; Wulder, M.A. Investigating the agreement between global canopy height maps and airborne LiDAR derived height estimates over Canada. Can. J. Remote Sens. 2013, 39, S139–S151. [Google Scholar] [CrossRef]
  15. Hilker, T.; Wulder, M.A.; Coops, N.C. Update of forest inventory data with LiDAR and high spatial resolution satellite imagery. Can. J. Remote Sens. 2008, 34, 5–12. [Google Scholar] [CrossRef]
  16. Hansen, M.C.; Potapov, P.V.; Goetz, S.J.; Turubanova, S.; Tyukavina, A.; Krylov, A.; Kommareddy, A.; Egorov, A. Mapping tree height distributions in Sub-Saharan Africa using Landsat 7 and 8 data. Remote Sens. Environ. 2016, 185, 221–232. [Google Scholar] [CrossRef] [Green Version]
  17. Garestier, F.; Dubois-Fernandez, P.C.; Guyon, D.; Le Toan, T. Forest biophysical parameter estimation using L- and P-band polarimetric SAR data. IEEE Trans. Geosci. Remote Sens. 2009, 47, 3379–3388. [Google Scholar] [CrossRef]
  18. Karjalainen, M.; Kankare, V.; Vastaranta, M.; Holopainen, M.; Hyyppä, J. Prediction of plot-level forest variables using TerraSAR-X stereo SAR data. Remote Sens. Environ. 2012, 117, 338–347. [Google Scholar] [CrossRef]
  19. Capaldo, P.; Nascetti, A.; Porfiri, M.; Pieralice, F.; Fratarcangeli, F.; Crespi, M.; Toutin, T. Evaluation and comparison of different radargrammetric approaches for Digital Surface Models generation from COSMO-SkyMed, TerraSAR-X, RADARSAT-2 imagery: Analysis of Beauport (Canada) test site. ISPRS-J. Photogramm. Remote Sens. 2015, 100, 60–70. [Google Scholar] [CrossRef]
  20. Solberg, S.; Astrup, R.; Breidenbach, J.; Nilsen, B.; Weydahl, D. Monitoring spruce volume and biomass with InSAR data from TanDEM-X. Remote Sens. Environ. 2013, 139, 60–67. [Google Scholar] [CrossRef]
  21. Zhang, Y.; He, C.; Xu, X.; Chen, D. Forest vertical parameter estimation using PolInSAR imagery based on radiometric correction. ISPRS Int. J. Geo-Inf. 2016, 5, 186. [Google Scholar] [CrossRef]
  22. Neuenschwander, A.; Pitts, K. The ATL08 land and vegetation product for the ICESat-2 Mission. Remote Sens. Environ. 2019, 221, 247–259. [Google Scholar] [CrossRef]
  23. Qi, W.; Dubayah, R.O. Combining Tandem-X InSAR and simulated GEDI lidar observations for forest structure mapping. Remote Sens. Environ. 2016, 187, 253–266. [Google Scholar] [CrossRef]
  24. Le Toan, T.; Quegan, S.; Davidson, M.W.J.; Balzter, H.; Paillou, P.; Papathanassiou, K.; Plummer, S.; Rocca, F.; Saatchi, S.; Shugart, H.; et al. The BIOMASS mission: Mapping global forest biomass to better understand the terrestrial carbon cycle. Remote Sens. Environ. 2011, 115, 2850–2860. [Google Scholar] [CrossRef] [Green Version]
  25. Moreira, A.; Krieger, G.; Hajnsek, I.; Papathanassiou, K.; Younis, M.; Lopez-Dekker, P.; Huber, S.; Villano, M.; Pardini, M.; Eineder, M.; et al. Tandem-L: A Highly Innovative Bistatic SAR Mission for Global Observation of Dynamic Processes on the Earth’s Surface. IEEE Geosci. Remote Sens. Mag. 2015, 3, 8–23. [Google Scholar] [CrossRef]
  26. White, J.C.; Wulder, M.A.; Vastaranta, M.; Coops, N.C.; Pitt, D.; Woods, M. The utility of image-based point clouds for forest inventory: A comparison with airborne laser scanning. Forests 2013, 4, 518–536. [Google Scholar] [CrossRef]
  27. Aguilar, M.A.; Del Mar Saldaña, M.; Aguilar, F.J. Generation and quality assessment of stereo-extracted DSM from GeoEye-1 and WorldView-2 imagery. IEEE Trans. Geosci. Remote Sens. 2014, 52, 1259–1271. [Google Scholar] [CrossRef]
  28. Immitzer, M.; Stepper, C.; Böck, S.; Straub, C.; Atzberger, C. Use of WorldView-2 stereo imagery and National Forest Inventory data for wall-to-wall mapping of growing stock. For. Ecol. Manag. 2016, 359, 232–246. [Google Scholar] [CrossRef]
  29. Jayathunga, S.; Owari, T.; Tsuyuki, S. Evaluating the Performance of Photogrammetric Products Using Fixed-Wing UAV Imagery over a Mixed Conifer-Broadleaf Forest: Comparison with Airborne Laser Scanning. Remote Sens. 2018, 10, 187. [Google Scholar] [CrossRef]
  30. Véga, C.; St-Onge, B. Height growth reconstruction of a boreal forest canopy over a period of 58 years using a combination of photogrammetric and lidar models. Remote Sens. Environ. 2008, 112, 1784–1794. [Google Scholar] [CrossRef] [Green Version]
  31. Järnstedt, J.; Pekkarinen, A.; Tuominen, S.; Ginzler, C.; Holopainen, M.; Viitala, R. Forest variable estimation using a high-resolution digital surface model. ISPRS-J. Photogramm. Remote Sens. 2012, 74, 78–84. [Google Scholar] [CrossRef]
  32. Pearse, G.D.; Dash, J.P.; Persson, H.J.; Watt, M.S. Comparison of high-density LiDAR and satellite photogrammetry for forest inventory. ISPRS-J. Photogramm. Remote Sens. 2018, 142, 257–267. [Google Scholar] [CrossRef]
  33. Gan, J.; Ge, J.; Liu, Y.; Wang, Z.; Wu, X.; Ding, S. Forest ecosystem quality change over ten years (2000–2010) in the Three Gorges Dalaoling Nature Reserve. Plant Sci. J. 2015, 33, 766–774. [Google Scholar]
  34. Wang, T.; Zhang, G.; Li, D.; Tang, X.; Jiang, Y.; Pan, H.; Zhu, X.; Fang, C. Geometric accuracy validation for ZY-3 satellite imagery. IEEE Geosci. Remote Sens. Lett. 2014, 11, 1168–1171. [Google Scholar] [CrossRef]
  35. Li, D. China’s first civilian three-line-array stereo mapping satellite ZY-3. Acta Geodaetica et Cartographica Sinica 2012, 41, 317–322. [Google Scholar]
  36. Pan, H.; Zhang, G.; Tang, X.; Wang, X.; Zhou, P.; Xu, M.; Li, D. Accuracy analysis and verification of ZY-3 products. Acta Geodaetica et Cartographica Sinica 2013, 42, 738–744. [Google Scholar]
  37. Ni, W.; Sun, G.; Ranson, K.J.; Pang, Y.; Zhang, Z.; Yao, W. Extraction of ground surface elevation from ZY-3 winter stereo imagery over deciduous forested areas. Remote Sens. Environ. 2015, 159, 194–202. [Google Scholar] [CrossRef]
  38. Farr, T.G.; Kobrick, M. Shuttle radar topography mission produces a wealth of data. Eos 2000, 81, 583–585. [Google Scholar] [CrossRef]
  39. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The Shuttle Radar Topography Mission. Rev. Geophys. 2007, 45, 1–43. [Google Scholar] [CrossRef]
  40. Kobrick, M. On the toes of giants-How SRTM was born. Photogramm. Eng. Remote Sens. 2006, 72, 206–210. [Google Scholar]
  41. Nikolakopoulos, K.G.; Kamaratakis, E.K.; Chrysoulakis, N. SRTM vs ASTER elevation products. Comparison for two regions in Crete, Greece. Int. J. Remote Sens. 2006, 27, 4819–4838. [Google Scholar] [CrossRef]
  42. Rosen, P.A.; Hensley, S.; Joughin, I.R.; Li, F.K.; Madsen, S.N.; Rodríguez, E.; Goldstein, R.M. Synthetic aperture radar interferometry. Proc. IEEE 2000, 88, 333–380. [Google Scholar] [CrossRef]
  43. Hu, Z.; Peng, J.; Hou, Y.; Shan, J. Evaluation of recently released open global digital elevation models of Hubei, China. Remote Sens. 2017, 9, 262. [Google Scholar] [CrossRef]
  44. Krieger, G.; Moreira, A.; Fiedler, H.; Hajnsek, I.; Werner, M.; Younis, M.; Zink, M. TanDEM-X: A satellite formation for high-resolution SAR interferometry. IEEE Trans. Geosci. Remote Sens. 2007, 45, 3317–3341. [Google Scholar] [CrossRef]
  45. De Oliveira, C.G.; Paradella, W.R.; da Silva, A.D.Q. Assessment of radargrammetric DSMs from TerraSAR-X Stripmap images in a mountainous relief area of the Amazon region. ISPRS-J. Photogramm. Remote Sens. 2011, 66, 67–72. [Google Scholar] [CrossRef]
  46. Su, Y.; Guo, Q. A practical method for SRTM DEM correction over vegetated mountain areas. ISPRS-J. Photogramm. Remote Sens. 2014, 87, 216–228. [Google Scholar] [CrossRef]
  47. Vermote, E.; Justice, C.; Claverie, M.; Franch, B. Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product. Remote Sens. Environ. 2016, 185, 46–56. [Google Scholar] [CrossRef]
  48. Zhu, Z.; Woodcock, C.E. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sens. Environ. 2012, 118, 83–94. [Google Scholar] [CrossRef]
  49. Wang, Y.; Li, G.; Ding, J.; Guo, Z.; Tang, S.; Wang, C.; Huang, Q.; Liu, R.; Chen, J.M. A combined GLAS and MODIS estimation of the global distribution of mean forest canopy height. Remote Sens. Environ. 2016, 174, 24–43. [Google Scholar] [CrossRef] [Green Version]
  50. Teillet, P.M.; Guindon, B.; Goodenough, D.G. On the slope-aspect correction of multispectral scanner data. Can. J. Remote Sens. 1982, 8, 84–106. [Google Scholar] [CrossRef]
  51. Breiman, L. Random forest. Mach. Learn. 2001, 45, 5–32. [Google Scholar] [CrossRef]
  52. Congalton, R.G. A review of assessing the accuracy of classifications of remotely sensed data. Remote Sens. Environ. 1991, 37, 35–46. [Google Scholar] [CrossRef]
Figure 1. Location of the study area.
Figure 1. Location of the study area.
Forests 10 00105 g001
Figure 2. Elevation (SRTMGL1) of the study area.
Figure 2. Elevation (SRTMGL1) of the study area.
Forests 10 00105 g002
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. 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.
Forests 10 00105 g003
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. 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.
Forests 10 00105 g004
Figure 5. The crude canopy height model (CHM) extracted from the ZY-3 DSM and SRTMGL1.
Figure 5. The crude canopy height model (CHM) extracted from the ZY-3 DSM and SRTMGL1.
Forests 10 00105 g005
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. 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.
Forests 10 00105 g006
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.
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 10 00105 g007
Figure 8. Refined CHM samples on the elevation layer.
Figure 8. Refined CHM samples on the elevation layer.
Forests 10 00105 g008
Figure 9. Classification result of forest types.
Figure 9. Classification result of forest types.
Forests 10 00105 g009
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. 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.
Forests 10 00105 g010
Figure 11. Variable importance of the Random Forest model.
Figure 11. Variable importance of the Random Forest model.
Forests 10 00105 g011
Figure 12. Out-of-bag estimation of the Random Forest model.
Figure 12. Out-of-bag estimation of the Random Forest model.
Forests 10 00105 g012
Figure 13. Forest canopy height predicted by the Random Forest model.
Figure 13. Forest canopy height predicted by the Random Forest model.
Forests 10 00105 g013
Figure 14. Spatial distribution of the survey plots.
Figure 14. Spatial distribution of the survey plots.
Forests 10 00105 g014
Figure 15. Independent validation of the predicted forest canopy height.
Figure 15. Independent validation of the predicted forest canopy height.
Forests 10 00105 g015

Share and Cite

MDPI and ACS Style

Liu, M.; Cao, C.; Dang, Y.; Ni, X. Mapping Forest Canopy Height in Mountainous Areas Using ZiYuan-3 Stereo Images and Landsat Data. Forests 2019, 10, 105. https://doi.org/10.3390/f10020105

AMA Style

Liu M, Cao C, Dang Y, Ni X. Mapping Forest Canopy Height in Mountainous Areas Using ZiYuan-3 Stereo Images and Landsat Data. Forests. 2019; 10(2):105. https://doi.org/10.3390/f10020105

Chicago/Turabian Style

Liu, Mingbo, Chunxiang Cao, Yongfeng Dang, and Xiliang Ni. 2019. "Mapping Forest Canopy Height in Mountainous Areas Using ZiYuan-3 Stereo Images and Landsat Data" Forests 10, no. 2: 105. https://doi.org/10.3390/f10020105

APA Style

Liu, M., Cao, C., Dang, Y., & Ni, X. (2019). Mapping Forest Canopy Height in Mountainous Areas Using ZiYuan-3 Stereo Images and Landsat Data. Forests, 10(2), 105. https://doi.org/10.3390/f10020105

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop