ICESat-2 ALTIMETRY AS GEODETIC CONTROL

Digital elevation models (DEMs) are of fundamental importance for a large variety of scientific and commercial applications. Many geoscience studies require the most precise and current information about the Earth’s topography. Independent quality assessments of these DEMs are crucial to their appropriate use in land process studies, as inputs to models, and for detection of topographic change. The Ice, Cloud and land Elevation Satellite (ICESat) provided globally-distributed elevation data of high accuracy that demonstrated to be well-suited for evaluating continental DEMs after appropriate editing (Carabajal and Harding, 2005; Carabajal and Harding, 2006; Carabajal et al., 2010 and 2011; Carabajal and Boy, 2016). ICESat-2, launched on September 15th, 2018, provides an opportunity to develop a dataset suitable for Geodetic Ground Control. With increased coverage, ICESat-2/ATLAS features 6 laser beams with 532 nm wavelength, using photon counting technologies. With a nearly polar orbit, altimetry from ICESat-2 is available for latitudes reaching up to 88 degrees, on a 91-day repeat orbit with monthly sub-cycles. ICESat-2’s footprint size is ~17m, at 10 kHz pulse repetition frequency, or 0.75 m along track. Its pointing control is 45 m, with a pointing knowledge of 6.5 m, and a single photon precision of 800 ps. Sophisticated data processing techniques on the ground, optimized by surface type, produce high quality estimates of topography. We illustrate the use of ICESat-2 altimetry to assess DEM’s accuracy using ATL08 release 002 elevations (Land and Vegetation) products (Neuenschwander and Pitts, 2019), showing comparable results to those using ICESat-derived Geodetic Ground Control.


