An Alternative Approach for Registration of High-Resolution Satellite Optical Imagery and ICESat Laser Altimetry Data

Satellite optical images and altimetry data are two major data sources used in Antarctic research. The integration use of these two datasets is expected to provide more accurate and higher quality products, during which data registration is the first issue that needs to be solved. This paper presents an alternative approach for the registration of high-resolution satellite optical images and ICESat (Ice, Cloud, and land Elevation Satellite) laser altimetry data. Due to the sparse distribution characteristic of the ICESat laser point data, it is difficult and even impossible to find same-type conjugate features between ICESat data and satellite optical images. The method is implemented in a direct way to correct the point-to-line inconsistency in image space through 2D transformation between the projected terrain feature points and the corresponding 2D image lines, which is simpler than discrepancy correction in object space that requires stereo images for 3D model construction, and easier than the indirect way of image orientation correction via photogrammetric bundle adjustment. The correction parameters are further incorporated into imaging model through RPCs (Rational Polynomial Coefficients) generation/regeneration for the convenience of photogrammetric applications. The experimental results by using the ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) images and ZY-3 (Ziyuan-3 satellite) images for registration with ICESat data showed that sub-pixel level registration accuracies were achieved after registration, which have validated the feasibility and effectiveness of the presented approach.


Introduction
Remote sensing observation of the Antarctic ice sheets, compared to traditional in situ measurements, is an important and relatively efficient method to acquire the spatial information of the ice-sheet topography as well as its change, and provides essential information for the study of global climate change and sea level rise [1][2][3][4]. Satellite optical images and altimetry observations are the two main data sources used in Antarctic remote sensing. A number of radar altimetry satellites, such as Seasat (Sea Satellite), Geosat (Geodetic Satellite), ERS (European Remote Sensing Satellite), and Envisat (Environmental Satellite), had found applications in Antarctic research since its emergence. However, the radar altimetry data has accuracy challenge over rough topographic areas and small mountain glaciers or ice caps [5]. Relatively, laser altimetry can provide higher accuracy observation, e.g., the norminal accuracy of the ICESat (Ice, Cloud, and land Elevation Satellite) altimetry data is as high as 0.15 m [4]. However, the sparse distribution of the ICESat points leaves quite a large Following the detailed methodology of the approach, comprehensive experiments using the ASTER images and ZY-3 images for registration with ICESat altimetry data were conducted and the registration results are demonstrated and discussed, leading to the conclusion in the final section. Figure 1 illustrates the framework of the methodology for the registration of satellite optical images and ICESat altimetry data. The method consists of three main steps: (1) extraction of the matching primitives including 3D feature points automatically extracted from ICESat data and corresponding 2D feature lines interactively extracted from optical images; (2) projection of the 3D feature points to image space and transformation between the projected points and the corresponding feature lines based on point-to-line distance minimization constraint; and (3) regeneration of the satellite imaging model by calculating the bias-compensated RPCs (Rational Polynomial Coefficients) incorporating the transformation parameters determined in step 2.

Point-to-Line Registration of Satellite Optical Images and ICESat Altimetry Data
Sensors 2016, 16,2008 3 of 17 direct discrepancy correction in image space through 2D transformation between the projected terrain feature points and the corresponding 2D image feature lines. The method can work for mono images and requires no conjugate points, which facilitates the registration of sparse points and the optical images. Following the detailed methodology of the approach, comprehensive experiments using the ASTER images and ZY-3 images for registration with ICESat altimetry data were conducted and the registration results are demonstrated and discussed, leading to the conclusion in the final section. Figure 1 illustrates the framework of the methodology for the registration of satellite optical images and ICESat altimetry data. The method consists of three main steps: (1) extraction of the matching primitives including 3D feature points automatically extracted from ICESat data and corresponding 2D feature lines interactively extracted from optical images; (2) projection of the 3D feature points to image space and transformation between the projected points and the corresponding feature lines based on point-to-line distance minimization constraint; and (3) regeneration of the satellite imaging model by calculating the bias-compensated RPCs (Rational Polynomial Coefficients) incorporating the transformation parameters determined in step 2.

