Mapping Canopy Heights of Poplar Plantations in Plain Areas Using ZY 302 Stereo and Multispectral Data

Forest canopy height plays an important role in forest management and ecosystem modeling. There are a variety of techniques employed to map forest height using remote sensing data but it is still necessary to explore the use of new data and methods. In this study, we demonstrate an approach for mapping canopy heights of poplar plantations in plain areas through a combination of stereo and multispectral data from China’s latest civilian stereo mapping satellite ZY3-02. First, a digital surface model (DSM) was extracted using photogrammetry methods. Then, canopy samples and ground samples were selected through manual interpretation. Canopy height samples were obtained by calculating the DSM elevation differences between the canopy samples and ground samples. A regression model was used to correlate the reflectance of a ZY3-02 multispectral image with the canopy height samples, in which the red band and green band reflectance were selected as predictors. Finally, the model was extrapolated to the entire study area and a wall-to-wall forest canopy height map was obtained. The validation of the predicted canopy height map reported a coefficient of determination (R2) of 0.72 and a root mean square error (RMSE) of 1.58 m. This study demonstrates the capacity of ZY3-02 data for mapping the canopy height of pure plantations in plain areas.


Introduction
Forest canopy height is a critical parameter for forest management and ecosystem modeling.Field measurements are essential but labor intensive and costly, especially in extensive and remote areas [1,2].A variety of methods for mapping forest height using remote sensing data have been developed in recent years: the light detection and ranging (LiDAR) technique uses point cloud or waveform data from laser pulses to detect the vertical structure of a forest.Distance is determined by measuring the time taken for the laser pulse to travel between the sensor and target.In recent years, the airborne laser scanning (ALS) technique has made rapid progress and is now widely used in forestry applications.It has proven to be a reliable LiDAR method which can measure forest height with sub-meter accuracy [3,4].The Geoscience Laser Altimeter System (GLAS) onboard NASA's Ice, Cloud, and land Elevation Satellite (ICESat) collected data from 2003 to 2009.It offered an unprecedented opportunity for estimating forest height at a global scale [5].Lefsky and Simard et al. have produced global forest height maps by combining the GLAS data with Moderate Resolution Imaging Spectroradiometer (MODIS) data [6,7].Recent developments include new terrain correction methods [8] and the combination of GLAS data with Landsat time-series data [9].It is expected that the application of spaceborne LiDAR in large-scale forest height mapping will be further developed when data from the ICESat-2 [10] and GEDI (Global Ecosystem Dynamics Investigation, a LiDAR system mounted on the International Space Station) [11] become available.
The synthetic aperture radar (SAR) techniques used for forest height mapping include radargrammetry, interferometry SAR (InSAR), and polarimetric interferometric SAR (PolInSAR).A major advantage of SAR is its high temporal resolution, because SAR has the ability to acquire data under various illumination and weather conditions.Radargrammetry is based on SAR stereo images.It has recently received attention because of the emergence of high spatial resolution SAR images [12].Radargrammetry needs a digital terrain model (DTM) for deriving forest height [13].InSAR is based on the phase differences between two complex SAR images, which also need a DTM for deriving forest height.In addition, the time interval between the two images should be short to avoid decoherence [14,15].PolInSAR includes phase difference methods and model-based methods [16,17].Based on the coherence optimization method, Cloude and Papathanassiou were the first to use PolInSAR in forest height mapping [18].Yamada et al. have used the ESPRIT method to separate phase centers of the canopy and ground [19].The random volume over ground (RVoG) model is the most commonly used model among the model-based methods.Papathanassiou et al. have used the six-dimension non-liner iteration method for the inversion of the RVoG model [20].Cloude and Papathanassiou have proposed the three-stage inversion method, which has improved the accuracy and reduced the complexity of the inversion process [21].The most commonly used SAR data include the TerraSAR-X (X-band), Radarsat-2 (C-band), and ALOS-2 (L-band).The data from the Biomass (P-band) (due for launch in 2021) and Tandem-L (L-band) (due for launch in 2022) missions will alleviate the shortage of L-band and P-band SAR data suitable for the extraction of forest structural parameters [22,23].
Since the emergence of advanced sensors has improved the radiometric and spatial resolution, and advances in computing technology have made complex algorithms for image matching practical, the application of digital photogrammetry (DP) in forestry has received more and more attention in recent years [24,25].DP uses a rigorous or rational function model to characterize the geometry of the acquisition system and uses parallax differences to compute the coordinates and elevations of the matched pixels.Airborne DP systems and several spaceborne DP systems, such as the WorldView series, provide sub-meter resolution stereo images from which a high-quality digital surface model (DSM) can be derived [26].The forest canopy height can be obtained by calculating the difference between the DSM and the ALS-derived or field-surveyed DTM that provides the bare earth elevation [27][28][29].Because the cost of space-borne DP data is usually less, it is ideal for repetitive forest canopy height surveys over a wide area.In addition, due to the long history of the photogrammetry technique, there are a considerable number of photographs in several countries.By digitizing these photographs, researchers may be able to study changes in forest height over a long period of time [29].
However, compared with LiDAR-based techniques, few studies have been carried out on mapping forest canopy heights using surging stereo images.This is partly due to the lack of DTMs.One of the objectives of this study was to evaluate an alternative method for obtaining canopy height samples.In this study, we extracted a DSM from the stereo images of the ZY3-02 satellite.Canopy height samples of poplar plantations were derived from the DSM based on manually selected canopy and ground samples.In order to obtain a canopy height map with complete horizontal coverage, we combined canopy height samples with multispectral image.There is theoretical and empirical evidence indicating that both biochemical properties (e.g.water and pigment content) and the structures of forests are the driving forces regarding the response of multispectral reflectance to canopy height [30,31].For example, higher forests usually have higher chlorophyll content per unit area, resulting in stronger reflections in the green band.At the same time, higher forests usually have more complex structures that create more shadows, causing less light to be reflected into the sensor.The dominant effect depends on the specific spectral band and forest type.Many researchers have utilized the correlation between the reflectance of multispectral images and forest height in their studies [6,9,32].In this study, we correlated the reflectance of a ZY3-02 multispectral image with canopy height samples through a regression model and produced a wall-to-wall forest canopy height map of the entire study area after extrapolation.Hence, this study is also an evaluation of the potential of ZY3-02 data in forest canopy height mapping.In addition, we propose methods for setting ground control points (GCPs) and extracting the GCPs' elevations from free terrain data.

