Estimates of Forest Canopy Height Using a Combination of ICESat-2 / ATLAS Data and Stereo-Photogrammetry

: Forest canopy height is an indispensable forest vertical structure parameter for understanding the carbon cycle and forest ecosystem services. A variety of studies based on spaceborne Lidar, such as ICESat, ICESat-2 and airborne Lidar, were conducted to estimate forest canopy height at multiple scales. However, while a few studies have been conducted based on ICESat-2 simulated data from airborne Lidar data, few studies have analyzed ATL08 and ATL03 products derived from the ATLAS sensor onboard ICESat-2 for regional vegetation canopy height mapping. It is necessary and promising to explore how data obtained by ICESat-2 can be applied to estimate forest canopy height. This study proposes a new means to estimate forest canopy height, deﬁned as the mean height of trees within a given forest area, using a combination of ICESat-2 ATL08 and ATL03 data and ZY-3 satellite stereo images. Five procedures were used to estimate the forest canopy height of the city of Nanning in China: (1) Processing ground photons in a 30 m × 30 m grid; (2) Extracting a digital surface model (DSM) using ZY-3 stereo images; (3) Calculating a discontinuous canopy height model (CHM) dataset; (4) Validating the DSM and ground photon height using GEDI data; (5) Estimating the regional wall-to-wall forest canopy height product based on the backpropagation artiﬁcial neural network (BP-ANN) model and Landsat 8 vegetation indices and independent accuracy assessments with ﬁeld measured plots. The validation shows a root mean square error (RMSE) of 3.34 m to 3.47 m and a coe ﬃ cient of determination R 2 = 0.51. The new method shows promise and can be used for large-scale forest canopy height mapping at various resolutions or in combination with other data, such as SAR images. Finally, this study analyzes resolutions and how to ﬁlter e ﬀ ective data when ATL08 data are directly used to generate regional or global vegetation height products, which will be the focus of future research. X.L.; X.L.; B.B.; visualization, B.B. and Y.D.; supervision, C.C.; M.X.,


Introduction
As essential terrestrial ecosystems, forests form an important part of the earth's carbon cycle [1,2]. Forests not only play an indispensable role in various ecological and social services but also help maintain carbon balance and mitigate climate change, as they are considered carbon sinks [3][4][5][6][7]. The total carbon storage of forest is equivalent to 85% of total terrestrial stocks and 75% of total Montesano explored the abilities and uncertainties of ICESat-2 photon counting data in estimating coniferous forest biomass according to different biomass gradients and sampling distances along the orbit [34]. The authors found the 50 m scale to be the optimal horizontal resolution of ICESat-2 model fitting, and the biomass level was 20 Mg/ha with the estimation error of 20% to 50%. Gwenzi evaluated the potential of ATLAS data use for canopy height extraction from savanna ecosystems [14]. The authors found the correlation between canopy height extracted from the MATLAS simulator, which generates ATLAS-like data, and that from discrete return Lidar (DRL) to be weaker than that from Multiple Altimeter Beam Experimental Lidar (MABEL, an airborne photon counting Lidar sensor) data and predicted that the number of signal photons in spaceborne ATLAS data would be greatly reduced, causing the accuracy of canopy height estimation to decrease. Narine used simulated ICESat-2 photon counting Lidar data to generate point cloud data with the same along-orbit segmentation length as the ICESat-2 land vegetation product (ATL08), 100 m [15]. The relationship between photon counting vegetation products and forest AGB and canopy coverage obtained by airborne Lidar was analyzed in no-noise scenes, daytime scenes and night scenes [15]. Neuenschwander presented ATLAS data and verified the accuracy of its ground elevation [35][36][37]. The RMSE value of ground elevation data obtained from the ATLAS is 0.85 m, which is lower than that of the Space Shuttle Radar topographic survey mission (SRTM).
However, few studies have used ICESat-2 vegetation height products to extrapolate and map regional vegetation canopy height [38,39], and most of them are based on ICESat-2 simulation data of airborne Lidar. The purpose of this paper is to explore ways to use the ATL08 and ATL03 products to map regional forest canopy heights with a 30 m resolution. There are many definitions of forest canopy height. Forest canopy height in this study refers to the mean height of trees within a given forest area. In this study, ICESat-2 ATLAS photon data and ZY-3 stereo image pair data are combined to generate a discontinuous CHM dataset. GEDI data are used to verify the accuracy of ground photon elevation and the DSM. The CHM is used to create a training sample, and the BP-ANN model is used to generate a forest canopy height map of the city of Nanning. Finally, this study analyzes resolution and how to filter effective data that emerge when ATL08 data are directly used to generate regional or global vegetation height products, which will be the focus of future research.

