POINT-TO-SURFACE MATCHING FOR DEM CORRECTION USING ICESAT DATA

Vendor-provided rational polynomial coefficients (RPCs) are commonly used to generate digital elevation models (DEMs) from highresolution satellite images. This results in a level of accuracy that can be improved using ground control points (GCPs). It is well known that due to the inherent bias of sensor orientation the generated DEM is distorted. In the traditional way of working, the bias is corrected by integrating GCPs into the standard processing chain. This involves additional effort, since the provision of GCPs and the measurement of their image coordinates are required. In this paper, we examine whether and how the data recorded by NASA's ICESat (Ice, Cloud, and Land Elevation Satellite) mission can be used as GCPs without measuring image coordinates. The starting point are DEMs that are generated by image matching from KOMPSAT-3 satellite images with given RPCs. We developed a point-to-surface matching method that matches the ICESat points to the DEM in order to correct the DEM and improve its precision. For the experimental investigations a KOMPSAT 3 image pair is used that covers an area of 20 by 16 km in the Yangsan city regions. The generated DEM has a height accuracy of about 9 m. The pointto-surface algorithm with 505 ICESat points led to an improvement of the DEM height accuracy to about 2 m.


INTORDUCTION
Photogrammetric processing of high-resolution satellite images for products like digital elevation models (DEM) requires rigorous sensor modelling or highly accurate approximations to it in the form of rational polynomial coefficients (RPCs) which are widely free from sensor orientation biases (Dial, Grodecki, 2005). The RPC is used widely for convenient geo-positioning. As demonstrated for many high-resolution imaging satellites using RPC geo-positioning, there is ample experimental evidence that the relative accuracy to meter level is attainable. These methods require one or more ground control points (GCPs). The provision of GCPs can be time consuming, particularly if fieldwork, such as GPS surveying, is required. In traditional approaches, the GCPs need to be measured in the images to make them usable for these methods. Ground point determination to 0.3 pixels is possible under ideal conditions of well defined and recognizable ground points, precise image measurements, and the provision of sensor calibration data (Ebner et al., 1996). Accuracies between 0.5 and 2 pixels are often found in practical work. Based on rigorous models, concepts are being pursued without relying on conventional GCPs. Bouillon et al. (2006) performed bundle block adjustments using tie points without GCPs to improve the quality of DEM generation from SPOT-5 HRS stereo images. The horizontal and vertical DEM accuracy for 90% of the points were 8.4 m and 4.5 m, respectively. The model reported by Toutin et al. (2012) with virtual control points led to elevation errors with 2.6 m and 2.1 m (68% confidence level) for the World-View-1 and 2 stereo images. Other approaches to improve the quality of generated DEMs have their origin in point determination using DEMs as the control information (Rosenholm, Torlegård, 1988;Ebner, Strunz, 2017). Given that the orientation difference between the points and 1 Corresponding author surface is small, the temporary pairing of the points with points on the surface is based on the same horizontal position. By undertaking pairing and minimization of the height differences iteratively, every iteration brings the points closer to the surface. This method is also called the least-height difference (LHD) algorithm because the sum of the squared height differences is minimized to estimate the transformation parameters. Kim and Jeong (2011) examined the suitability of the LHD algorithm for the precise mapping of push broom images. They reported that the 3D similarity transformation is evident when errors occur only in the form of the time-invariant position and attitude biases of the coordinate frame of the push broom images. Chen et al. (2017) employed a shuttle radar topography mission (SRTM) DEM for an RPC correction of images from TianHui-1 satellite based on the LHD algorithm. Cao et al. (2019) evaluated the reliability of three types of existing global DEM to improve the accuracy of the sensor model of numerous images of the ZY-3 satellite based on the bundleadjustment with the LHD algorithm. In this paper we investigate how the data recorded by the NASA satellite mission ICESat (Ice, Cloud, and Land Elevation Satellite) can be used as GCPs. The ICESat data are freely accessible and can avoid collecting GCPs. The ICESat was launched in January 2002 and collected laser altimeter data until 2009 (NASA, 2020). The sole instrument on ICESat was the Geoscience Laser Altimeter System (GLAS). GLAS produces a series of approximately 70 m diameter laser spots that are separated by approximately 170 m along the spacecraft's ground track. The primary objective of ICESat is a determination of the interannual and long-term changes in polar ice-sheet volume to sufficient accuracy to assess their impact on global sea level (Zwally et al., 2002). The accuracy of the elevation measurement over the polar sea-ice was less than ±14 cm (Shuman et al., 2006;Kurtz et al., 2008). The secondary objectives include a The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B4-2020, 2020 XXIV ISPRS Congress (2020 edition) measurement of cloud and aerosol height profiles, and important for our investigation, the measurements of land elevation and vegetation cover. The location data derived from the laser footprint of ICESat can separate canopies and ground reflections. Hence, the data are also used to identify the vegetation height of a forest (Lefsky et al., 2005). ICESat was decommissioned in 2010, but a follow-on mission, ICESat-2, was developed with the launch of the satellite in 2018.
In this research, we use a point-to-surface matching algorithm that we developed for the correction of DEMs generated from KOMPSAT satellite images with given RPCs (Lee and Hahn, 2019). With this algorithm we match the ICESat points to the DEM with the idea to correct the DEM and improve its precision. The process is fully automated -no interaction with an operator is necessary. Because the number of ICESat points is comparatively small, matching is not very time-consuming even with many iterations. The matching result between the ICESat points and the DEM is represented by the estimated parameters of the geometric transformation. The DEM is then corrected using the transformation parameters.
The use of ICESat points should allow an improvement in DEM accuracy. However, it must be ensured that only those GCPs are included in the point-to-surface adaptation where vegetation or other influences do not lead to deviations of the DEM from the surface of the terrain. For this reason, we also analyze the use of a land cover map to investigate the influence of ICESat data on the point-to-surface matching.

