1. Introduction
The spatial structure of forested terrain is listed as an important indicator for monitoring carbon stocks by the International Union of Forest Research Organizations (IUFRO) [
1,
2]. Assessing ground topography in forested terrain is a prerequisite to accurately determine the forest spatial structure; therefore, high spatial resolution modeling is necessary to characterize forest ecosystems [
3,
4]. Most optical remote sensing systems can measure ground topography in forested terrain; however, they have poor measurement accuracy (elevation difference = 2.9 m to 4.9 m) [
5,
6]. Spaceborne [
7], airborne [
8,
9], and terrestrial [
10,
11] Light Detection And Ranging (LiDAR) systems have shown great potential for acquiring accurate topographic information in this field. Although airborne and terrestrial LiDAR can accurately quantify ground topography in forested terrain, these methods remain largely impractical at large spatial scales due to high data acquisition costs [
8,
12,
13,
14]. Spaceborne LiDAR is unique since it comes with low acquisition costs and provides a synoptic perspective of certain plot-level details from orbit [
15].
The Ice, Cloud, and land Elevation satellite-1 (ICESat-1) [
16], the Global Ecosystem Dynamics Investigation (GEDI) [
17], and the Ice, Cloud, and land Elevation satellite-2 (ICESat-2) [
18] are typical spaceborne LiDAR systems. ICESat-1 and GEDI carry large-footprint and waveform LiDAR systems. The Geoscience Laser Altimeter System (GLAS) instrument aboard ICESat-1 was launched in 2003 and decommissioned in 2009 [
19]. GLAS is the first spaceborne LiDAR instrument designed to make global observations. GLAS operates a single laser beam from a ~600 km orbit at 40 Hz and has a 70 m diameter footprint and a ~170 m sampling rate along track.
GLAS waveform data has been successfully used to estimate the vertical structure of forest terrain, including ground topography and canopy heights [
19,
20,
21]. Harding et al. [
20] found that the waveform was an accurate representation of the canopy height distribution within a GLAS footprint. Lefsky et al. [
21] observed that the models combing GLAS waveforms and Shuttle Radar Topography Mission (SRTM) could explain ~59%–68% of the variance in the field-measured forest canopy height (root mean squared error––RMSE = 4.85–12.66 m); however, sloped ground in forested terrain reduced the canopy height accuracy by using waveform data. Chen [
22] found that the ground topography in forested terrain was the critical factor affecting the accurate measurement of canopy height using waveform data, and with increasing forest terrain complexity, the accuracy of estimating forest canopy height decreased. Fang et al. [
23] found that in forested terrain with complex ground topography, the GLAS waveform was characterized by multiple energy peaks, in which the ground topography might be broadened and mixed, making the extraction of canopy height difficult. In order to quantify the influence of ground topography on canopy height estimation using GLAS waveform data, Lee et al. [
24] found that without slope correction, the canopy height could be overestimated by 3 m over a 15 degree slope. Removing the ground topography in forested terrain from large LiDAR footprint could improve the accuracy of canopy height estimates. Claudia et al. [
25] revealed that GLAS height estimates were accurate for areas with a slope up to 10 degrees, whereas the waveform results for areas above 15 degrees were problematic. Ten-to-fifteen degree slopes have been found to be a critical crossover point. The aforementioned studies demonstrated that it was feasible to extract ground topography in forested terrain and canopy height from spaceborne waveform data at stand level; however, the accuracy of canopy height estimation was largely determined by the ground topography, and extracting canopy height across a large LiDAR footprint using waveform data over hilly or mountainous regions is a great challenge. The GEDI was launched on 5 December 2018; however, the GEDI spaceborne data has just recently been released, and no related study was found [
26].
The Advanced Topographic Laser Altimeter System (ATLAS) instrument aboard ICESat-2 was launched on 15 September 2018, and data was released on 30 May 2019. ATLAS is the first spaceborne photon-counting LiDAR instrument designed for continuous global observation of Earth [
27,
28,
29]. Different from the GLAS waveform-digitizing LiDAR system, ATLAS only responds to the presence of return signals and records the time tags with an output of 0 or 1; however, it does not record the return waveform [
30,
31,
32]. ATLAS operates six laser beams from a ~600 km orbit at 10 kHz and has a footprint (17 m in diameter) sampling rate of ~0.7 m along-track [
33,
34]. The center-to-center spacing along a track for ATLAS is narrower than that of GLAS (170 m). The high repetition rate enables ATLAS to obtain nearly continuous tracking information, which is necessary to measure the ground topography in forested terrain. While the GLAS LiDAR system uses a laser beam, the ATLAS configuration uses a diffractive optical element to split the laser into six beams arranged as three beam pairs, each of which consists of a strong and weak energy beam at a 4:1 ratio, allowing for local slope determination between each beam pair as well as compensation for varying surface reflectance [
27,
33,
34]. The travel time of each detected photon is used to determine a unique XYZ location on the Earth’s surface [
35,
36]. After ATLAS data was released, Neuenschwander found good correlations between matching Digital Terrain Model (DTM) from airborne LiDAR data and ATLAS data (R
2 = 0.99, RMSE = 0.85 m) [
37]. Wang et al. found that the overall mean difference and RMSE values between the ground elevations retrieved from the ICESat-2 data and the airborne LiDAR-derived ground elevations are −0.61 m and 1.96 m, respectively [
38]. However, he primarily examined the retrieved canopy height accuracy from the ICESat-2 strong beam and did not analyze the accuracy of the ICESat-2 weak beam. Under the same orbital conditions, ATLAS can acquire more continuous photon cloud data using the six-beam instrument with different laser pointing angles and laser intensities. The measurement accuracy of the different ATLAS channels remains to be quantified [
39]. To the authors’ knowledge, only a few studies have been carried out to analyze the multi-beam geometrical features for measuring ground topography in forested terrain from photon-counting data onboard ICESat-2. Therefore, the effective quantifying of the ground topography in forested terrain using the six-beam photon-counting data is essential to quantify the performance of the unique photon-counting instrument onboard ICESat-2.
The objective of this study is to assess the performance of the ICESat-2/ATLAS multi-channel photon data for estimating ground topography in forested terrain by comparing the derived ground topography from different ATLAS beam photon-counting data with that from Goddard’s LiDAR, Hyperspectral and Thermal imager (G-LiHT) data. The paper also analyzes the correlation between laser intensity parameters, laser pointing angle parameters, and estimated ground topography error in forested terrain.
2. Materials and Methods
2.1. Study Area
The study area (33.564°N, 81°684′W) is a forested area within the City of Aiken, South Carolina, USA (
Figure 1). Vegetation footprint types in the study area include cultivated land (0.04%), forest (88.52%), shrubland (0.51%), wetland (6.97%), and artificial surfaces (3.95%) [
39]. The upland forest has many tree species, including sand post oak (
Quercus margaretta), loblolly pine, water oak (
Quercus nigra), hickory (
Carya), and turkey oak (
Quercus laevis) [
40]. The elevation of the study area ranges from 91 m to 164 m. Vegetation coverage in the study area ranges from 25% to 66%.
2.2. Data
The ICESat-2 mission produces along-track ground topography in forested terrain that includes telemetry data (ATL00), reformatted telemetry (ATL01), science unit converted telemetry (ATL02), global geolocated photon data (ATL03), land vegetation elevation (ATL08), and a land/canopy grid (ATL18) [
37,
38,
39,
40,
41]. ATL00, ATL01, and ATL02 are original photon data sets without scientific algorithms. ATL03 is the geolocated photon cloud and serves as the input data for each of the higher-level data products such as ATL08 and ATL18. The ATL08 algorithm was developed specifically for the extraction of terrain and canopy heights from the ATL03 photon cloud data, and the ATL08 geophysical data product has a 100 m step size in the along-track direction [
33]. The ATL03 product not only includes latitude, longitude, height, and signal photon confidence level of each received photon, but also includes tx_pulse_energy, tx_pulse_skew_est, tx_pulse_width_lower, and tx_pulse_width_upper parameters, which may be related to laser intensity and laser pointing angle [
37]. All ICESat-2 data products were acquired from
https://search.earthdata.nasa.gov.
Here the ATL03 product parameters were used, including lat_ph, lon_ph, h_ph, geoid, delta_time, signal_conf_ph, sc_orient, tx_pulse_energy, tx_pulse_skew_est, tx_pulse_width_lower, and tx_pulse_width_upper. The names and corresponding descriptions of product parameters are listed in
Table 1 [
37,
38,
39]. In order to reduce the influence of noise photons in forested terrain ground topography measurement, a signal_conf_ph of 4 was used as the signal photon parameter evaluation standard.
A diagram of ICESat-2 for estimating ground topography in forested terrain is illustrated in
Figure 2. The forward orientation (sc_orient=1) corresponds to ATLAS traveling along the +x direction in the ATLAS instrument reference frame [
41,
42,
43]. The ATLAS signal photon shown in the yellow square represents the photons detected from the gt3r laser channel. The ATLAS footprint shown in the red square represents the photons detected from the gt3l laser channel. The number of photons in the gt3l channel is less than in the gt3r channel, which is due to the backward orientation of ATLAS [
43]. The photon level in the along-track direction was selected to calculate the ground topography in forested terrain. In the right figures, the two laser beams are 90 m apart. The ground topography measured by the two laser beams is similar, and the terrain elevation ranges from 130 m–145 m.
Due to the influence of sunlight as well as atmospheric and system noise, a large number of noise photons are present in the ATLAS data, which seriously reduce the ground elevation measurement accuracy. In order to improve the estimation accuracy of ATLAS photon data, NASA proposed a Differential, Regressive, and Gaussian Adaptive Nearest Neighbor (DRAGANN) method and ATL08 data classified algorithm to filter out noise photon data and classify ground photons [
41,
42,
43]. In order to explore the estimation accuracy of forested terrain from ATLAS data, this contribution chose to associate the ATL08 classified label with the ATL03 photon data and used the ground signal photons flag mentioned in ATL08 as ground photons to establish an ATLAS-based DTM (
Table 2) [
37].
To assess the accuracy of the six beam-ATLAS DTM, the DTM obtained from the ATLAS data was compared with airborne discrete-return LiDAR data, collected for the same longitude and latitude using the multi-sensor instrument G-LiHT [
44]. G-LiHT provides distributed laser pulses for measuring ground topography and canopy heights (
Table 3) [
45].
Both Level 2 and Level 3 products along the flight transects were generated from airborne LiDAR data from the G-LiHT science team. The DTM has a 1 m-resolution and was released as a Tag Image File Format (TIFF) profile. The DTM was assessed to validate the ground topography accuracy [
44,
45]. The trajectory of the G-LiHT KML (Keyhole Markup Language) data (blue line) and ATLAS data (green line) illustrates the location of the study area of the NASA EARTHDATA (
Figure 3). This illustration also includes two G-LiHT DTM profiles used in the study, AMIGACarb_Augusta_FIA_Sep2011_l16s597_DTM.tif and AMIGACarb_Augusta_FIA_ Sep2011_ l40s557_DTM.tif, respectively.
2.3. Methodology
The primary challenge of this contribution centers on reducing the influence of noise photons on the ground elevation data derived from ICESat-2 data, distinguishing canopy signal photons and ground signal photons, and matching the ATL03, ATL08, and G-LIHT data. Although ATLAS data has more noise photons, the NASA official team used the DRAGANN and an algorithm for determining Land Vegetation along-track to provide classification labels (classed_pc_flag) for the photon data [
43]. In this contribution, the ground signal photon classification label ATL08 is used for ground photons, and the DTM
ATLAS will be established based on ATL03 data and ATL08 label. This contribution presents a quantitative assessment of the ground topography in forested terrain using ATL03 and ATL08 data compared to airborne G-LiHT LiDAR data.
The geolocation between the ATLAS and G-LiHT data is not completely along orbit; therefore, this paper proposes an approach based on the ATL03 profile to match these two datasets (
Figure 3). To clearly illustrate the proposed methodology, an overview of the major steps is exhibited in
Figure 4 and described as follows:
(1) Identifying study site. In this step, we combine ATL03 lat_ph, lon_ph, and G-LiHT KML profile and identify the study site.
(2) Obtain parameters from the ATL03 and ATL08 data by matching different channels under the same orbit conditions using time tags (delta_time in ATL03 HDF5 profile). To extract the photon’s height (h_ph which is relative to the WGS-84 ellipsoid), latitude (lat_ph) and longitude (lon_ph), signal_conf_ph, sc_orient, tx_pulse_energy, tx_pulse_skew_est, tx_pulse_width_lower, and tx_pulse_width_upper parameters from the ATL03 HDF5 profile, combine the geoid and h_ph by interpolating. Extract the photon classification parameters (classed_pc_flag), classed_pc_indx, and ph_segment_id parameters from the ATL08 HDF5 profile.
(3) Establishing the relationship between ATL03 and ATL08 data photon classification parameters by classed_pc_indx, ph_segment_id and applying each photon classification label from ATL08 to each photon data from ATL03.
(4) Establishing the DTMATLAS. The photons with a signal confidence flag from high confidence (signal_confidence = 4 in ATL03 HDF5 profile) and photon classification parameter (classed_pc_flag=1 in ATL08 HDF5 profile) were used to establish the DTMATL03.
(5) Obtaining the DTMG-LiHT. In this step, we extract the latitude-longitude information from DTMATL03 to match the DTMG-LiHT profile corresponding position generated from G-LiHT and generated DTMG-LiHT with the corresponding ATLAS footprint latitude-longitude. If the absolute difference between the elevation of ATLAS ground photons and the corresponding elevation of the DTMG-LiHT is more than 20 m, this photon was classified as an invalid ground photon.
(6) Assessing the performance of ICESat-2/ATLAS multiple channels photon data. In this final step, we compare the DTMATL03 profile with the corresponding DTMG-LiHT profile, compute and analyze the evaluating indicator from different channels. In order to quantify the influence of different laser intensity parameters and laser pointing angle parameters on the estimation accuracy of ground elevation, corresponding four parameters as follow: tx_pulse_energy, tx_pulse_skew_est, tx_pulse_width_lower and tx_pulse_width_upper are extracted and analyzed the relationship between the four parameters and elevation errors.
The
Figure 5 shows the DTM files of G-LiHT and ATLAS photon data corresponding to the two tracks in the study area. The green dots are ATLAS footprints.
2.4. Accuracy Evaluation
DTM data derived from airborne G-LiHT LiDAR data were utilized to assess the accuracy of the ATLAS-derived ground elevations. The ground elevation errors in the ATLAS data mentioned in this contribution are calculated by subtracting the G-LiHT elevations from the ATLAS elevations [
38]. Seven statistical variables, including root mean squared error (RMSE), mean absolute error (MAE), coefficient of determination (R
2), mean error (ME), Pearson correlation coefficient, Spearman correlation coefficient, and Kendall correlation coefficient between the ground elevations and the corresponding G-LiHT’s DTM values were calculated to quantitatively evaluate the accuracy of the ATLAS-derived ground elevations. Three statistical variables, Pearson correlation coefficient, Spearman correlation coefficient, and Kendall correlation coefficient between the elevation errors and the corresponding laser intensity parameter and laser pointing angle parameters were calculated to quantitatively evaluate the correlation of laser intensities and laser pointing angles with the elevation error in estimating ground topography in forested terrain.
4. Discussion
In order to study the influence of laser intensity and laser pointing angle on ground elevation estimation accuracy, the ATLAS data errors under different laser intensities and laser pointing angles are analyzed respectively, and the correlation between elevation errors and corresponding laser parameters is examined.
4.1. Retrieved Ground Topography in Forested Terrain for Different Laser Intensities
The mean statistical indicators of the different laser intensities are listed in
Table 5 and the correlation coefficient statistics between tx_pulse_energy parameters and the ATLAS data elevation errors are listed in
Table 6.
The estimated ground topography for different intensity beams shows significant agreement with the reference DTM elevations. For the weak beam and strong beam, the mean R2 values are 1, the ME values are 0.33 m and 0.27 m, and the RMSE values are 0.76 m and 0.74 m, respectively. The correlation coefficient of all types is greater than 0.89.
For the varying laser intensities in this data set, the statistical indicators for the strong beam performed better than that of weak beam with the lower RMSE (RMSEstrong beam = 0.74 m and RMSEweak beam = 0.76 m), lower MAE (MAEstrong beam = 0.51 m and MAEweak beam = 0.54 m), lower ME (ME strong beam = 0.33 m and ME weak beam = 0.27 m), and higher correlation coefficient. A possible reason is that the weak beam channel has fewer signal photons compared to the strong channel, making measuring ground topography in forested terrain using the weak beam more difficult. Using a strong beam, ATLAS could produce more signal photons than under the weak beam. The laser intensity ratio of strong beam to weak beam is 4:1. Depending upon the surface reflectance and atmospheric conditions, up to 16 photons per outgoing shot could be detected for the strong beam, while the weak beam could detect only 4 photons. However, the strong and weak beams can both provide useable data for measuring ground topography in forested terrain.
To further explore the influence of laser intensities on elevation errors, we calculated three correlation coefficients between tx_pulse_energy parameters and elevation errors (
Table 6). For all the data, the correlation coefficients for the elevation errors and tx_pulse_energy parameters are greater than 0.5. In addition, the Spearman correlation coefficient values for the various laser intensities are greater than 0.74, indicating there is a significant correlation between the tx_pulse_energy parameters and elevation error. However, this contribution only explores the correlation between tx_pulse_energy and elevation error for a laser intensity ranging from 0.02 mJ to 0.09 mJ. Future studies with the tx_pulse_energy parameters will need to perform a more detailed analysis on the effect of tx_pulse_energy on elevation error.
Compared to the forested terrain ground topography estimation method in proposed by Neuenschwander et al. [
37], higher R
2 values and lower RMSE values were observed in the strong beam mode and weak beam mode. However, this study proposes a photon level which is different to the Neuenschwander et al. method [
37], which notes that the result of a photon subset using a strong beam and a weak beam can reasonably explain the ground topography in the forest study area.
4.2. Retrieved Ground Topography in Forested Terrain Elevation with Different Laser Pointing Angles
The mean statistical indicators for the different laser pointing angles are listed in
Table 7 and the correlation coefficient statistics between tx_pulse_skew_est, tx_pulse_width_lower, tx_pulse_width_upper parameters and the elevation errors of ATLAS data are listed in
Table 8.
The estimated ground topography in forested terrain using different laser pointing angles shows strong agreement with the reference DTMG-LiHT elevations, as demonstrated by the R2 values equaling 1.00 and the RMSE values less than 0.92 m.
Results indicated that the gt1 channel pointing (R
2gt1 channel = 1.00, RMSE
gt1 channel=0.62 m, MAE
gt1 channel = 0.4m) performed better than gt2 channel pointing (R
2gt2 channel = 1.00, RMSE
gt2 channel = 0.92 m, MAE
gt2 channel = 0.62m) and gt3 channel pointing (R
2gt3 channel = 1.00, RMSE
gt3 channel = 0.71 m, MAE
gt3 channel = 0.56m), which was due to several hardware reasons. On one hand, the gt1 channel could achieve more effective forested terrain signal photons than other channels in the study area. More signal photons can give a clearer depiction of ground topography in forested terrain. On the other hand, the photon rates of the gt1 and gt3 channels are higher than the gt2 channel, which is consistent with the description of the different laser pointing angles [
46].
To further explore the influence of laser pointing angles on elevation error, we calculated three correlation coefficients between tx_pulse_skew_est, tx_pulse_width_lower, tx_pulse_width_upper parameters and elevation error. In the ATL03 Algorithm Theoretical Basis Document (ATBD), these parameters may be related to the laser pointing angles. The quantitative results of the correlation coefficients are summarized in
Table 8. For all the data, the mean correlation coefficients for the elevation errors are less than 0.20. There is no significant correlation between the tx_pulse_skew_est, tx_pulse_width_lower, tx_pulse_width_upper parameters and the elevation error. However, this contribution only explores the correlation between tx_pulse_skew_est, tx_pulse_width_lower, tx_pulse_width_upper and the elevation error. Future studies needed to analyze other laser pointing angles parameters’ relative elevation errors.
The results of this contribution performed better than that proposed by Neuenschwander et al. [
37] (R
2 = 0.99, RMSE = 0.85), which notes that the result of a photon subset using different laser pointing angles can reasonably explain the ground topography in the study area.
4.3. Retrieved Ground Topography in Forested Terrain Elevation with ATLAS
Most optical remote sensing systems could provide images of the horizontal distribution of ground topography, and the product generally follows the uppermost surface elevation (i.e., representing a digital surface model, DSM). However, the optical remote sensing images do not provide detailed information on the vertical distribution of ground topography in forested terrain, without regard to whether the surface is comprised of forest or not [
5,
6,
47]. In contrast, the LiDAR photon counting signature from ICESat-2/ATLAS could provide a direct depiction for ground topography in forested terrain. In this contribution, the close correspondence between the ATLAS and G-LiHT (R
2 = 1.00, RMSE = 0.75 m) confirms that the received photon data are an accurate representation of the ground topography forested terrain elevation within the ATLAS footprints. As a comparison, in areas of low relief (slope ≤ 5°) and middle dense tree cover (tree cover = 20%–40%), the mean and standard deviation of elevation differences between the ICESat/GLAS centroid and SRTM are −2.48 ± 4.04 m [
48]. Thus, the estimation of ground topography for forest-covered areas is able to be accomplished with ICESat-2/ATLAS.
The results from this contribution indicate that the ground topography in forested terrain elevation can be estimated using photon data from ICESat-2/ATLAS multi-channel. We were able to retrieve terrain elevation successfully in forest-covered areas. Prior work showing the correlation between spaceborne LiDAR-measured canopy height and ground topography in forested terrain [
21,
22,
23], which provide confidence that ICESat-2/ATLAS photon data in combination with GEDI data can substantially contribute to a global inventory of forest biomass. The work also provides insights for future work to improve the accuracy of the canopy height estimations.
5. Conclusions
In this contribution, ICESat-2 data is used to measure ground topography in forested terrain using different channels. The retrieved ground topography was validated by experiments with G-LiHT airborne data at different laser pointing angles and laser intensity types on the same route. Based on the results, the following conclusions can be drawn:
(1) Both qualitative and quantitative results indicate that at all laser intensities and laser pointing types resulted in a mean R2 = 1.00 and mean RMSE = 0.75 m, highlighting the ability of the ATL03 and ATL08 data to retrieve ground elevations.
(2) A significant correlation exists between the tx_pulse_energy parameters and elevation error. There is no significant correlation between tx_pulse_skew_est, tx_pulse_width_lower, tx_pulse_width_upper parameters and elevation error.
These conclusions give valuable insight into the ground topography in forested terrain using different ATLAS channels. Nevertheless, there are still many issues to be addressed in the future. Since ATLAS data is still in the research stage, we only considered the effects of laser pointing angles and laser intensity on retrieving ground topography. Other factors (e.g., canopy height, canopy cover, etc.) influencing the results were not considered. Therefore, the effects of other factors on retrieving ground topography over forested terrain using ATLAS data should be thoroughly examined in the future.