Extraction of Matching Primitives
Because the ICESat laser points are sparse, it is difficult and even impossible to identify the individual conjugate points from both the ICESat data and the optical images. However, both datasets contain terrain information, thus it is possible to obtain terrain features from both of them [43]. Terrain feature points such as abrupt slope changes, can be extracted from the elevation profiles along ICESat tracks. These feature points often correspond to ridge crests or foot boundary points, as illustrated in Figure 2a. The terrain feature points are calculated through curve fitting using the adjacent several laser points, thus they are more robust and accurate than direct laser observation points [43]. This is implemented automatically by setting a threshold value of slope change, and all the terrain feature points with a slope change value larger than the threshold will be automatically extracted. Then these terrain feature points are projected into the image space providing initial

Extraction of Matching Primitives
Because the ICESat laser points are sparse, it is difficult and even impossible to identify the individual conjugate points from both the ICESat data and the optical images. However, both datasets contain terrain information, thus it is possible to obtain terrain features from both of them [43]. Terrain feature points such as abrupt slope changes, can be extracted from the elevation profiles along ICESat tracks. These feature points often correspond to ridge crests or foot boundary points, as illustrated in Figure 2a. The terrain feature points are calculated through curve fitting using the adjacent several laser points, thus they are more robust and accurate than direct laser observation points [43]. This is implemented automatically by setting a threshold value of slope change, and all the terrain feature points with a slope change value larger than the threshold will be automatically extracted. Then these terrain feature points are projected into the image space providing initial locations for guidance of corresponding image feature lines extraction. The corresponding image feature lines are extracted interactively with the aid of edge detection algorithms such as the Canny algorithm [44]. Due to the difference of radiation intensity caused by terrain fluctuation and surface material variation, some feature lines may be not distinct and of high uncertainty. These feature lines will be abandoned and only those distinct and of high certainty will be kept, as illustrated in Figure 2b. These feature lines correspond to ridge crest lines or foot boundary lines. The length of the extracted feature lines should be long enough that after transformation, the transformed feature points will coincide with the corresponding feature lines within the expected error limit.
It should be noted that the ground features, especially those ice-rock boundary lines, may have moved due to the glacier change between the ICESAT measurements and the image acquisitions. In order to minimize the registration error introduced by possible glacier change between the acquisition times of different datasets, we chose the ICESat observations which were obtained close to the time or in the same season of the satellite image acquisition for the extraction of conjugate terrain features. In addition, a good distribution of the terrain features in the full image extent contributes to a robust and higher registration accuracy. locations for guidance of corresponding image feature lines extraction. The corresponding image feature lines are extracted interactively with the aid of edge detection algorithms such as the Canny algorithm [44]. Due to the difference of radiation intensity caused by terrain fluctuation and surface material variation, some feature lines may be not distinct and of high uncertainty. These feature lines will be abandoned and only those distinct and of high certainty will be kept, as illustrated in Figure 2b. These feature lines correspond to ridge crest lines or foot boundary lines. The length of the extracted feature lines should be long enough that after transformation, the transformed feature points will coincide with the corresponding feature lines within the expected error limit. It should be noted that the ground features, especially those ice-rock boundary lines, may have moved due to the glacier change between the ICESAT measurements and the image acquisitions. In order to minimize the registration error introduced by possible glacier change between the acquisition times of different datasets, we chose the ICESat observations which were obtained close to the time or in the same season of the satellite image acquisition for the extraction of conjugate terrain features. In addition, a good distribution of the terrain features in the full image extent contributes to a robust and higher registration accuracy.