POINT-TO-SURFACE MATCHING MODEL
The objective of the matching approach is to find the transformation between the points and the surface, which consists of a 3D rotation and a 3D translation, as well as a scale parameter, such that the sum of the squared distances between the points and surface is minimized. For the LHD algorithm, the height difference between point and point * on the DEM is used as the measure of distance. Instead of the height difference, we use the shortest distance, , from point to the surface at point , as shown in Figure 1. The determination of the shortest distance, , from point to the surface can be formulated as an optimization problem of finding the corresponding point . If the surface is described by a plane (as sketched in Figure 1), the process is straightforward because any point on the plane (e.g. * ) can be used to calculate the distance, , of point from the plane according to where indicates the i th point; is a normal vector of length, 1; is the location of the GCP, and is the location of the corresponding point on the surface. For non-planar surfaces, patches, like bilinear surfaces, the optimization can be solved iteratively by approximating the tangential planes, for which the surface point that best matches can be found iteratively.
The 3D similarity transformation model is as follows: * = + (2) * = + ∆ is the location of the transformed GCP; = ( , , ) is the 3D rotation matrix; is a 3D translation vector, and s is a scale parameter. The least-squares approach for point-to-surface matching follows directly by taking points outside the surface (in the present case, the GCPs) into account. The GCPs substitutes the surface points, and replaces 0 . If the surface normal is normalized to length 1 ( /‖ ‖), Equation (1) calculates the signed Euclidean distance to the tangential plane. For the least-squares approach, suitable approximate values are needed for the parameters of the geometric transformation. In matching the GCPs to a DEM, a rotation matrix = and a translation vector = are suitable starting values. Using the notation * for iteratively transformed GCPs.
The linearized model of Equation (3)  , results in the standard least-squares estimates of the parameters, ̂= ( ) − . The weights shall be determined appropriately. If the horizontal displacement between the DEM and the GCPs is greater than the DEM grid width, the iterative estimation of the parameters causes the transformed GCPs to move across the adjacent grid cells. The normal vector may differ significantly between the adjacent DEM grid cells. An experimentally observed consequence was that the estimates were obtained that transform the GCPs into adjacent grid cells during a single iteration and transforms them back in the following iteration. To mitigate this influence, only one-third of the iteratively estimated values (̂/3 ) was taken into consideration when updating the parameters. This smooths the iteration behaviour within the least-squares estimation process but causes an increase in the number of iterations. After a successful match, elevation residuals remain between ICESat and DEM. In high-relief areas with a forest-canopy, the elevation difference increases between ICESat and SRTM DEM, whereas the elevation difference decreases in areas of low-relief and sparse tree-cover because of the SRTM data produced by radar scattering over the canopies Bhang et al., 2007). A land cover classified image was used to avoid the effects of the tree and slope in a mountain area while matching. Figure 2 gives an overview of the concept and processing steps underlying the experimental investigations with the KOMPSAT-3 imagery.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B4-2020, 2020 XXIV ISPRS Congress (2020 edition) Figure 2. Concept and workflow of using ICESat data for correcting KOMPSAT-3 DEM