Study Area
Our study area was located on the southern edge of Horqin Sandy Land, eastern Inner Mongolia, China, with a total area of 2500 km 2 (42 • 17'11"-42 • 48'44"N, 119 • 34'6"-120 • 18'15"E) (Figure 1).The area belongs to the mid-temperate sub-arid zone, with an annual precipitation between 310 and 460 mm.The altitude there ranges from 396 to 855 m.The terrain is relatively flat with an average slope of 3.4 • .Most of the forests in the region are poplar plantations, aged between five and 25 years.These plantations are grid-and patch-shaped.They were planted for the Three-North shelter forest program, a large-scale afforestation project in China which began in 1978, and play an important role in blocking sandstorms, preventing soil erosion, and improving the ecological environment.
ISPRS Int.J. Geo-Inf.2019, 8 FOR PEER REVIEW 3 canopy height samples through a regression model and produced a wall-to-wall forest canopy height map of the entire study area after extrapolation.Hence, this study is also an evaluation of the potential of ZY3-02 data in forest canopy height mapping.In addition, we propose methods for setting ground control points (GCPs) and extracting the GCPs' elevations from free terrain data.

Study Area
Our study area was located on the southern edge of Horqin Sandy Land, eastern Inner Mongolia, China, with a total area of 2500 km 2 (42°17'11"-42°48'44"N, 119°34'6"-120°18'15"E) (Figure 1).The area belongs to the mid-temperate sub-arid zone, with an annual precipitation between 310 and 460 mm.The altitude there ranges from 396 to 855 m.The terrain is relatively flat with an average slope of 3.4°.Most of the forests in the region are poplar plantations, aged between five and 25 years.These plantations are grid-and patch-shaped.They were planted for the Three-North shelter forest program, a large-scale afforestation project in China which began in 1978, and play an important role in blocking sandstorms, preventing soil erosion, and improving the ecological environment.

ZY3-02 Data
The ZY3-02 satellite was launched on May 30, 2016, with a design life of five years.It is an upgraded successor to China's first civilian stereo mapping satellite, ZY-3, which was launched in 2012.The ZY3-02 satellite operates in a sun-synchronous near-polar orbit at an altitude of 505 km.It carries a multispectral camera and three panchromatic cameras pointing forward, nadir, and backward.The multispectral camera consists of four channels: blue (450-520 nm), green (520-590 nm), red (630-690 nm), and near infrared (770-890 nm), with a ground sample distance (GSD) of 5.8 m.The forward and backward cameras are arranged at an inclination of ±22° from nadir to realize a base-to-height (B/H) ratio of 0.88.The GSD of the nadir panchromatic camera is 2.1 m.The GSD of the forward and backward panchromatic cameras have been improved from 3.5 m in ZY-3 to 2.7 m. ZY3-02 images have a radiometric resolution of 10 bits and a swath width of 50 km.According to an assessment made by Xu et al. [33], the planimetric and vertical accuracies of the ZY3-02 sensorcorrected products are better than 2.5 m and 2 m, respectively, with a few GCPs.
The ZY3-02 sensor-corrected products used in this study were acquired on September 18, 2017 and were provided by the China Centre for Resources Satellite Data and Application (CRESDA,

ZY3-02 Data
The ZY3-02 satellite was launched on May 30, 2016, with a design life of five years.It is an upgraded successor to China's first civilian stereo mapping satellite, ZY-3, which was launched in 2012.The ZY3-02 satellite operates in a sun-synchronous near-polar orbit at an altitude of 505 km.It carries a multispectral camera and three panchromatic cameras pointing forward, nadir, and backward.The multispectral camera consists of four channels: blue (450-520 nm), green (520-590 nm), red (630-690 nm), and near infrared (770-890 nm), with a ground sample distance (GSD) of 5.8 m.The forward and backward cameras are arranged at an inclination of ±22 • from nadir to realize a base-to-height (B/H) ratio of 0.88.The GSD of the nadir panchromatic camera is 2.1 m.The GSD of the forward and backward panchromatic cameras have been improved from 3.5 m in ZY-3 to 2.7 m. ZY3-02 images have a radiometric resolution of 10 bits and a swath width of 50 km.According to an assessment made by Xu et al. [33], the planimetric and vertical accuracies of the ZY3-02 sensor-corrected products are better than 2.5 m and 2 m, respectively, with a few GCPs.
The ZY3-02 sensor-corrected products used in this study were acquired on September 18 2017 and were provided by the China Centre for Resources Satellite Data and Application (CRESDA, http://www.cresda.com/CN/).We used the forward/backward stereo images to extract the DSM.The multispectral image was used to classify the landcover and predict forest canopy heights.Both the stereo and multispectral images were able to completely cover the study area.At the time of data acquisition, there was no cloud over the study area.

SRTMGL1 Data
The Shuttle Radar Topography Mission (SRTM) collected data with a C-band radar interferometry system onboard the space shuttle Endeavour from February 11 to February 22, 2000.The SRTM product provides the elevation of the land between 60 • north and 56 • south latitude (covering more than 80% of Earth's total landmass) [34].The SRTM product was first released in 2003 and has been updated several times since.The SRTM global 1 arc second data (SRTMGL1) for China have been available since July 2015 and have a spatial resolution of about 27 m for our study area.The vertical accuracy of SRTMGL1 in plain areas is about 1.9 m according to an assessment made by Hu et al. [35].We downloaded the void-filled SRTMGL1 Version 3 product from NASA's website (https://earthdata.nasa.gov/).This data was used to determine the elevation of GCPs.

Methods
In this study we first extracted a DSM from the ZY3-02 forward/backward stereo images using photogrammetry methods.Then, the canopy samples and ground samples were selected through the manual interpretation of a ZY3-02 multispectral image.The canopy height samples were obtained by calculating the DSM elevation difference for each pair of canopy/ground samples.After that, the canopy height samples were divided into a training set and validation set.A multiple linear regression model was established to correlate the reflectance of the multispectral image with the canopy heights of the training samples.Finally, the model was extrapolated to all forest regions in the study area and a wall-to-wall forest canopy height map was obtained.The map was validated independently by the validation samples (Figure 2).http://www.cresda.com/CN/).We used the forward/backward stereo images to extract the DSM.The multispectral image was used to classify the landcover and predict forest canopy heights.Both the stereo and multispectral images were able to completely cover the study area.At the time of data acquisition, there was no cloud over the study area.

SRTMGL1 Data
The Shuttle Radar Topography Mission (SRTM) collected data with a C-band radar interferometry system onboard the space shuttle Endeavour from February 11 to February 22, 2000.
The SRTM product provides the elevation of the land between 60° north and 56° south latitude (covering more than 80% of Earth's total landmass) [34].The SRTM product was first released in 2003 and has been updated several times since.The SRTM global 1 arc second data (SRTMGL1) for China have been available since July 2015 and have a spatial resolution of about 27 m for our study area.The vertical accuracy of SRTMGL1 in plain areas is about 1.9 m according to an assessment made by Hu et al. [35].We downloaded the void-filled SRTMGL1 Version 3 product from NASA's website (https://earthdata.nasa.gov/).This data was used to determine the elevation of GCPs.

Methods
In this study we first extracted a DSM from the ZY3-02 forward/backward stereo images using photogrammetry methods.Then, the canopy samples and ground samples were selected through the manual interpretation of a ZY3-02 multispectral image.The canopy height samples were obtained by calculating the DSM elevation difference for each pair of canopy/ground samples.After that, the canopy height samples were divided into a training set and validation set.A multiple linear regression model was established to correlate the reflectance of the multispectral image with the canopy heights of the training samples.Finally, the model was extrapolated to all forest regions in the study area and a wall-to-wall forest canopy height map was obtained.The map was validated independently by the validation samples (Figure 2).

DSM Extraction
We extracted DSM using the OrthoEngine module of the PCI Geomatica software (PCI Geomatics Enterprises, Inc., Canada).The module uses the polynomial coefficients, GCPs, and tie

DSM Extraction
We extracted DSM using the OrthoEngine module of the PCI Geomatica software (PCI Geomatics Enterprises, Inc., Canada).The module uses the polynomial coefficients, GCPs, and tie points (TIEs) to compute a math model that relates the rows and columns of the matched pixels with ground coordinates and elevations.The polynomial coefficients were provided in the rational polynomial coefficients (RPC) files distributed with the ZY3-02 data.GCPs were set evenly throughout the study area referencing the ZY3-02 multispectral image and ESRI's online World Imagery (Environmental Systems Research Institute, Inc., United States).In this study, we extracted the elevation of GCPs from the SRTMGL1 to minimize fieldwork.However, this brought some difficulties to locating the GCPs.Most of the surface features in the study area were located in valleys, buildings, and at the intersections of forest belts, where elevation changes occur.The values of the SRTMGL1 at these locations may be affected by nearby elevation changes, since the spatial resolution of the SRTMGL1 is about 27 m in the study area.This may introduce bias to the GCPs.To avoid this problem, the GCPs were placed in flat regions and their locations were determined by the intersections of the lines connecting surface features.Figure 3 illustrates the scheme for setting GCPs on ESRI's online World Imagery layer.Figure 4 shows the locations of all 49 GCPs on the true color composite ZY3-02 multispectral image.In the extraction of the DSM, we first rediscovered the four surface features for each GCP in the forward and backward view images, then located the GCPs by connecting the features.
ISPRS Int.J. Geo-Inf.2019, 8 FOR PEER REVIEW 5 points (TIEs) to compute a math model that relates the rows and columns of the matched pixels with ground coordinates and elevations.The polynomial coefficients were provided in the rational polynomial coefficients (RPC) files distributed with the ZY3-02 data.GCPs were set evenly throughout the study area referencing the ZY3-02 multispectral image and ESRI's online World Imagery (Environmental Systems Research Institute, Inc., United States).In this study, we extracted the elevation of GCPs from the SRTMGL1 to minimize fieldwork.However, this brought some difficulties to locating the GCPs.Most of the surface features in the study area were located in valleys, buildings, and at the intersections of forest belts, where elevation changes occur.The values of the SRTMGL1 at these locations may be affected by nearby elevation changes, since the spatial resolution of the SRTMGL1 is about 27 m in the study area.This may introduce bias to the GCPs.To avoid this problem, the GCPs were placed in flat regions and their locations were determined by the intersections of the lines connecting surface features.Figure 3 illustrates the scheme for setting GCPs on ESRI's online World Imagery layer.Figure 4 shows the locations of all 49 GCPs on the true color composite ZY3-02 multispectral image.In the extraction of the DSM, we first rediscovered the four surface features for each GCP in the forward and backward view images, then located the GCPs by connecting the features.In addition, we collected a total of 226 TIEs interactively in the forward and backward view images.The Semi-Global Matching (SGM) method [36] was used to match pixels.In order to avoid points (TIEs) to compute a math model that relates the rows and columns of the matched pixels with ground coordinates and elevations.The polynomial coefficients were provided in the rational polynomial coefficients (RPC) files distributed with the ZY3-02 data.GCPs were set evenly throughout the study area referencing the ZY3-02 multispectral image and ESRI's online World Imagery (Environmental Systems Research Institute, Inc., United States).In this study, we extracted the elevation of GCPs from the SRTMGL1 to minimize fieldwork.However, this brought some difficulties to locating the GCPs.Most of the surface features in the study area were located in valleys, buildings, and at the intersections of forest belts, where elevation changes occur.The values of the SRTMGL1 at these locations may be affected by nearby elevation changes, since the spatial resolution of the SRTMGL1 is about 27 m in the study area.This may introduce bias to the GCPs.To avoid this problem, the GCPs were placed in flat regions and their locations were determined by the intersections of the lines connecting surface features.Figure 3 illustrates the scheme for setting GCPs on ESRI's online World Imagery layer.Figure 4 shows the locations of all 49 GCPs on the true color composite ZY3-02 multispectral image.In the extraction of the DSM, we first rediscovered the four surface features for each GCP in the forward and backward view images, then located the GCPs by connecting the features.In addition, we collected a total of 226 TIEs interactively in the forward and backward view images.The Semi-Global Matching (SGM) method [36] was used to match pixels.In order to avoid  In addition, we collected a total of 226 TIEs interactively in the forward and backward view images.The Semi-Global Matching (SGM) method [36] was used to match pixels.In order to avoid the loss of precision, the sampling distance of the output DSM was set to 2 m.In order to preserve all the details, we did not filter the DSM.

Sample Selection
We selected 51 pairs of canopy/ground samples in flat areas through the manual interpretation of the ZY3-02 multispectral image (Figure 5a,c).Areas affected by shadows and mismatches during the DSM extraction process were excluded from the sample selection, referencing the ZY3-02 DSM displayed with hill shade effect (Figure 5b,d).Canopy samples were set in relatively uniform and closed forests.The elevation undulation of the surrounding bare ground needed to be within 2 m (i.e., the vertical accuracy of ZY3-02) so that the sampled forests were unlikely to be on slopes.The ground samples were set in flat areas without trees, buildings, or ditches and as close as possible to the canopy samples, to make the elevation of the ground samples close to the elevation of the ground beneath the sampled forest canopies.The attributes of the sampling areas are summarized in Table 1.We averaged the DSM in each sampling area and obtained the canopy height samples by calculating the difference of the averaged DSM elevation for each pair of canopy/ground samples.The canopy height samples had a minimum height of 0.5 m, a maximum height of 14.3 m, and an average height of 6.1 m.
the loss of precision, the sampling distance of the output DSM was set to 2 m.In order to preserve all the details, we did not filter the DSM.

Sample Selection
We selected 51 pairs of canopy/ground samples in flat areas through the manual interpretation of the ZY3-02 multispectral image (Figure 5a,c).Areas affected by shadows and mismatches during the DSM extraction process were excluded from the sample selection, referencing the ZY3-02 DSM displayed with hill shade effect (Figure 5b,d).Canopy samples were set in relatively uniform and closed forests.The elevation undulation of the surrounding bare ground needed to be within 2 m (i.e. the vertical accuracy of ZY3-02) so that the sampled forests were unlikely to be on slopes.The ground samples were set in flat areas without trees, buildings, or ditches and as close as possible to the canopy samples, to make the elevation of the ground samples close to the elevation of the ground beneath the sampled forest canopies.The attributes of the sampling areas are summarized in Table 1.We averaged the DSM in each sampling area and obtained the canopy height samples by calculating the difference of the averaged DSM elevation for each pair of canopy/ground samples.The canopy height samples had a minimum height of 0.5 m, a maximum height of 14.3 m, and an average height of 6.1 m.

Canopy Height Modeling and Extrapolation
The canopy height samples were randomly divided into a training set and a validation set, with 34 samples in the training set and 17 samples in the validation set (Figure 5e).We first performed geometric correction to the ZY3-02 multispectral image using the ZY3-02 DSM.Then, we performed radiometric correction to obtain the surface reflectance of each band through the radiometric correction workflow in PCI Geomatica software (PCI Geomatics Enterprises, Inc., Canada).The atmospheric condition of the input image was estimated based on the season when the image was captured and the location at which the image was taken.The workflow works with a database of atmospheric correction functions stored in lookup tables for different altitude profiles of pressure,

Canopy Height Modeling and Extrapolation
The canopy height samples were randomly divided into a training set and a validation set, with 34 samples in the training set and 17 samples in the validation set (Figure 5e).We first performed geometric correction to the ZY3-02 multispectral image using the ZY3-02 DSM.Then, we performed radiometric correction to obtain the surface reflectance of each band through the radiometric correction workflow in PCI Geomatica software (PCI Geomatics Enterprises, Inc., Canada).The atmospheric condition of the input image was estimated based on the season when the image was captured and the location at which the image was taken.The workflow works with a database of atmospheric correction functions stored in lookup tables for different altitude profiles of pressure, humidity, and aerosol type [37].We built a multiple linear regression model to correlate the mean surface reflectance of the multispectral image in sampling areas with the canopy heights derived from the ZY3-02 DSM.The modeling was implemented in SPSS software (International Business Machines Corp., United States) and a stepwise method was used to select variables [32].
We classified the ZY3-02 multispectral image using the maximum likelihood supervised classification algorithm.The algorithm assumes that the statistics for each class in each band are normally distributed and calculates the probability that a given pixel belongs to a specific class.Each pixel is assigned to the class that has the highest probability (i.e., the maximum likelihood) based on the probability density functions estimated from the training samples.The classification was implemented through classification workflow using ENVI software (Harris Geospatial Solutions, Inc., United States).Training samples were selected in the multispectral image through manual interpretation.The number of samples for forest, water, bare surfaces, and farmland were 80, 20, 64, and 58, respectively.The classification results were exported as raster files without smoothing or aggregation.
After that, the surface reflectance of the multispectral image in forest regions was input into the regression model and a wall-to-wall forest canopy height map for the forest regions throughout the study area was obtained.

Results
Figure 6 shows the comparison of the ZY3-02 DSM and SRTMGL1 at the locations of the GCPs.The coefficient of determination (R 2 ), mean error (µ), and standard deviation of errors (σ) were 0.9992, 0.04 m, and 1.91 m, respectively.
ISPRS Int.J. Geo-Inf.2019, 8 FOR PEER REVIEW 7 humidity, and aerosol type [37].We built a multiple linear regression model to correlate the mean surface reflectance of the multispectral image in sampling areas with the canopy heights derived from the ZY3-02 DSM.The modeling was implemented in SPSS software (International Business Machines Corp., United States) and a stepwise method was used to select variables [32].We classified the ZY3-02 multispectral image using the maximum likelihood supervised classification algorithm.The algorithm assumes that the statistics for each class in each band are normally distributed and calculates the probability that a given pixel belongs to a specific class.Each pixel is assigned to the class that has the highest probability (i.e., the maximum likelihood) based on the probability density functions estimated from the training samples.The classification was implemented through classification workflow using ENVI software (Harris Geospatial Solutions, Inc., United States).Training samples were selected in the multispectral image through manual interpretation.The number of samples for forest, water, bare surfaces, and farmland were 80, 20, 64, and 58, respectively.The classification results were exported as raster files without smoothing or aggregation.
After that, the surface reflectance of the multispectral image in forest regions was input into the regression model and a wall-to-wall forest canopy height map for the forest regions throughout the study area was obtained.

Results
Figure 6 shows the comparison of the ZY3-02 DSM and SRTMGL1 at the locations of the GCPs.The coefficient of determination (R 2 ), mean error (μ), and standard deviation of errors (σ) were 0.9992, 0.04 m, and 1.91 m, respectively.After the stepwise process, the red band and green band were selected as the predictors of the optimal model.The model was of the form where H is the forest canopy height, Red is the surface reflectance of the red band, and Green is the surface reflectance of the green band.The model and predictors were statistically significant, with the p-values of the F-test (for the model) and t-test (for the predictors) being less than 0.01.The model had an R 2 of 0.67 and a root mean square error (RMSE) of 1.55 m (Figure 7).After the stepwise process, the red band and green band were selected as the predictors of the optimal model.The model was of the form H = 1.3219 − 0.1013 * Red + 0.0713 * Green (1) where H is the forest canopy height, Red is the surface reflectance of the red band, and Green is the surface reflectance of the green band.The model and predictors were statistically significant, with the p-values of the F-test (for the model) and t-test (for the predictors) being less than 0.01.The model had an R 2 of 0.67 and a root mean square error (RMSE) of 1.55 m (Figure 7).The accuracy of the landcover classification was evaluated by 500 points randomly set in the study area.The reference categories of each point were determined by manual interpretation and were used to validate the results of the supervised classification.The confusion matrix showed an overall accuracy of 84%.The predicted forest canopy height map of the entire study area is shown in Figure 8a.To validate predicted forest canopy height, we calculated the mean heights in areas of the validation samples and compared them with the DSM-derived canopy heights.The validation reported an R 2 of 0.72 and a RMSE of 1.58 m (Figure 9).The accuracy of the landcover classification was evaluated by 500 points randomly set in the study area.The reference categories of each point were determined by manual interpretation and were used to validate the results of the supervised classification.The confusion matrix showed an overall accuracy of 84%.The predicted forest canopy height map of the entire study area is shown in Figure 8a.The accuracy of the landcover classification was evaluated by 500 points randomly set in the study area.The reference categories of each point were determined by manual interpretation and were used to validate the results of the supervised classification.The confusion matrix showed an overall accuracy of 84%.The predicted forest canopy height map of the entire study area is shown in Figure 8a.To validate the predicted forest canopy height, we calculated the mean heights in areas of the validation samples and compared them with the DSM-derived canopy heights.The validation reported an R 2 of 0.72 and a RMSE of 1.58 m (Figure 9).To validate the predicted forest canopy height, we calculated the mean heights in areas of the validation samples and compared them with the DSM-derived canopy heights.The validation reported an R 2 of 0.72 and a RMSE of 1.58 m (Figure 9).

Discussion
The use of the intersections of the lines connecting surface features showed good performance in determining the locations of the GCPs.The standard deviation of 1.91 m in Figure 6 was slightly better than the vertical accuracy of 2 m for ZY3-02 reported by Xu et al. [33].The canopy samples and ground samples were selected through manual interpretation to ensure the quality of canopy height samples.After selecting the samples, it was necessary to correlate the reflectance of the multispectral image with the DSM-derived canopy heights through a statistical model in order to obtain a wall-towall canopy height map.Since the forests in the study area were relatively short, a linear model would have been able to effectively characterize the relationship between the reflectance and canopy height [32].However, in tall forests (e.g.taller than 30 m), the relationship between reflectance and canopy height may deviate from linear form and a nonlinear model may be required in that case [9].The reflectance of the red band contributed most in the estimation of canopy height, which was consistent with the study of Hansen et al. [9].
We did not use field data to validate our results.However, the accuracy of the DSM-derived canopy height samples should be consistent with the vertical accuracy of ZY3-02, which is better than 2 m based on the evaluations conducted by Xu et al. [33].In addition, we used the average DSM elevation in sampling areas to calculate the canopy height samples.This further reduced the likelihood of error and made the accuracy of the DSM-derived canopy heights closer to that of the field measurements [38].
The results in Figure 9 were acceptable considering the planimetric and vertical accuracies of the stereo images and the limited bands of the multispectral image.However, in order to better meet the requirements of forestry applications, further improvements would be necessary in future.Finer and better stereo images are worth trying to obtain, if available.For example, China's Gaofen-7 stereo mapping satellite, due for launch in 2019, provides stereo images with sub-meter spatial resolution and can carry out topographic mapping on a scale of 1:10,000.It would also be helpful to use more spectral bands such as the short-wave near-infrared band, which has a strong response to the water content of canopies.With more bands, various vegetation indices can be used in the modeling, which will help to reduce the effects of illumination conditions and canopy surface undulations.Finally, machine learning algorithms may be useful because they usually do not require assumptions about the form of the model, and this will also be the direction of our future work.

Conclusions
In this study, we mapped the canopy heights of poplar plantations in plain areas using ZY3-02 data by combining photogrammetry with statistical modeling.We found that using the intersections of the lines connecting surface features to determine the locations of the ground control points

Discussion
The use of the intersections of the lines connecting surface features showed good performance in determining the locations of the GCPs.The standard deviation of 1.91 m in Figure 6 was slightly better than the vertical accuracy of 2 m for ZY3-02 reported by Xu et al. [33].The canopy samples and ground samples were selected through manual interpretation to ensure the quality of canopy height samples.After selecting the samples, it was necessary to correlate the reflectance of the multispectral image with the DSM-derived canopy heights through a statistical model in order to obtain a wall-to-wall canopy height map.Since the forests in the study area were relatively short, a linear model would have been able to effectively characterize the relationship between the reflectance and canopy height [32].However, in tall forests (e.g.taller than 30 m), the relationship between reflectance and canopy height may deviate from linear form and a nonlinear model may be required in that case [9].The reflectance of the red band contributed most in the estimation of canopy height, which was consistent with the study of Hansen et al. [9].
We did not use field data to validate our results.However, the accuracy of the DSM-derived canopy height samples should be consistent with the vertical accuracy of ZY3-02, which is better than 2 m based on the evaluations conducted by Xu et al. [33].In addition, we used the average DSM elevation in sampling areas to calculate the canopy height samples.This further reduced the likelihood of error and made the accuracy of the DSM-derived canopy heights closer to that of the field measurements [38].
The results in Figure 9 were acceptable considering the planimetric and vertical accuracies of the stereo images and the limited bands of the multispectral image.However, in order to better meet the requirements of forestry applications, further improvements would be necessary in future.Finer and better stereo images are worth trying to obtain, if available.For example, China's Gaofen-7 stereo mapping satellite, due for launch in 2019, provides stereo images with sub-meter spatial resolution and can carry out topographic mapping on a scale of 1:10,000.It would also be helpful to use more spectral bands such as the short-wave near-infrared band, which has a strong response to the water content of canopies.With more bands, various vegetation indices can be used in the modeling, which will help to reduce the effects of illumination conditions and canopy surface undulations.Finally, machine learning algorithms may be useful because they usually do not require assumptions about the form of the model, and this will also be the direction of our future work.

Conclusions
In this study, we mapped the canopy heights of poplar plantations in plain areas using ZY3-02 data by combining photogrammetry with statistical modeling.We found that using the intersections of the lines connecting surface features to determine the locations of the ground control points achieved good performance.It avoided the interferences from elevation changes at the surface features when extracting the elevation of GCPs from Shuttle Radar Topography Mission global 1 arc second data.It was a simple and effective way to obtain canopy height samples by selecting canopy samples and ground samples through manual interpretation and calculating their elevation differences.The red band and green band were selected by a stepwise method in the establishing of a multiple linear regression model.The validation of the model-predicted canopy height indicated a coefficient of determination of 0.72 and a root mean square error of 1.58 m.The proposed method extends the application of ZY3-02 data to mapping canopy heights of pure plantations in plain areas.

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

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

Figure 3 .
Figure 3. Scheme for setting GCPs on ESRI's online World Imagery layer (Environmental Systems Research Institute, Inc., United States).Subgraphs (a-c) show the locations of GCPs determined by intersections of forest belts, ditches in farmlands, and edges of buildings and valleys, respectively.

Figure 4 .
Figure 4. Locations of all 49 GCPs on the true color composite ZY3-02 multispectral image.

Figure 3 .
Figure 3. Scheme for setting GCPs on ESRI's online World Imagery layer (Environmental Systems Research Institute, Inc., United States).Subgraphs (a-c) show the locations of GCPs determined by intersections of forest belts, ditches in farmlands, and edges of buildings and valleys, respectively.

Figure 3 .
Figure 3. Scheme for setting GCPs on ESRI's online World Imagery layer (Environmental Systems Research Institute, Inc., United States).Subgraphs (a-c) show the locations of GCPs determined by intersections of forest belts, ditches in farmlands, and edges of buildings and valleys, respectively.

Figure 4 .
Figure 4. Locations of all 49 GCPs on the true color composite ZY3-02 multispectral image.

Figure 4 .
Figure 4. Locations of all 49 GCPs on the true color composite ZY3-02 multispectral image.

Figure 5 .
Figure 5. Scheme of sample selection.Subgraphs (a-d) illustrate the setting of canopy samples and ground samples referencing the false color composite ZY3-02 multispectral image and the ZY3-02 DSM displayed with hill shade effect.Subgraph (e) shows the locations of all 51 canopy height samples (which were also the locations of the canopy samples) and marks the training samples and validation samples.

Figure 5 .
Figure 5. Scheme of sample selection.Subgraphs (a-d) illustrate the setting of canopy samples and ground samples referencing the false color composite ZY3-02 multispectral image and the ZY3-02 DSM displayed with hill shade effect.Subgraph (e) shows the locations of all 51 canopy height samples (which were also the locations of the canopy samples) and marks the training samples and validation samples.

Figure 6 .
Figure 6.Comparison of the ZY3-02 DSM and SRTMGL1 at the locations of the GCPs.

Figure 6 .
Figure 6.Comparison of the ZY3-02 DSM and SRTMGL1 at the locations of the GCPs.

Figure 7 .
Figure 7.Comparison of the model-predicted heights and the DSM-derived heights for the training samples.Legend: R 2 , coefficient of determination; RMSE, root mean square error.
Figure 8b,c show the canopy height maps of two regions corresponding to the subgraphs (a,b) and (c,d) in Figure 5.

Figure 8 .
Figure 8. Subgraph (a) shows the predicted forest canopy height map of the entire study area.Subgraphs (b,c) show the canopy height maps of two regions corresponding to the subgraphs (a,b) and (c,d) in Figure 5.

7 .
Comparison of the model-predicted heights and the DSM-derived heights for the training samples.Legend: R 2 , coefficient of determination; RMSE, root mean square error.
Figure 8b,c show the canopy height maps of two regions corresponding to the subgraphs (a,b) and (c,d) in Figure 5.

Figure 7 .
Figure 7.Comparison of the model-predicted heights and the DSM-derived heights for the training samples.Legend: R 2 , coefficient of determination; RMSE, root mean square error.
Figure 8b,c show the canopy height maps of two regions corresponding to the subgraphs (a,b) and (c,d) in Figure 5.

Figure 8 .
Figure 8. Subgraph (a) shows the predicted forest canopy height map of the entire study area.Subgraphs (b,c) show the canopy height maps of two regions corresponding to the subgraphs (a,b) and (c,d) in Figure 5.

Figure 8 .
Figure 8. Subgraph (a) shows the predicted forest canopy height map of the entire study area.Subgraphs (b,c) show the canopy height maps of two regions corresponding to the subgraphs (a,b) and (c,d) in Figure 5.

Figure 9 .
Figure 9.Comparison of the model-predicted heights and the DSM-derived heights for the validation samples.

Figure 9 .
Figure 9.Comparison of the model-predicted heights and the DSM-derived heights for the validation samples.

Table 1 .
Area, DSM pixel counts, and mean DSM elevation of the sampling areas.

Table 1 .
Area, DSM pixel counts, and mean DSM elevation of the sampling areas.