Discrepancy Correction and Data Registration
The registration of ICESat laser points and optical images is based on the constraint that the projected terrain feature points extracted from ICESat data should coincide with the corresponding 2D feature lines measured from optical images. The 3D terrain feature points are first projected into image space via the vendor-provided imaging model-e.g., the RPCs provided along with the images. Due to the platform diversity and the errors in satellite orbit and attitude determination as well as inaccuracy of sensor calibration, there exist discrepancies between the projected terrain feature points and the measured feature lines in image space.
Generally, there are two ways to eliminate the geometric discrepancies. One is to correct the resulting discrepancies directly in image space or object space, which can be referred to as the direct method. The other is to correct the imaging model parameters through photogrammetric bundle adjustment to achieve geometric consistency between the two types of datasets, which can be referred to as the indirect method. In the latter method, the parameters and observations which cause the discrepancies are adjusted, including satellite attitude and orbit observations (referred to as exterior orientation) and sensor calibration parameters (referred to as interior orientation). The indirect method requires an initial rigorous sensor model and is of high computational complexity. It is impracticable if only the RPCs are available. The direct method, on the other hand, is simpler for comprehension and easier for computation, and is applicable for both a rigorous sensor model and RPC model. According to the space where the discrepancies considered and corrected, there are two schemes for the direct method. One is to correct the discrepancies in image space and the other is in object space, the latter one was adopted in the research in [43]. The object-space-based scheme requires stereo images to retrieve the 3D information and may fail to work with mono-view images, while it is of no problem for the image-space-based scheme. Discrepancy correction in image space performs much better than that in object space in the case of multisource image data [45].
The discrepancies between the projected terrain feature points and the measured feature lines in image space can be modeled and compensated by polynomial transformation with constraint of minimizing the distance between points and corresponding lines. Let us say that (x i , y i ) are the image coordinates of the ith projected terrain feature point, (x ci , y ci ) are the image coordinates after transformation, and a i x + b i y + c i = 0 is the corresponding feature line measured in image space.
We have where k xj and k yj (j = 0,1,2) are the polynomial transformation parameters. According to the selection of the parameters, we have different transformation models as follows. (a) , translation and scale model with four parameters.
x i y i , affine model with six parameters.
Generally, higher order polynomials are not necessary and the affine model is efficient enough to model and compensate the discrepancies [46,47].
For each pair of matched features, we have distance observation equation: In Equation (2), a i , b i , c i are the known coefficients of the extracted image feature lines. Combine Equations (1) and (2) to replace (x ci , y ci ) with image coordinates (x i , y i ) of projected terrain feature point and the unknown transformation parameters k xj and k yj . The transformation parameters can be solved using the least-square minimization constraint min( i ), and n should be more than the number of model parameters for a steady and accurate solution.
For n pairs of matched features, the distance observation equations can be written in the form of: where X represents the vector of the to-be-solved transformation parameters k's, A is the corresponding coefficient matrix, L is the vector of the remaining invariants, and V represents the vector of residuals.
If we have n points, and k's are the affine model parameters, the constructions of A, X and L are as follow: a n a 2 n + b 2 n a n x n a 2 n + b 2 n a n y n The least-square solution of Equation (3) is: In order to remove the possible outliers for a robust estimation, 3-sigma (3 σ) criterion is applied and those features with residuals more than 3 σ (σ = V T V n − t , t is the number of transformation parameters) will be abandoned [48].
With application of the transformation model, the inconsistency between the two datasets can be minimized within an expected error to achieve data registration.

Imaging Model Regeneration for Optical Images
With the discrepancies being compensated by direct correction in image space, registration of the two datasets is achieved. However, the correction model always has to be provided along with the image, which is not convenient for applications. In order to incorporate the correction parameters, the original imaging model has to be updated.
Compared with the rigorous sensor models, the RFM (Rational Function Model, also called RPC model) is more popular and supported by most satellite image vendors and photogrammetric software. Therefore, bias correction can be incorporated into the photogrammetric model by generation/regeneration of a new set of bias-free RPCs. There a lot of studies have been published on RFM-based geometric performance evaluation and bias correction for various high-resolution satellite images [49][50][51][52][53]. The generation of a new set of bias-free RPCs includes the following steps [54]. (1) Establishment of a three-dimensional (3D) lattice in ground space. The lattice contains m × n × k (m and n represent the number of rows and columns in horizontal plane and k is the number of elevation layers) virtual object points. The dimension of the lattice should cover the range of the 3D terrain surface within the full extent of the image. To achieve a robust and accurate solution, the number of elevation layers should be greater than three, and the total number of object points should be adequate enough (e.g., 100) considering 78 to-be-solved RPC parameters.

Data Used
The  [56,57]. Considering that the vertical and geolocation accuracy of the ICESat data is higher than that of the satellite images, the ICESat data is treated as reference, and the satellite images is going to be co-registered to it. The satellite images overlayed with the ICESat data are shown in Figures 3 and 4, which reveal a quite sparse distribution of the ICESat data, with large blank areas lacking altimetry observation between ICESat tracks.