EXPERIMENT AND RESULTS
The test site used is Yangsan City, which is located in the east region of the Korean Peninsula ( Figure 3). The KOMPSAT-3 images cover an area of approximately 20 km by 16 km in the Yangsan City region. The test site offers diverse possibilities to determine the suitability of the ICESat data. The ICESat global elevation product type GLAH06 was used in the experiment without additional waveform analysis for vegetation or ground reflection. The elevation data obtained in February, March, June, September, October, and November from 2003 to 2008 were used for the experiment. Therefore, it can be estimated that the data collected during February and March were affected slightly by the canopy and vegetation compared to the other seasons. Here, the data in the winter season (February and March) is half a percent. A total of 505 ICESat points were available for the point-to-surface matching investigations ( Figure  4). In this region, 13 national control points (NCPs) were available for use in the experiment to obtain an independent estimate of the points to the DEM transformation parameters. For the accuracy assessment of the corrected DEM, they further serve as ground truth. The NCPs were established throughout South Korea by the Korean National Geographic Information Institute (NGII, 2020). They result from geodetic measurements with GNSS and geodetic network processing. Moreover, 15 checkpoints were recorded by a GPS survey within the region of the DEM for ground truth investigations. The NCPs and GPS points are used as independent checkpoints for accuracy assessment of the matching. These points are well recognizable in the KOMPSAT-3 image pair. Figure 4 presents the ICESat points, NCPs, and GPS points located in the project area.  A DEM with a 5 m grid spacing is generated from the KOMPSAT-3 image pair and the vendor-provided RPCs using the ERDAS IMAGINE Photogrammetry Suite in the test site ( Figure 7). The ERDAS image matching algorithm produces a DEM that does not represent the ground in some regions but rather approximates a surface above the ground. Since the laser measurements have the potential to detect the ground below the tree cover in winter, the height of the ICESat can differ from the ERDAS DEM at the same horizontal location, which is visible in mountain areas via the height of the canopy. Even with data outside the winter season, it is difficult to analyze whether the height relates to the canopy or the terrain surface. To avoid the effects of tree cover and topographic relief in a mountain area, a land-cover classified image was produced to select only the points on the ground. Figure 5 shows a classification image using ERDAS S/W. The classification was performed using the supervised maximumlikelihood method with a sampling forest, water, urban, and ground area. ICESat points related to the ground were selected by the land cover classified image. Figure 6 shows that the selected 76 ICESat points are located in low-elevation areas of a surrounding ground area. In contrast, most of the points that are in exposed locations in the mountains have been removed. To apply point-to-surface matching, an experiment was performed with the following two cases of transformation parameters: Case 1: Parameters using all ICESat points Case 2: Parameters using 76 ICESat points by applying the classified image The iteration performance of the matching procedure in both cases was investigated (Figure 8). Figure 8 shows the shifts and rotation angles resulting from the estimation process that are updated at each iteration. In Case 1, all parameters converged at the 15 th iteration, while the rotation angles of Case 2 required a few more iterations until convergence of the algorithm was reached. In particular, the movement of kappa ( ) changed greatly during the iteration process compared to other parameters in two cases. Figure 8 (3 rd column) shows the convergence of the process, in which the mean distance of the transformed points to the surface decreases from 7.8 m and 7.0 m at the beginning to less than 5.5 m and 2.1 m at the end of the iterations, respectively, for Cases 1 and 2. The error in Case 1 was obviously larger than that of Case 2.
(a) Case 1 (b) Case 2 Figure 8. Convergence behaviour of the matching algorithm Table 1 lists the estimated parameters of point-to-surface matching for both cases. For the estimated shifts and rotation angles there are much smaller standard deviations in Case 1 than in Case 2, which is due to the different number of points.

Parameters
Case 1 Table 1. Transformation parameters and the standard deviations of point-to-surface matching Figure 9 shows the 3D distances between 505 ICESat points and the surface points before and after matching for Cases 1 and 2. Overall, the distance error is reduced significantly after matching both cases compared to that before matching. The distance errors between 505 ICESat points and the surface points were calculated using the estimated transformation parameters in Table 2. The RMSE values in both cases after matching were approximately 3 m smaller than those before matching in the distance. In both cases, the RMSE and Max values showed similar results. All points for the accuracy comparison were mostly mountain regions, which involves components of forest-canopy heights and high-relief.
(a) Before matching (b) Case 1 (c) Case 2 Figure 9. Visualization of the distances between 505 ICESat points (red circles) and the surface points (blue triangles) Distance error (m) The points in the KOMPSAT-3 image pair can therefore be easily inspected. If matching performed well, there should be almost no deviation in the distance between the 28 points and the surface points of the DEM. Figure 10 plots the distance errors for the 28 points. Table 3 lists the RMSE and the maximum errors of the distance between the NCPs-GPS points and DEM.
The distance values of the NCPs-GPS points to the surface at the beginning of the point-to-surface matching were 9.1 m with a maximum distance of 12.9 m. At the end of the iterations, the RMSE was reduced to 4.2 m and 1.6 m, and the maximum distance was 7.1 m and 3.0 m, respectively, for Cases 1 and 2. This proves that the parameters using the selected ICESat points on the ground from the land cover image were more reliable than those not using the land cover image. Because half of the dataset was collected in winter, there is the possibility that the terrain surface height under the canopy can be obtained in mountainous regions at this time. In contrast, the DEM of ERDAS is produced in the height of trees in the mountainous regions. Therefore, it seems that a height difference occurred, which affects the accuracy of the transformation parameters obtained from Case 1.

CONCLUSIONS
This paper proposes a point-to-surface matching procedure to correct a DEM using ICESat points, which is generated from KOMPSAT-3 satellite images with vendor-provided RPCs. The point-to-surface matching algorithm minimizes the sum of the squared distances between the points and the tangential planes at the correspondence points. The points used as GCPs are the ICESat location dataset provided by NASA.
In the experimental investigations, the land classified image was used in the matching algorithm to correct the DEM. The accuracy of matching DEM and ICESat points led to a more accurate result by exploiting the classification image.
The DEM before the correction deviated from the ground truth points by about 9 m in height. Using the proposed method, an improvement could be achieved by reducing the displacements to an height accuracy level of 2 m. This accuracy is achieved by including only those GCPs on the ground, in which vegetation or other influences do not lead to deviations of the DEM from the terrain surface.