ICESat-2 Mission and the ATLAS Instrument
The Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) mission (Markus et al., 2017;Neuman et al., 2019b) is the successor to the ICESat mission (Zwally et al., 2002;Schutz et al., 2005). ICESat-2 was launched on September 15th, 2018, into a 92° inclination, 500-km altitude, near-circular orbit. The mission aims primarily to monitor changes in the cryosphere, quantifying the contributions to sea-level change from glaciers, and ice sheets, and processes driving them, characterizing annual changes in thickness of sea ice to examine ice/ocean/atmosphere exchanges of energy, mass and moisture, and collecting valuable data to feed predictive models. In addition, ICESat-2 is collecting valuable scientific data globally, quantitatively characterizing topography and vegetation to measure vegetation canopy height as a basis for estimating large-scale biomass and biomass change, monitoring inland water, sea level changes, densely sampling the Earth's surface and the atmosphere. The mission's scientific objectives require tight vertical and horizontal accuracy of the measurements.
The Advanced Topographic Laser Altimeter System (ATLAS) onboard ICESat-2 is capable of detecting single photon level reflections (Neumann et al., 2019b;. ATLAS has a repetition rate of 10 kHz, using a 532 nm wavelength laser, with 6-beams, separated by ~2.5 to 3.5 km on the surface, grouped into 3 "tracks" of strong/weak beam pairs. This beam configuration allows for measurement of the surface slope in both the along-and across-track directions. The beam energy ratio (strong to weak) is 41, while the beam energy per pulse is 175 J  17 J for the strong and 45 J  5 J for the weak spots. The pulse width is about 1.5 ns FWHM (full width half max), with a spot diameter on the Earth surface of about 17 m. The receiver field-of-view diameter is ~45 m. The instrument precisely records the roundtrip travel time of each returned photon, determining times of flight (TOF) with a precision of 800 ps. The instrument is sensitive to solar background noise, making it a challenge to distinguish the surface photons (signal) from those scattered within the atmosphere (noise). Solar background photon rates recorded by ATLAS vary mainly with the sun angle, also with the atmospheric and Earth reflectivity at 532 nm at the specific location, reaching ~10 MHz in regions with high solar angle and reflectance in the ATLAS field of view under clear skies. Markus et al. (2017), Neumann et al. (2019b) and  discuss the ICESat-2 mission and ATLAS instrument characteristics in detail.

ICESat-2 Data Products and DEM Used
For this investigation, we have worked with the ICESat-2 ATL08 (ATLAS/ICESat-2 L3A Land and Vegetation Height) products . Their distribution is global, and also include ice sheets. The ATL08 products process as input medium-high confidence classified photons from the ATL03 (ATLAS/ICESat-2 L2A Global Geolocated Photon Data) products (Neumann et al., 2019a), as identified by the signal_conf_ph flags with values of 3 and 4. The photon cloud is then processed using the DRAGANN processing scheme, described in the ATL08 Algorithm Theoretical Basis Document (ATBD) . Output products include canopy and terrain metrics for a fixed segment distance of 100 m along track. Data included in this study was collected between. October 14th, 2018 and November 15th, 2019.
We used the Release-3 of the void-filled 3-arcsecond Shuttle Radar Topography Mission (SRTM) DEM (SRTMGL3 v003) provided by the National Aeronautics and Space Administration (NASA) and the National Geospatial-Intelligence Agency (NGA) (https://lpdaac.usgs.gov/products/srtmgl3v003/, https://doi.org/10.5066/F7F76B1X) to evaluate the differences with respect to the ATL08 elevation products. We selected this product because of its spatial resolution (~90 m), which is equivalent to the resolution of the ATL08 land products.

DEM Accuracy Assessments with Laser Altimetry from Space
ICESat altimetry has been extensively use to evaluate the quality of global elevation models Harding, 2005 and2006). A dataset was developed for topographic Ground Control Points (GCPs), for which a very stringent editing criteria was used to identify high quality global GCPs for land applications. For its development, various ways to eliminate the possibility of including data contaminated with effects that result in degradation in the quality of the elevation products, such as return saturation, and cloud contamination, and data from low energy returns were explored. The GCPs from ICESat have subdecimeter vertical accuracy and better than 10 m horizontal accuracy. The development of this database is documented in Carabajal et al. (2010 and, and the data include laser returns that satisfy high accuracy elevation requirements. Their accuracy estimates have been supported by accuracy estimates from rigorous analysis of instrument calibration and validation schemes using ocean scan maneuvers and cross-overs (Carabajal et al., 2011). This high quality laser altimetry estimates of land elevations have been used in many validation studies to evaluate the quality of Digital Elevation Models (DEMs) like those produced by the SRTM mission (Farr and Kobrick, 2000), described in Harding (2005 and2006) and Carabajal et al. (2010), evaluations of GMTED2010  described in Carabajal et al. (2011), and as part of the validation efforts for various versions of ASTER GDEM, shown for ASTER GDEM V3 (NASA/METI/AIST/Japan Spacesystems, and U.S./Japan ASTER Science Team, 2019) as in Carabajal and Boy (2016). These studies have looked at the distribution of elevation data quality as a function of terrain elevation and relief, roughness, and vegetation cover. The ICESat-2 elevation products present us with the opportunity to develop a similar high-quality dataset that can be used for these types of studies, after identifying high quality data suitable for this purpose.

ICESat-2 Elevation Uncertainty
ATLAS accuracy is a composite of ranging precision of the instrument, radial orbital uncertainty, geolocation knowledge, forward scattering in the atmosphere, and tropospheric path delay uncertainty. The ranging precision is a function of the laser pulse width, the surface area potentially illuminated by the laser, and the uncertainty in the timing electronics. The requirement on radial orbital uncertainty is specified to be less than 4 cm and tropospheric path delay uncertainty is estimated to be 3 cm. The ranging precision for flat surfaces, is expected to have a standard deviation of approximately 25 cm. The composite of each of the errors can also be thought of as the spread of photons about a surface, the point spread function or Znoise. The estimates of radial orbital uncertainty, geolocation knowledge, forward scattering in the atmosphere, and tropospheric path delay uncertainty for a photon will be represented on the ATL03 data product as the final geolocated accuracy in the X, Y, and Z (or height) direction ( Z). These parameters have different temporal and spatial scales, and vary over time. In the ATL03 products, Z represents the best uncertainty for each geolocated photon, but it does not incorporate the uncertainty associated with local slope of the topography. The slope component to the geolocation uncertainty is a function of both, the geolocation knowledge of the pointing (which is required to be less than 6.5 m) multiplied by the tangent of the surface slope. For less than 1-degree slope (flat topography) this uncertainty is  0.25 m. For a 10-degree surface slope, it can reach 1.19 m. When combined with Z, the uncertainty associated with the local slope will produce the sigma_atlas_land. Ultimately, the uncertainty reported on ATL08 includes the sigma_atlas_land, and the local rms values of heights computed within each data parameter 100 m segment.
The uncertainty of the terrain height for a 100 m segment in the elevations reported in the ATL08 products is described in Equation 1.4 of the ATL08 Algorithm Theoretical Document (ATBD) . The uncertainty of the mean terrain height for the segment is given by the h_te_uncertainty parameter. It incorporates all systematic uncertainties (e.g. timing, orbits, geolocation, etc.) as well as uncertainty from errors of identified photons. This parameter is described in Section 1, Equation 1.4. When there are not a sufficient number of ground photons in the point cloud classified as canopy or ground in the ATL08 processing, an invalid value is reported and no interpolation will be done to compute an elevation. The parameter h_te_std is the standard deviations of terrain points about the interpolated ground surface within the segment, and provides an indication of surface roughness.
The pre-launch best estimate of the ICESat-2's expected horizontal accuracy is 4.9 m at 1-sigma, while the single photon horizontal geolocation requirement is 6.5 m at 1-sigma. On-orbit estimates, needed to understand the actual performance using several months of post-launch calibration and validation efforts, have resulted in the current estimates of accuracy. They consider current estimates of ranging, timing, positioning and pointing. The Precision Orbit Determination (POD) team has performed pointing calibration solutions using data from all planned roundthe-world scan calibration maneuvers (Scott Luthcke, personal communication). Calibration pointing biases include timevarying mean roll/pitch biases, roll/pitch bias orbital variation, including any necessary reference frame corrections to the data when necessary.
At the time of this analysis, the POD team's estimated orbit performance exceeded 3 cm radial RMS accuracy. High elevation independent SLR residuals indicate 1.65 cm RMS radial orbit accuracy. Long-wavelength (~1700 km) pointing bias (time varying orbital variation and bias) 1-sigma geolocation error contribution was ~1.7 m. Preliminary range bias calibrations showed trends of 0.29 mm/day over all spots, with trends ranging from 0.16 to 0.35 mm/day across spots. Roll and Pitch error after calibrations showed a spread 1 arcsecond or ~2 m on the ground; 1-sigma is 0.3 arcsecond or 0.6 m error on the ground. Pointing calibration of orbital variation and time varying bias has significantly improved geolocation to meet mission requirements. Roll pointing error is significantly subarcsecond, meeting mission requirements. These estimates will be largely improved with the inclusion of crossovers.

Editing Strategy
For this paper, we chose to compare ATL08 terrain elevations to the elevations in SRTM 90 m model. Statistics for all beams have been computed after appropriate editing was performed. In this section, we show statistics and plots only for the Strong beam 3, the closest to nadir, as a representative example.  Table 1. Statistics for elevation Differences between SRTM 90 m and ICESat-2 ATL08 h_te_interp (m) shown in Figure 1. Statistics have been computed before editing, and after editing the data for h_te_uncertainty  7.5 m.
Data identified as "water" using the segment_watermask was excluded when equal to 1 (0 = no water). This parameter is referenced from the Global Raster Water Mask at 250 m spatial resolution (Carroll et al, 2009; available online at http://glcf.umd.edu/data/watermask/). The ATL08 parameter h_te_uncertainty, the total uncertainty of ground height estimates which includes all known sources of uncertainty in geolocation, pointing angle, timing, radial orbit errors, etc. and the uncertainty due to local slope was used to edit the data before all comparisons were made. Only data with h_te_uncertainty less or equal to 7.5 m was included. About 3 % of the data is excluded using this total uncertainty threshold.
We looked at the differences between SRTM 90 m and ATL08 terrain with respect to the standard deviation of ground height (h_te_std) with respect to the interpolated ground surface (h_te_interp), which provides an indication of the surface roughness. We also examined them with respect to the terrain slope within each segment (terrain_slope), which is computed from a linear fit of the terrain photons. Their calculations are shown in . Before editing, means, standard deviations and RMSEs increase with increasing roughness and slope, as expected. For Australia, higher relief, and slope regions are generally correlated with regions that are more vegetated. Therefore, after editing for tree cover as shown in the subsequent sections, to only include relatively unvegetated terrain in these comparisons, the importance of showing the statistics of the differences with relief and slope became less significant, and did not include them in this analysis, or used them for further editing. However, for other geographic regions, looking at the relationship with these parameters may be important in defining and adequate editing criteria.
We examined the possibility of further editing the ICESat-2 data taking into consideration some of the parameters derived from the atmospheric products. We considered only keeping data from clear skies, indicated by a Cloud Confidence Flag (cloud_flag_asr) equal to 0, and when the Multiple Scattering Warning Flag (msw_flag) was equal to 0, indicating no presence of multiple scattering due to clouds. The information distributed by the NSIDC website regarding known issues with ICESat-2 products indicates that the cloud_flag_asr parameter has shown to work well over Antarctica and Greenland, and over the oceans, but it has issues over land, where it tends to underestimate cloud cover, making it a less reliable indicator of cloud contamination. We did not investigate the surface_h_dens parameter, very infrequently defined, indicating the possible presence of low clouds (< ~200 m) which could be possibly misidentified as the surface. Editing with these additional cloud parameters reduced the data available for this analysis to less than 10 % of all data with reasonable h_te_uncertainty values, and precluded us from having a sensible geographic distribution throughout the continent with which to perform the evaluations.

Differences with Respect to Landcover
We used the segment_landcover parameter, based on the 0.5 km MODIS land cover product from 2012 (Channan et al, 2015;Friedl et al, 2010, available at http://glcf.umd.edu/data/lc/index.shtml). After editing based on the h_te_uncertainty less or equal to 7.5 m, the most abundant land cover in Australia is represented by Open Shrubland (L-Class = 7), followed by Savannas, Woody Savannas, Grasslands and Croplands (L-Class = 9, 8, 10 and 12, respectively).
The statistics of the SRTM 90 m differences with h_te_interp and landcover are shown in Table 2 and Figure 2, which also shows the geographic distribution of land covers represented. Elevation differences for Open Shrublands are largely Gaussian, and have means of 2.28 m, with a median of 2.31 m, and standard deviation and RMSE of 2.54 m and 3.41 m, respectively (See Table 2). The Skewness and Kurtosis of the distribution is 0.17 and 52.53. The tails of the differences range from -113.43 m to 137.20 m.

Elevation Differences with Respect to Tree Cover
For the data selected with total h_te_uncertainty of less or equal to 7.5 m, we further identified selected data for relatively bare ground cover by using the "landsat_perc" parameter, which represents the average percentage value of the valid (value 100) Landsat Tree Cover Continuous Fields product for each 100 m segment. Statistics were computed for categories starting at 0% Tree cover, in 5% increments (Figure 3 and Table 3). The majority of the data falls within the less than 5% Tree cover, and less than 20% Tree cover, with a geographic distribution show on the top of Figure 3, mostly in the central part of the continent. Higher tree cover concentrations are distributed in the coastal regions of the South West and South East of the continent. Even though the number of points for categories with larger tree cover decreases. Those data represent terrain elevations computed from photon clouds that also include the signal from the canopy after classification is performed. Means and standard deviations for those categories increase with increased percent tree cover.  The SRTM 90m -ATL08 differences increase with Tree cover, illustrating that the SRTM phase center elevation is typically located within the canopy (Carabajal and Harding, 2006). The ATL08 terrain elevations are more comparable with SRTM elevations corresponding to low tree cover regions. Therefore, in the following sections, we will illustrate comparisons where the h_te_uncertainty less or equal to 7.5 m and the percent tree cover is less or equal to 5%. About 20 % of the data is excluded using this total uncertainty threshold for relatively low tree cover.

Elevation Differences with Respect to Signal to Noise
The Signal to Noise Ratio of geolocated photons (snr parameter) is determined by the ratio of the superset of ATL03 signal and DRAGANN found signal photons used for processing the ATL08 segments to the background photons (i.e., noise) within the same ATL08 segments. Table 4 shows the statistics when data with total h_te_uncertainty of less or equal to 7.5 m in combination with landcover less than 5 % Tree cover is used for editing. Figure 4 shows the geographic distribution of snr, and the elevation differences statistics with respect to SRTM 90 m to look at the differences for relatively bare ground cover. Table 4. Statistics for elevation differences SRTM90 -h_te_interp (in m) for h_te_uncertainty less or equal to 7.5 m and %Tree (%T) less or equal to 5% with respect to Signal to Noise Ratio (snr) corresponding to Figure 4. Bins incremented by 5. Figure 4. Differences between SRTM 90 m and ATL08 h_te_interp (in m) with respect to Signal to Noise Ratio (snr) for h_te_uncertainty less or equal 7.5 m and %Tree less than 5%. Mean, Median, Standard Deviations and RMSE are shown for snr using bin increments of 5.
For relatively bare cover, 20% of the data is being edited (going from ~14.4 million returns to ~ 10 million returns). The distributions do not include as many outliers. Means and Medians The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII- B3-2020, 2020XXIV ISPRS Congress (2020 are reduced by ~0.30 m ranging from 1.92 m to 2.23 m, and standard deviations and RMSEs are also reduced by ~0.70 m, ranging from 2.31 m to 2.91m and 3.10 m to 3.54 m, respectively.
There is no clear relationship between the elevation differences and the signal to noise for the 100 m segments after editing is performed.

Elevation Differences with Respect to Apparent Surface
Reflectance.
The asr parameter represents the apparent surface reflectance computed in the ATL09 atmospheric processing, reported as valid in the ATL08 segment when reported in the cloud products. For relatively bare ground cover, there is no clear relationship between the apparent surface reflectance and biases in SRTM. Figure 5. Differences between SRTM 90 m and ATL08 h_te_interp (in m) with respect to Apparent Surface Reflectance (asr) for h_te_uncertainty less or equal to 7.5 m and %Tree cover less or equal to 5%. Mean, Median, Standard Deviations and RMSE are shown for asr using 0.1 bin increments. Table 5. Statistics for elevation differences SRTM90 -h_te_interp (in m) for h_te_uncertainty less or equal to 7.5 m and %Tree cover less or equal to 5%, with respect to Apparent Surface Reflectance (asr) corresponding to Figure 5 above, using 0.1 bin increments.

Elevation Differences with Respect to Number of Ground Photons in the Segment
We looked at the differences with respect to the Number of Ground Photons identified in the 100 m segment (n_te_photons). After editing for the h_te_uncertainty less or equal to 7.5 m, the majority of the data includes terrain segments sampled with up to 500 photons or less, and mostly sampled with 100 to 200 photons (7,872,854). The means of the differences between SRTM 90 m and ATL08 terrain range from 1.59 m to 3.45 m. The smallest mean, when the number of photons per segment are between 400 and 500, is 1.59 m, with a 1.68 m median, and a standard deviation and RMSE of 3.43 m and 3.82 m, respectively. These regions are geographically sampling the interior of Australia ( Figure 6). Categories when the number of photons are more than 500 per segment as less represented, and generally show larger mean differences, ranging from 2.36 m to 3.29 m, although the medians are more stable. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) Table 7. Statistics for elevation differences SRTM90 -h_te_interp (in m) for h_te_uncertainty less or equal to 7.5 m and %Tree cover less or equal to 5% with respect to Number of Ground Photons per Segment (n_te_photons) corresponding to Figure 6 above, using 100 photons bin increments.
For low tree cover (Table 7 and Figure 6), the statistics show that the differences between SRTM 90 m and ATL08 terrain for up to 500 photons per segment have means that vary between 1.32 m to 2.20 m, the standard deviations and RMSEs decrease from 4.01 m and 4.38 m, to 2.78 and 3.09 m, respectively. The more Gaussian distributions are seen for segments sampled with between 100 and 400 photons. The mean differences seem to indicate that mean positive biases in SRTM 90 m are slightly less than 2 m for unvegetated terrain. Figure 6. Differences between SRTM 90 m and ATL08 h_te_interp (in m) with respect to Number of Ground Photons per Segment (n_te_photons) for h_te_uncertainty less or equal to 7.5 m and %Tree cover less or equal to 5%. Mean, Median, Standard Deviations and RMSE are shown using 100 photons bin increments.

Elevation Differences for Data from Strong and Weak Beams
We computed the statistics for the differences for the edited data for all 6 laser beams, and segregated the data based on day and night (Table 8). During daytime conditions, the ICESat-2 data is collected during larger solar background rate conditions, challenging the algorithms used to classify signal photons and an unbiased estimation of terrain elevations. The observations made by the Strong beams are more than twice more abundant than those obtained by the Weak beams. When only looking at the differences when editing only by the h_te_uncertainty of less than or equal to 7.5 m, there is an ~0.20 m difference between Strong beams (1, 3, 5) and Weak beams (2, 4, 6), implying that elevations sampled by the Strong beams are 0.20 m lower than those derived from the Weak beams. In Table 8, for mostly bare ground, the differences between the Strong and Weak beams get 2-times smaller, and are ~0.10 m. Elevations sampled by the Strong beams are still lower than those derived from the Weak beam data. These mean differences in terrain elevations between Strong and Weak channels seems to indicated that the Weak spots may be biased when sampling the ground elevations, slightly overestimating the terrain heights. Table 8. Statistics for elevation differences SRTM90 -h_te_interp (in m) for h_te_uncertainty less or equal to 7.5 m and less or equal 5% tree cover for all laser beams, Strong (1, 3 and 5) and Weak (2, 4, and 6), and for Day (night_flag = 0) and Night (night_flag = 1).

SUMMARY AND CONCLUSIONS
We show that ATLAS/ICESat-2's ATL08 ground elevation products are of sufficient quality to produce Global Ground Control data to evaluate topographic datasets quality. Editing based on total uncertainty which considers POD and PPD geolocation accuracy, when combined with relatively bare ground cover , further discriminate data suitable for Ground Control. Apparent surface reflectance and signal-to-noise ratio do not appear to have a strong influence on the precision of the ATL08 ground product. To a certain extent, the accuracy of the ATL08 ground elevation improves when the number of ground photons increases. We also see a slight bias of ~0.10 m between the Weak and Strong beams, although there are 2-times more observations for the Strong beams compared to the Weak beams. Further investigation of this bias is required.
Comparisons with elevation differences between ICESat-derived Ground Control against SRTM DEMs in Australia for bare ground elevations are very similar, showing similar biases estimates for the DEMs, and equivalent relationships with land cover and relief (Carabajal et al, 2011). When vegetation is present, the SRTM DEM is sampling heights within the canopy, and shown by the positive differences with the ATL08 ground elevations in the edited data. The abundance of ICESat-2 data compared to its predecessor provides vast topographic information to contribute to the development and assessment of topographic assets.