Matching Primitives Extraction
The matching primitives, including the terrain feature points and the corresponding feature lines, were extracted from the ICESat altimetry data and the optical images respectively following the principle described in Section 2.1. As the ASTER images were acquired in January 2005, and the ZY-3 images were acquired in March 2014, the terrain feature points used as control for registration were mainly extracted from the ICESat laser points which were acquired in the same season close to January for ASTER and March for ZY-3. In the experiments, for the registration of ASTER images and ICESat data, nine pairs of terrain feature points and corresponding feature lines were extracted as control for each image, and another eight pairs were extracted as check for independent accuracy evaluation for each image. For the registration of ZY-3 images and ICESat data, the number of control features and check features are 10 and 8, respectively. The two images in each stereo pair share the same ground points for both the control features and check features. Figures 5 and 6 show the distribution of the control features and the check features for the ASTER images and ZY-3 images, respectively.

Matching Primitives Extraction
The matching primitives, including the terrain feature points and the corresponding feature lines, were extracted from the ICESat altimetry data and the optical images respectively following the principle described in Section 2.1. As the ASTER images were acquired in January 2005, and the ZY-3 images were acquired in March 2014, the terrain feature points used as control for registration were mainly extracted from the ICESat laser points which were acquired in the same season close to January for ASTER and March for ZY-3. In the experiments, for the registration of ASTER images and ICESat data, nine pairs of terrain feature points and corresponding feature lines were extracted as control for each image, and another eight pairs were extracted as check for independent accuracy evaluation for each image. For the registration of ZY-3 images and ICESat data, the number of control features and check features are 10 and 8, respectively. The two images in each stereo pair share the same ground points for both the control features and check features. Figures 5 and 6 show the distribution of the control features and the check features for the ASTER images and ZY-3 images, respectively.

Registration Results
By using the extracted control features, the transformation model parameters can be determined for each image by least squares adjustment as introduced in Section 2.2. Different scenarios of parameterization were tested to compare the effectiveness of the four transformation models indicated in Section 2.2 and the results are listed in Table 1, from which we can see that before registration, the average discrepancies are more than four pixels between the ASTER images and the ICESat data. Among the four transformation models, the affine model produced the highest accuracy to sub pixel level after registration. The ZY-3 images showed relatively smaller inconsistency with the ICESat data at around 1-2 pixels before registration, and it was further improved to sub-pixel level after registration by using the similarity or affine transformation model. Figure 7 illustrates an example showing the geometric relation between an ICESat track and the corresponding ASTER image feature lines before and after registration, which reveals that a higher consistency between the ICESat terrain feature points and the corresponding image feature lines is achieved after registration.  (2)) between feature lines and projected feature points before and after registration using different transformation models (unit: pixel).