Description of the Study Area
Nanning, located in the south-central area of Guangxi, is the capital of the Guangxi Zhuang Autonomous Region. Geographically, Nanning is located between 22 • 13 and 23 • 32 north and 107 • 45 and 108 • 51 east with a total land area of 22,100 km 2 ( Figure 1). Nanning consists of eight administrative regions (municipal districts and the Yongning, Mashan, Binyang, Longan, Wuming, Shanglin, and Heng counties). Elevation in the study area ranges from 70 m to 600 m, and flat land is the main geomorphic type, accounting for 57.78% of the total area. The climate in Nanning is a subtropical monsoon with average annual precipitation levels ranging from 1241 mm to 1753 mm, and the average annual temperature is approximately 21.6 • C [40]. Additionally, the presence of abundant water and heat resources render Nanning rich in forest resources. The main tree species present are eucalyptus, masson pine and cunninghamia lanceolata. According to the results of the fifth Guangxi forest resources survey, in 2016, the forested area of Nanning reached 1,050,000 hectares, and the forest volume reached 51.99 million cubic meters (http://search.nanning.gov.cn/).

ICESat-2 ATLAS
The Advanced Topographic Laser Altimeter System (ATLAS) aboard ICESat-2 emits three pairs of beams with a spacing of 3.3 km and with a distance of approximately 90 m between the two beams. A pair consists of a strong and a weak beam. The diameter of the footprint is 14 m to 17 m and, the sampling interval along the track is set to 70 cm to facilitate intensive sampling and better capture elevation changes [29]. ATL01 to ATL21 contain three levels of ICESat-2 mission products. The Level-2 product of global geolocated photons (ATL03) measures the latitude, longitude, time and absolute height above the World Geodetic System 1984 (WGS84) for each photon [41]. ATL03 provides photon information for all Level-3A products such as land, ice elevation (ATL06), Arctic/Antarctic sea ice elevation (ATL07), land, water and vegetation elevation (ATL08) data [42]. The corresponding data of three pairs of laser beams are numbered gt1l-gt3r, representing the six ground tracks (GT) from left to right. ATL08 provides terrain, canopy height, and canopy cover parameters at an along-track segment of 100 m, while the 100 m segment consists of 5 sequential 20 m segments labeled with segment IDs in ATL03 and ATL08 [43]. The canopy height was calculated by the 98% height of all individual relative canopy heights, which is derived from a surface fitting line and vegetation surface fitting line within a 100 m segment. ATL08 provides a photon classification label (Noise, Ground, Canopy, and Top of Canopy), and the labeling of ground photons is crucial for this study. In this study, ATL03 and ATL08 data were utilized in combination where ATL03 provided height and coordinate information while ATL08 provided photon classification information. ATL03 and ATL08 data were downloaded from the National Snow and Ice Data Center (NSIDC, https://nsidc.org/data/ATL03/versions/1,

ICESat-2 ATLAS
The Advanced Topographic Laser Altimeter System (ATLAS) aboard ICESat-2 emits three pairs of beams with a spacing of 3.3 km and with a distance of approximately 90 m between the two beams. A pair consists of a strong and a weak beam. The diameter of the footprint is 14 m to 17 m and, the sampling interval along the track is set to 70 cm to facilitate intensive sampling and better capture elevation changes [29].
ATL01 to ATL21 contain three levels of ICESat-2 mission products. The Level-2 product of global geolocated photons (ATL03) measures the latitude, longitude, time and absolute height above the World Geodetic System 1984 (WGS84) for each photon [41]. ATL03 provides photon information for all Level-3A products such as land, ice elevation (ATL06), Arctic/Antarctic sea ice elevation (ATL07), land, water and vegetation elevation (ATL08) data [42]. The corresponding data of three pairs of laser beams are numbered gt1l-gt3r, representing the six ground tracks (GT) from left to right. ATL08 provides terrain, canopy height, and canopy cover parameters at an along-track segment of 100 m, while the 100 m segment consists of 5 sequential 20 m segments labeled with segment IDs in ATL03 and ATL08 [43]. The canopy height was calculated by the 98% height of all individual relative canopy heights, which is derived from a surface fitting line and vegetation surface fitting line within a 100 m segment. ATL08 provides a photon classification label (Noise, Ground, Canopy, and Top of Canopy), and the labeling of ground photons is crucial for this study. In this study, ATL03 and ATL08 data were utilized in combination where ATL03 provided height and coordinate information while ATL08 provided photon classification information. ATL03 and ATL08 data were downloaded from the National Snow and Ice Data Center (NSIDC, https://nsidc.org/data/ATL03/versions/1, https://nsidc.org/data/ATL08/versions/1) (Tables 1 and 2). ICESat-2 ATL03 and ATL08 data for the whole study area were acquired from October  30, 2018 to April 9, 2019, and data overlap with the locations of ZY-3 data for October 30, 2018 to  November 03, 2018 (Table 1, Figure 2).

ZY-3 Data
The ZY-3 satellite, the first civil cartographic satellite designed for stereo photogrammetry in China launched in 2012, can measure ground heights and is mainly used for 1:50,000 topographic maps. Applying ZY-3 stereo image data for forestry research can generate a large area continuous digital surface model (DSM), which can effectively prevent errors caused by the spatial interpolation of discrete points. ZY-3 carries three panchromatic sensors pointing in the nadir, backward, and forward directions with resolutions of 2.1 m × 2.1 m, two 3.5 m × 3.5 m, respectively, and one multispectral sensor that consists of four bands (near infrared, red, green, and blue) with a resolution of 5.8 m × 5.8 m [44]. The ZY-3 02 satellite launched in 2016 is an upgraded successor of ZY-3 that have improved resolution of backward and forward panchromatic sensors from 3.5 m × 3.5 m to 2.7  The ZY-3 satellite, the first civil cartographic satellite designed for stereo photogrammetry in China launched in 2012, can measure ground heights and is mainly used for 1:50,000 topographic maps. Applying ZY-3 stereo image data for forestry research can generate a large area continuous digital surface model (DSM), which can effectively prevent errors caused by the spatial interpolation of discrete points. ZY-3 carries three panchromatic sensors pointing in the nadir, backward, and forward directions with resolutions of 2.1 m × 2.1 m, two 3.5 m × 3.5 m, respectively, and one multispectral sensor that consists of four bands (near infrared, red, green, and blue) with a resolution of 5.8 m × 5.8 m [44]. The ZY-3 02 satellite launched in 2016 is an upgraded successor of ZY-3 that have improved resolution of backward and forward panchromatic sensors from 3.5 m × 3.5 m to 2.7 m × 2.7 m [45]. ZY-3 and ZY-3 02 have the same ground swath of approximately 50 km × 50 km, and the same radiometric resolution of 10 bits [46].
In this study, the ZY-3 and ZY-3 02 image data were downloaded from the China Centre for Resources Satellite Data and Application (CRESDA, http://www.cresda.com/CN/) with the acquisition period set to 10 March 2018 to 5 October 2018 considering data availability and cloud cover (Table 1). Due to the presence of considerably more mountain shadows in forward view scenes, nadir and backward view panchromatic scenes were utilized to extract the DSM while using multispectral scenes to extract the forest area.

Landsat 8 Operational Land Imager Data
Landsat 8 with the Operational Land Imager (OLI) sensor provides satellite images with wavelengths of 0.43 µm to 1.38 µm. OLI images were obtained from the United States Geological Survey (USGS, https://glovis.usgs.gov/) on 31 October 2018 and 6 October 2018 ( Table 1). The three scenes show less than 5% cloud cover. Bands 1-7, which contain four visible bands, one near infrared band and two shortwave infrared bands with a resolution of 30 m × 30 m, were utilized to extract the forest area by maximum likelihood classification (MLC) to calculate the vegetation indices and to estimate the forest canopy height of the study area. Preprocessing operations used included geometric correction, radiometric calibration, and atmospheric correction by the Fast Line-of-Sight Atmospheric Analysis of Hypercubes (FLAASH) module in ENVI to eliminate atmospheric effects.
This study refers to the existing literature, and a total of ten indices are selected for analysis [44,47,48], including MSR, SLAVI, NDVI, and others (Table 3).

GEDI
The Global Ecosystem Dynamics Investigation (GEDI), the first spaceborne Lidar system designed specifically for studies of forest structures and launched in 2018, is an Earth Venture Instrument (EVI) developed by NASA and equipped with three lasers and full waveform recording Lidar instrument aboard the International Space Station (ISS) [56,57]. The GEDI operates at 1064 nm and produces eight beam ground tracks spaced approximately 600 m apart in the cross-track direction and sampling approximately 60 m along a track from a 4.2 km wide swath. Each footprint reflects Lidar energy with a diameter of approximately 30 m and a vertical measurement accuracy level of approximately 2-3 cm [57,58]. GEDI L2A product data with ground elevation and relative height (RH) metrics were taken from NASA's Land Processes Distributed Active Archive Center (LPDAAC, https://e4ftl01.cr.usgs.gov/GEDI/GEDI02_A.001/) using GEDI Finder to locate GEDI orbits (files) that intersect with the study area. GEDI data for 29 April 2019 to 1 October 2019 were acquired as the only GEDI data available for the study area when this study was conducted ( Table 2). The GEDI ground tracks are displayed in Figure 3.

Field Plots and Auxiliary Data
A total of 66 field plots for 15 to 21 January 2018 were obtained to implement an independent accuracy assessment of estimated canopy height. Every field plot covers a 30 m × 30 m square area with information on measured mean canopy height and other forest structure parameters ( Figure 1). Considering the economic costs and availability of field survey data, some of the field measurement points are scattered across the study area and the rest are concentrated within the Gaofeng forest farm. All trees with a diameter at breast height (DBH) of greater than 5 cm in each plot were recorded. The forest canopy height of each tree was derived from distance and angle measurements. Distance was measured by a hand-held laser range finder, and angles were measured by a mechanical clinometer. The forest canopy height of the recorded trees was used to calculate the mean canopy height of each square. The coordinates of the southwest and northeast corners of each plot were recorded by GPS, and the center points of the two recorded points were calculated and was taken as the plot coordinates. SRTMGL1 data were used in collecting ground control points (GCPs) during the extraction of the DSM with ZY-3 data and were taken from NASA. The georeferenced image provides accurate coordinate information used for ZY-3 and Landsat 8 in the geometric correction.

Methods
A new canopy height estimation method was developed by combining the ICESat-2 ATLAS data with ZY-3 stereo images. Figure 4 shows the flowchart for this study which consists of five sections: acquired as the only GEDI data available for the study area when this study was conducted ( Table  2). The GEDI ground tracks are displayed in Figure 3.

Field Plots and Auxiliary Data
A total of 66 field plots for 15 to 21 January 2018 were obtained to implement an independent accuracy assessment of estimated canopy height. Every field plot covers a 30 m × 30 m square area with information on measured mean canopy height and other forest structure parameters (Figure 1). Considering the economic costs and availability of field survey data, some of the field measurement points are scattered across the study area and the rest are concentrated within the Gaofeng forest farm. All trees with a diameter at breast height (DBH) of greater than 5 cm in each plot were recorded. The forest canopy height of each tree was derived from distance and angle measurements. Distance was measured by a hand-held laser range finder, and angles were measured by a mechanical clinometer. The forest canopy height of the recorded trees was used to calculate the mean canopy height of each square. The coordinates of the southwest and northeast corners of each plot were recorded by GPS, and the center points of the two recorded points were calculated and was taken as the plot coordinates. SRTMGL1 data were used in collecting ground control points (GCPs) during the extraction of the DSM with ZY-3 data and were taken from NASA. The georeferenced image provides accurate coordinate information used for ZY-3 and Landsat 8 in the geometric correction.

Processing ICESat-2/ATLAS ATL08 and ATL03 Data
Light beams from the ATLAS laser hits the earth's surface, and a handful of photons are then reflected back to the laser. The number of returning photons depends on the outgoing energy, solar and atmospheric conditions and on surface reflectance. The number of ground photons is significantly greater than that of vegetation because the reflectance of the terrain surface is typically approximately 0.3, which is higher than 0.1 of vegetation [43]. The ATL03 data product not only provides basic information on photons (latitude, longitude, and height of the WGS84 ellipsoid) but also extracts signal photons from noise photons. In producing the ATL08 data product, further filtering techniques to remove noise photons, ground finding filter, and canopy top filter were conducted orderly to label all the signal photons in ground, canopy and noise, respectively.
In this study, ground photons were used as the DEM to provide land surface height information. Densely distributed ground photons were resampled into a 30 m × 30 m grid with an average ground photon height ( Figure 5). As a result, a total of 409,515 ground photons were processed into 22,437 new points. Then, if the number of ground photons in a 30 m × 30 m grid was less than 4, the new point was deemed ineffective and was removed since the expected number of signal photons per laser for the vegetated surface is 0 to 4 photons [43]. For example, "p5" shown in Figure 5c was removed.
A new canopy height estimation method was developed by combining the ICESat-2 ATLAS data with ZY-3 stereo images. Figure 4 shows the flowchart for this study which consists of five sections: (1) Calculating the elevation of ground photons in 30 m × 30 m pixels; (2) Extracting the ZY-3 DSM; (3) Producing a discontinuous canopy height model (CHM) dataset after validation; (4) Estimating the forest canopy height by the BP-ANN model; and (5) Assessing the accuracy of the estimated canopy height map. Light beams from the ATLAS laser hits the earth's surface, and a handful of photons are then reflected back to the laser. The number of returning photons depends on the outgoing energy, solar and atmospheric conditions and on surface reflectance. The number of ground photons is significantly greater than that of vegetation because the reflectance of the terrain surface is typically approximately 0.3, which is higher than 0.1 of vegetation [43]. The ATL03 data product not only provides basic information on photons (latitude, longitude, and height of the WGS84 ellipsoid) but also extracts signal photons from noise photons. In producing the ATL08 data product, further filtering techniques to remove noise photons, ground finding filter, and canopy top filter were conducted orderly to label all the signal photons in ground, canopy and noise, respectively.
In this study, ground photons were used as the DEM to provide land surface height information. Densely distributed ground photons were resampled into a 30 m × 30 m grid with an average ground

Extracting the DSM
The digital surface model (DSM) was derived from ZY-3 and ZY-3 02 stereo images using PCI Geomatica Banff 2019 software with the OrthoEngine module. A nadir and backward-looking stereo pair with an inclination angle of ±22 • -23.5 • between the sensors was used for processing due to the presence of fewer mountain shadows. The OrthoEngine module's main purpose is to compute a math model that associates columns and rows of the matched pixels with ground coordinates (X, Y) and elevations (Z) using rational polynomial coefficients (RPCs), GCPs and tie points (TPs) [45]. The RPC file was distributed with the ZY-3 data file. The Semi-Global Matching (SGM) algorithm was used in the module to match pixels [46]. A total of 92 GCPs and 62 TPs were collected automatically at a variety of elevations with less than two pixel residual errors, meeting the recommended minimum value requirements [59]. The resolution of the output DSM was set to 2 m × 2 m to maintain surface details and preserve precision [44,45].
A two-step mask extraction operation was implemented to obtain the effective DSM: (1) The score band of PCI output can specify the reliability of DSM values, and values of 90 to 100 denote a good pixel; and (2) The forest area of ZY-3 coverage was classified by MLC using ZY-3 multispectral images.
To meet the 30 m × 30 m resolution mapping requirement, the average value of the 30 m × 30 m grid was calculated for subsequent processing.

Calculating a Discontinuous CHM Dataset
In a large number of studies, a CHM is obtained by subtracting DSM and DEM layers. However, the existing high-precision and high-resolution DEM is mostly obtained via ground point cloud fitting using airborne Lidar, and the DSM is obtained from the vegetation surface fitted by stereo image pairs or vegetation point clouds. In this study, the ground photon height of spaceborne Lidar is used as the DEM, making it possible to use the point sampling DEM of spaceborne Lidar, which could be applied for large-scale and high-efficiency CHM sampling and regional or global canopy height mapping.
The DSM derived from ZY-3 data and by subtracting calculated ground photons' average height values provides a discontinuous CHM dataset as canopy height samples. Figure 6 shows the profile of the discontinuous CHM dataset. The heights of all the data were based on the WGS 84 ellipsoid. Red, purple, and back dots represent the ZY-3 DSM, ground photon DEM, and original ground photons, respectively. Distances between red and purple dots denote canopy heights. Similarly, ZY-3 and Landsat 8 forest areas were used as masks to extract effective CHM samples. A total of 1229 CHM points were derived after filtering.

Validating ZY-3 DSM and Ground Photon Values via GEDI
The evaluation of ZY-3 DSM and ground photon average value accuracy is necessary for further analysis. The waveform for each footprint was processed to determine terrain elevations and canopy heights relative to the WGS84 ellipsoid. The elevation of the center of the lowest mode from the Gauss fitting curve of return energy was used to validate the accuracy of the ground photons' average height values while the elevation of the highest detected return was used to validate that of the ZY-3 DSM.
We assume the following: if the grid of 30 m resolution Landsat data in which GEDI footprint points are located is the same as the grid of newly calculated ground photon values, they are considered to be intersected, that is, the same point. To obtain as many intersection points as possible, all GEDI data for April to October 2019 that are currently open access were downloaded. The GEDI DEM corresponding to these intersection points was used to verify the average height values of ground photons. In total, 61 intersection points were extracted (Figure 3). values provides a discontinuous CHM dataset as canopy height samples. Figure 6 shows the profile of the discontinuous CHM dataset. The heights of all the data were based on the WGS 84 ellipsoid. Red, purple, and back dots represent the ZY-3 DSM, ground photon DEM, and original ground photons, respectively. Distances between red and purple dots denote canopy heights. Similarly, ZY-3 and Landsat 8 forest areas were used as masks to extract effective CHM samples. A total of 1229 CHM points were derived after filtering.

Validating ZY-3 DSM and Ground Photon Values via GEDI
The evaluation of ZY-3 DSM and ground photon average value accuracy is necessary for further analysis. The waveform for each footprint was processed to determine terrain elevations and canopy heights relative to the WGS84 ellipsoid. The elevation of the center of the lowest mode from the Gauss fitting curve of return energy was used to validate the accuracy of the ground photons' average height values while the elevation of the highest detected return was used to validate that of the ZY-3 DSM.
We assume the following: if the grid of 30 m resolution Landsat data in which GEDI footprint points are located is the same as the grid of newly calculated ground photon values, they are considered to be intersected, that is, the same point. To obtain as many intersection points as possible, all GEDI data for April to October 2019 that are currently open access were downloaded. The GEDI DEM corresponding to these intersection points was used to verify the average height values of ground photons. In total, 61 intersection points were extracted ( Figure 3).

BP-ANN Modeling, Extrapolation and Validation
The Backpropagation-Artificial Neural Network (BP-ANN) was used to model the relationship between canopy height and vegetation indices derived from Landsat data and to predict canopy height. More than 30 types of neural networks, such as black-box algorithm models, have been developed since the first prototype was proposed [60,61]. This model has been widely used in various fields due to its powerful prediction and forecasting capabilities [60,[62][63][64].
For model building, a three-layer BP neural network model with 10 input layer neurons, 11 hidden layer neurons and one output layer was constructed. The ten vegetation indices were used as

BP-ANN Modeling, Extrapolation and Validation
The Backpropagation-Artificial Neural Network (BP-ANN) was used to model the relationship between canopy height and vegetation indices derived from Landsat data and to predict canopy height. More than 30 types of neural networks, such as black-box algorithm models, have been developed since the first prototype was proposed [60,61]. This model has been widely used in various fields due to its powerful prediction and forecasting capabilities [60,[62][63][64].
For model building, a three-layer BP neural network model with 10 input layer neurons, 11 hidden layer neurons and one output layer was constructed. The ten vegetation indices were used as input layers. In total, 80% of the model training samples were randomly selected from the discontinuous CHM dataset for a total of 983 samples. The "scikit-learn" package in python was used to implement BP-ANN model training. A wall-to-wall forest canopy height map was derived by extrapolating the BP-ANN model to the whole study area with ten vegetation index bands.
For accuracy evaluation, the mapping accuracy of canopy heights was independently validated using three validation datasets that consist of the remaining 20% CHM dataset, a total of 66 field plots, and their combination. The root mean square error (RMSE) was calculated to quantify the error of the estimated forest canopy: where n represents the number of samples in the validation dataset, H i represents the actual canopy height value and true value, and the H i represents the predicted canopy height values. Figure 7 shows the photon classification results of some ATL08 samples. The ground photons shown in Figure 7 and the DSM are used to calculate the discontinuous CHM. Figure 8a (Table 4). This feature shows that the CHM dataset is representative for subsequent modeling and canopy height mapping.    (Table 4). This feature shows that the CHM dataset is representative for subsequent modeling and canopy height mapping.

Comparison of DSM and Ground Photon Values with GEDI Data
The 61 points shown in Figure 9a are the points at which the newly calculated average elevation of ground photons intersects with the GEDI footprint for the whole study area. The validation result for ground photon elevation with GEDI terrain elevation shows a coefficient of determination R 2 = 0.995 and an RMSE = 5.73 m. Even though the RMSE is high, which is explained further blow, more than 50% of the points in Figure 9a have the difference within ± 3 m between ground photon elevation and GEDI terrain elevation. The validation results of the ZY-3 DSM with GEDI land surface elevation show an R 2 = 0.991 and an RMSE = 6.59 m. The validation results for ground photon elevation with the SRTM DEM shows an R 2 = 0.997 and an RMSE = 6.76 m (Table 5).

Comparison of DSM and Ground Photon Values with GEDI Data
The 61 points shown in Figure 9a are the points at which the newly calculated average elevation of ground photons intersects with the GEDI footprint for the whole study area. The validation result  different from the capacity for ZY-3 stereo images to capture the top of vegetation and which is caused by the difference in sensor imaging. (c) The different sensors have varied effects. In this study, the 30 m × 30 m resolution ground photon elevation is calculated, and the spot diameter of the GEDI is exactly 30 m. (d) Terrain factors will significantly affect the reliability of waveform data. When data conditions permit, a more reliable verification result will be obtained by applying airborne Lidar in the study area.  Errors of the three precision verifications mentioned above mainly result from the following: (a) There is a deviation in the data acquisition time, as the GEDI only has data for 2019. (b) Figure 9b shows that the scatter fitting line is below 1:1 because the highest point of the first return Gaussian fitting peak of GEDI full waveform data represents the top of vegetation, which is significantly different from the capacity for ZY-3 stereo images to capture the top of vegetation and which is caused by the difference in sensor imaging. (c) The different sensors have varied effects. In this study, the 30 m × 30 m resolution ground photon elevation is calculated, and the spot diameter of the GEDI is exactly 30 m. (d) Terrain factors will significantly affect the reliability of waveform data. When data conditions permit, a more reliable verification result will be obtained by applying airborne Lidar in the study area. Figure 10 shows the forest canopy height estimated by the BP-ANN with training samples of the discontinuous CHM dataset. The estimated forest canopy height range of 3 m to 34 m is similar to the field measured forest canopy height range of 6.2 m to 29.2 m. Across the study area, forest canopy heights of 11 to 14 m accounted for the largest proportion. The spatial distribution of the estimated forest canopy height was found to be related to the elevation of the study area. The height of the forest canopy in mountainous areas of elevations of greater than 300 m is more than 11 m as opposed to that in bare area. Due to the effects of human activities and urban expansion, the spatial distribution of forests in plain areas is relatively fragmented while the distribution of forests in mountainous areas is relatively continuous. Most of the latter are evergreen coniferous forests with a long cutting period and evergreen broad-leaved forests with a short cutting period, such as eucalyptus forests. In addition, the unique karst landforms in Guangxi, especially in Longan County in western Nanning, are mostly distributed across natural forests because they cannot be used as cultivated land. The distribution of forest canopy height is random and is not as uniform as that of artificial forest.   Figure 11 shows scatter plots of independent accuracy validation for the estimated forest canopy height using 20%, a total of 257 CHM datasets, 66 field measured plots, and a combination of them. The elevation of field measured plots varies from 90 m to 402 m, with slopes of 1 • to 36 • . The age of the forest ranges from 2 years to 56 years, which was estimated by local forestry experts in the field. The average forest canopy height of field measured plots is 13.32 m. The validation results shown in Figure 11a-c present an R 2 = 0.51 and an RMSE of 3.34 m for part of the discontinuous CHM dataset, an R 2 = 0.51 and an RMSE of 3.47 m based on field measured plots, an R 2 = 0.51 and an RMSE for 3.38 m based on the combination of both, respectively, demonstrating the effectiveness of the BP-ANN model ( Table 6).

Large Scale Forest Canopy Height Mapping
We identified a new means to estimate forest canopy height using a combination of ICESat-2 ATLAS data and stereo-photogrammetry. This method uses spaceborne Lidar ICESat-2 ATLAS data to replace airborne Lidar data and in turn obtain ground elevation information. Its low cost, unlimited data acquisition range, and global sampling range allow it to be the ability to be applied to a wider range. Given our successful application of this method to a study area in Nanning, it can be used over large scales, such as regional and global forest canopy height mapping. In addition to the stereo image pairs provided by the ZY-3 satellite, China's first submeter high-resolution optical transmission stereo mapping satellite GF-7 can provide high-precision stereo image pair data for superimposing ICESat-2 data to estimate global forest canopy height. Although GF-2 was not designed as a

Large Scale Forest Canopy Height Mapping
We identified a new means to estimate forest canopy height using a combination of ICESat-2 ATLAS data and stereo-photogrammetry. This method uses spaceborne Lidar ICESat-2 ATLAS data to replace airborne Lidar data and in turn obtain ground elevation information. Its low cost, unlimited data acquisition range, and global sampling range allow it to be the ability to be applied to a wider range. Given our successful application of this method to a study area in Nanning, it can be used over large scales, such as regional and global forest canopy height mapping. In addition to the stereo image pairs provided by the ZY-3 satellite, China's first submeter high-resolution optical transmission stereo mapping satellite GF-7 can provide high-precision stereo image pair data for superimposing ICESat-2 data to estimate global forest canopy height. Although GF-2 was not designed as a photogrammetric system, some studies have made use of the different orbit data from GF-2 in repeated observation areas to realize stereo observation and extract DSM data [65].
Additionally, other techniques such as radargrammetry can produce a DSM by using SAR images to replace the DSM derived from ZY-3 in this study [44,45]. In addition to ICESat-2, the GEDI can provide surface elevation information as a DEM [32,57,58].
In addition, to obtain mapping results with the same resolution as those of Landsat, 30 m × 30 m grid is used to process ground photons. If we want to obtain a lower resolution result, such as those of the Moderate Resolution Imaging Spectroradiometer (MODIS), this is feasible, but for higher resolution mapping, we should consider enough ground photons to increase the accuracy and reliability of the DEM. Further analysis and research are needed.

How to Filter Effective ATL08 Data
The canopy height values provided by ATL08 data can be used as samples to directly generate a wall-to-wall forest canopy height map. Since the resolution of ATL08 data is 100 m, ATL08 data are not suitable for vegetation mapping at a 30 m resolution as with Landsat data, and should be used to map at resolutions of greater than 100 m to eliminate errors caused by position deviations. The selection of effective ATL08 data is a key problem. These two fitting lines are highly dependent on the accuracy of photon classification, creating inevitable errors in the canopy height values of ATL08.
Our numerical analysis shows that, for ATL08 data of the study area, 6% of canopy height values are greater than 50 m, which can be considered invalid values. These invalid values may originate from the low precision surface and vegetation surface curves fitted by sparse signal photons, or even attribute values obtained by interpolation. For the remaining 94% of the data, even if we use the rest of the ATL08 data fields (such as clouds, snow, and urban areas) to exclude the influence of land surface types and clouds, the reliability of canopy height data is affected by many other factors. For instance, (1) irregular topography increases photon classification errors, and (2) when the photon signal of a surface or vegetation is weak and the number of point clouds is small, the fitting error of the surface curve will increase.
Many uncertain factors make it difficult to automatically filter effective data when data are used in regional mapping. Figure 12 shows a frequency statistics histogram of canopy and ground photons. We suggest that canopy and ground photon numbers be used as a basis for screening effective ATL08 data, as the stronger the photon signal, the smaller the error becomes. However, determining the photon number threshold with universal or adaptive rules is an issue worthy of future research.

How to Filter Effective ATL08 Data
The canopy height values provided by ATL08 data can be used as samples to directly generate a wall-to-wall forest canopy height map. Since the resolution of ATL08 data is 100 m, ATL08 data are not suitable for vegetation mapping at a 30 m resolution as with Landsat data, and should be used to map at resolutions of greater than 100 m to eliminate errors caused by position deviations. The selection of effective ATL08 data is a key problem. These two fitting lines are highly dependent on the accuracy of photon classification, creating inevitable errors in the canopy height values of ATL08.
Our numerical analysis shows that, for ATL08 data of the study area, 6% of canopy height values are greater than 50 m, which can be considered invalid values. These invalid values may originate from the low precision surface and vegetation surface curves fitted by sparse signal photons, or even attribute values obtained by interpolation. For the remaining 94% of the data, even if we use the rest of the ATL08 data fields (such as clouds, snow, and urban areas) to exclude the influence of land surface types and clouds, the reliability of canopy height data is affected by many other factors. For instance, (1) irregular topography increases photon classification errors, and (2) when the photon signal of a surface or vegetation is weak and the number of point clouds is small, the fitting error of the surface curve will increase.
Many uncertain factors make it difficult to automatically filter effective data when data are used in regional mapping. Figure 12 shows a frequency statistics histogram of canopy and ground photons. We suggest that canopy and ground photon numbers be used as a basis for screening effective ATL08 data, as the stronger the photon signal, the smaller the error becomes. However, determining the photon number threshold with universal or adaptive rules is an issue worthy of future research.

Conclusions
This study proposes a new means with which to estimate forest canopy height, which used a combination of ICESat-2 ATLAS data and ZY-3 stereo images to extract a discontinuous CHM dataset as training samples and extrapolated the BP-ANN model to the whole study area with ten vegetation index bands from Landsat 8 images. Ground photons of ATL08 and ATL03 were recalculated to

Conclusions
This study proposes a new means with which to estimate forest canopy height, which used a combination of ICESat-2 ATLAS data and ZY-3 stereo images to extract a discontinuous CHM dataset as training samples and extrapolated the BP-ANN model to the whole study area with ten vegetation index bands from Landsat 8 images. Ground photons of ATL08 and ATL03 were recalculated to obtain the average terrain height value for a 30 m × 30 m grid. A discontinuous CHM dataset was derived from the ZY-3 DSM by subtracting new ground photon heights. The accuracy of the ZY-3 DSM and ground photons' average values was evaluated by GEDI data. The validation results show an R 2 of greater than 0.991 and an RMSE of approximately 6 m, and the source of the observed error was discussed in detail. Regional forest canopy height mapping with a resolution of 30 m was executed based on the BP-ANN model using the CHM dataset as a training sample combined with vegetation indices from Landsat 8 data. The independent accuracy validation for the estimated forest canopy height shows an R 2 = 0.51 and an RMSE from 3.34 m to 3.47 m based on part of the CHM dataset and field measured plots.
This paper presents a very promising means to map forest canopy heights at regional and global scales. It uses spaceborne Lidar data instead of airborne Lidar data to provide terrain information at a low cost and with global coverage. In the future, other types of data, such as radar image generating DSM and GEDI providing ground information can adopt the method used in this work, which will be conducive to the realization of multiple data sources that complement each other and work together to map regional forest canopy heights. Methods of selecting effective ATL08 data to extrapolate a wall-to-wall forest canopy height should be the focus of future research work.