Registration Results
By using the extracted control features, the transformation model parameters can be determined for each image by least squares adjustment as introduced in Section 2.2. Different scenarios of parameterization were tested to compare the effectiveness of the four transformation models indicated in Section 2.2 and the results are listed in Table 1, from which we can see that before registration, the average discrepancies are more than four pixels between the ASTER images and the ICESat data. Among the four transformation models, the affine model produced the highest accuracy to sub pixel level after registration. The ZY-3 images showed relatively smaller inconsistency with the ICESat data at around 1-2 pixels before registration, and it was further improved to sub-pixel level after registration by using the similarity or affine transformation model. Figure 7 illustrates an example showing the geometric relation between an ICESat track and the corresponding ASTER image feature lines before and after registration, which reveals that a higher consistency between the ICESat terrain feature points and the corresponding image feature lines is achieved after registration.  (2)) between feature lines and projected feature points before and after registration using different transformation models (unit: pixel).  As discussed in Section 2.3, the transformation model needs to be fused into the imaging model through RPCs generation/regeneration for the convenience of applications. Based on the original imaging models and the derived affine transformation model, RPCs were generated for each optical image. In order to validate the accuracy improvement of the regenerated RPCs in photogrammetric applications, two DEMs (Digital Elevation Models) were generated from the stereo images by aid of an ERDAS (Earth Resource Development Assessment System) software, one using the original imaging model and the other using the regenerated bias-free RPCs. The two DEMs were then compared with the ICESat observations for consistency assessment before and after registration. Figures 8 and 9 illustrate two elevation profiles along the ICESat tracks for the ASTER image and ZY-3 image, respectively, from which we can see that significant offset biases exist between the ICESat elevations and the photogrammetric elevations before registration, and a higher consistency was achieved between the ICESat data and the elevations derived from the optical images using the regenerated bias-free RPCs after registration. The mean and standard deviation of the elevation discrepancies before and after registration are listed in Table 2. From Table 2, we can see that the mean elevation discrepancies were significantly reduced from −61.4 m to 5.9 m for ASTER and from −28.3 m to −1.4 m for ZY-3, both better than the corresponding image resolutions of 15 m and 2.1 m/3.5 m, respectively. As regards the standard deviation, the improvement for ASTER was greater than for ZY-3, revealing that the geometric discrepancy between the ICESat and ASTER data is more complex than that between the ICESat and ZY-3 data before registration and more parameters are needed for the correction of ASTER data, which is in agreement with the results obtained in Table 1. As discussed in Section 2.3, the transformation model needs to be fused into the imaging model through RPCs generation/regeneration for the convenience of applications. Based on the original imaging models and the derived affine transformation model, RPCs were generated for each optical image. In order to validate the accuracy improvement of the regenerated RPCs in photogrammetric applications, two DEMs (Digital Elevation Models) were generated from the stereo images by aid of an ERDAS (Earth Resource Development Assessment System) software, one using the original imaging model and the other using the regenerated bias-free RPCs. The two DEMs were then compared with the ICESat observations for consistency assessment before and after registration. Figures 8 and 9 illustrate two elevation profiles along the ICESat tracks for the ASTER image and ZY-3 image, respectively, from which we can see that significant offset biases exist between the ICESat elevations and the photogrammetric elevations before registration, and a higher consistency was achieved between the ICESat data and the elevations derived from the optical images using the regenerated bias-free RPCs after registration. The mean and standard deviation of the elevation discrepancies before and after registration are listed in Table 2. From Table 2, we can see that the mean elevation discrepancies were significantly reduced from −61.4 m to 5.9 m for ASTER and from −28.3 m to −1.4 m for ZY-3, both better than the corresponding image resolutions of 15 m and 2.1 m/3.5 m, respectively. As regards the standard deviation, the improvement for ASTER was greater than for ZY-3, revealing that the geometric discrepancy between the ICESat and ASTER data is more complex than that between the ICESat and ZY-3 data before registration and more parameters are needed for the correction of ASTER data, which is in agreement with the results obtained in Table 1.

Discussion
Due to the difference in platform and systematic errors inherent in sensor calibration and orientation observations, geometric inconsistency exists between two different types of datasets, even between datasets obtained at different times from the same sensor, for which data registration is an essential procedure to eliminate the geometric inconsistency for integration use. Generally, translation could compensate a majority part of the inconsistency, which is corresponding to the discrepancies induced by constant sensor calibration errors and constant orbit and attitude measurement errors. However, there may also exist additional discrepancies after translation

Discussion
Due to the difference in platform and systematic errors inherent in sensor calibration and orientation observations, geometric inconsistency exists between two different types of datasets, even between datasets obtained at different times from the same sensor, for which data registration is an essential procedure to eliminate the geometric inconsistency for integration use. Generally, translation could compensate a majority part of the inconsistency, which is corresponding to the discrepancies induced by constant sensor calibration errors and constant orbit and attitude measurement errors. However, there may also exist additional discrepancies after translation

Discussion
Due to the difference in platform and systematic errors inherent in sensor calibration and orientation observations, geometric inconsistency exists between two different types of datasets, even between datasets obtained at different times from the same sensor, for which data registration is an essential procedure to eliminate the geometric inconsistency for integration use. Generally, translation could compensate a majority part of the inconsistency, which is corresponding to the discrepancies induced by constant sensor calibration errors and constant orbit and attitude measurement errors. However, there may also exist additional discrepancies after translation compensation, which could be attributed to the drift and/or higher errors with orbit and attitude measurements [45,49]. Therefore, to compensate the additional discrepancies, more parameters are needed for correction, such as translation and scales, similarity transformation, and affine transformation, as tested in the experiments in Section 3.3. Theoretically, the more parameters used in transformation, the smaller residuals it will produce, and the larger number of control features it requires. However, the transformation selected should be as close as possible to the physical model in order to avoid an over-fitting problem. Generally, if only the residuals meet the accuracy expectation, the number of parameters should be as smaller as possible to avoid over fitting. In our experiments, the expected registration accuracy is sub-pixel level. Therefore, for the registration of ASTER images and ICESat data, the affine transformation model is required to achieve the expected sub-pixel accuracy. For the registration of ZY-3 images and ICESat data, however, the similarity transformation model is adequate enough, although the accuracy could be further improved by using the affine model.
With the geometric inconsistency being corrected by the transformation model in image space, the registration of the two datasets is achieved. However, the use of additional transformation model is not convenient in application such as DEM generation. Therefore, the transformation parameters are incorporated into the imaging model through RPCs generation/regeneration. The results of DEM discrepancy analysis along ICESat tracks conducted in Section 3.3 validated the effectiveness of the updated imaging models.
A similar terrain-feature-based method was presented to register the historical aerial images to ICESat data in [43]. It was based on 3D similarity transformation with point-to-line distance minimization constraint in object space to transform the local photogrammetric system to ICESat reference system, requiring that the 3D model of terrain features be first generated stereoscopically from the aerial images. In the experiment, after a rigid body transformation adjustment with a total of 15 characteristic lines and points, the geometric inconsistency between the ICESat data and the 3D model generated from the aerial images were decreased from the original more than 100 m down to several meters. The method in [43] shared similar concept with the method proposed in [19] for the registration of images and lidar data based on linear features. They both start by manipulating the photogrammetric images to produce a 3D model including a set of linear features, and followed by establishment of the transformation between the two coordinate systems based on the control features. The advantage is that they can be applied for the co-registration of multiple three-dimensional datasets regardless of their origin. However, the disadvantage also lies in the requirement of three-dimension of both datasets being considered for co-registration-i.e., they require stereo images for 3D model construction and could not work for mono images. Moreover, the extraction of corresponding linear features from stereo images for 3D construction may encounter difficulty due to occlusions since the linear features are along object space discontinuities [19]. The alternative simplified point-to-line feature matching method presented in this paper is based on direct discrepancy correction in image space for the registration of satellite optical images and ICESat data, which is simpler and easier for implementation, and it can work for individual images, no matter stereo or mono. Meanwhile, discrepancy correction in image space also allows for investigation of the satellite stability-e.g., the affine correction model indicates the potential existence of gyro drift error and radial ephemeris error [46,49]. Because the methods discussed above were also not fully automatic, the implementation of these methods is therefore not conducted for comparison in the experiments since the major differences lie in conception and its limitation rather than the accuracies achievable.

Conclusions
Due to the sparse distribution of ICESat data, it is difficult and even impossible to find individual conjugate points between ICESat data and optical images, an alternative simplified registration method based on point-to-line discrepancy correction in image space is presented for the registration of the two types of datasets. Compared with the indirect way of inconsistency elimination by orientation parameter correction through photogrammetric bundle adjustment, the direct discrepancy correction method is simpler for comprehension and easier for computation. There are also two schemes for direct discrepancy correction according to the space in which it is implemented, among which discrepancy correction in image space is simpler and more generic than that in object space and can work for individual images, whether they are stereo or mono. Moreover, it is suggested that the correction result be incorporated into imaging model via RPCs generation/regeneration for the convenience of applications.
The experimental results of using the ASTER images and ZY-3 images for registration with ICESat data show that sub-pixel level accuracies were achieved after registration. By comparing the DEMs derived from the stereo optical images and the ICESat data through elevation profile analysis along ICESat tracks, the results show a significantly higher consistency by using the regenerated RPCs with the correction model being incorporated after registration than using the original imaging models.
At the current stage, the matching features were extracted interactively. In further research, more effort could be paid to the automation of feature extraction and correspondences matching to improve the efficiency for practical applications.