ICESAT-2 LASER ALTIMETRY DATA-GUIDED HIGH-ACCURACY POSITIONING OF SATELLITE STEREO IMAGES

: Positioning accuracy is a key indicator that restricts the accuracy of satellite image mapping. Laser altimetry acquires surface height in decimeters or higher accuracy, providing high-accuracy topographic information for satellite photogrammetry. This study proposes a four-step high-accuracy positioning method called laser altimetry data-guided stereo images positioning (LGSP): photon outlier removal, digital surface model (DSM) generation by a variant of the tube-shape semi-global matching (tSGM), fast point cloud alignment, and DSM and rational function model (RFM) refining. First, outliers of the laser altimetry data are filtered out via confidence selection and open access DSM constraints. Second, the DSM around high-accuracy photons is generated via the variant of tSGM. Third, fast point cloud alignment is performed to align the DSM with laser altimetry data. Finally, we propose a high-resolution stereo model refinement method based on alignment parameters. Experiments on two satellite stereo datasets with different loads, resolutions, and topographies show that LGSP improved the planimetry accuracy of satellite stereo images after adding laser altimetry data. What’s more, LGSP significantly improved the height accuracy of satellite images. Specifically, LGSP reduced the height errors from the original 5.99 m to 1.00 m for the Ziyuan-3 stereo dataset and from the original 20.39 m to 0.83 m for the Gaofen-7 stereo dataset, and slightly outperformed bundle adjustment in which four ground control points were deployed around the four corners of the stereo images. RFM refining based on alignment parameters has been made publicly available at https://mega.nz/file/EFQASJ6A#BKjoupHVshylQR4wCxGd9BDtg08bhy-_aY-c_rHnl-0.


INTRODUCTION
Satellite photogrammetry plays an influential role in the fields of the social economy, ecological protection, and national defense (Tu et al., 2021). The positioning accuracy of satellite images, which has a significant impact on the subsequent processing and application, is a key indicator of satellite photogrammetry (Zhang et al., 2021). Because of the positioning error on the satellite, it is often necessary to use external controls to improve the positioning accuracy of satellite images (Pan et al., 2021). Ground control points (GCPs), the high-resolution digital orthophoto map, and digital elevation model (DEM) are commonly used for external control, but collecting these data is time-consuming and has high acquisition costs (Toutin, 2004). Therefore, open access data such as OpenStreetMap, google earth, and laser altimetry data are becoming a new trend as external controls.
Laser altimetry acquires surface height in decimeters or higher accuracy, providing highly precise topographic information for satellite photogrammetry (Di et al., 2012). Depending on the application's purpose, satellite image positioning based on laser altimetry data can be divided into two categories: 1) laser altimetry data are used as vertical control to eliminate the systematic errors of digital surface model (DSM) or DEM generated by satellite images, and 2) combined adjustment of laser altimetry data and satellite images.
Co-registration of two or more DSMs or DEMs is an effective method for eliminating the systematic errors of DSM or DEM (Gruen and Akca, 2005). Specifically, laser altimetry data were commonly used for the reference data to align the DSM generated by the stereo images. Lin et al. developed a coregistration process to co-register the three digital terrain models from multiple sources (Lin et al., 2010). Di et al. employed surface matching for DEMs generated from stereo images and laser altimetry data separately to achieve coregistration of Chang'E-1 stereo images and laser altimeter data. Then, the image sensor model was refined by adding attitude angle bias correction and evaluated by matching points instead of ground validation (Di et al., 2012). While co-registration is extremely valuable for generating refined DSMs, improving the positioning accuracy of satellite images has a wider range of applications.
Unlike the refined DSM, the combined adjustment of laser altimetry data and satellite images focuses on improving the positioning accuracy of satellite images (Yoon and Shan, 2005). Wu et al. developed a combined block adjustment for Chang'E-1 imagery and laser altimetry data based on local surface constraints (Wu et al., 2011). Liu et al. directly aligned the topographic feature points on laser altimetry data with the image and then enhanced the positioning model of the satellite image (Liu et al., 2016). Wang et al. first filtered the laser altimetry data, then replaced the height of the 3D coordinates obtained from the intersection with the height of the laser points and refined the orientation parameters based on free net adjustment . The aforementioned strategy can improve the height accuracy of satellite images by using the height of laser altimetry data on flat terrain. Li et al. transformed the registration between the laser altimetry data and stereo images to image-to-image matching within the circle of the laser footprint point (Li et al., 2018). The above methods significantly improve the height accuracy of satellite images only aided by the high-accuracy photons. The Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2), launched on 2018, produces a high density of photons in the along-track direction, but these photons have a large number of outliers, which present a new challenge for the combined adjustment of laser altimetry data and satellite images.
In this paper, we propose a four-step high-accuracy positioning method called laser altimetry data-guided stereo images positioning (LGSP), which is robust to outliers in laser altimetry data. First, we perform photon outlier removal on laser altimetry data and then propose a variant of tube-shape semiglobal matching (tSGM) for regions of interest to obtain highaccuracy DSM. After aligning the DSM and laser altimetry data, the DSM generated by dense matching and the high-resolution satellite image model are refined based on the alignment parameters. On the one hand, LGSP is insensitive to the outliers in laser altimetry data, due to using photon outlier removal and farthest points removal in iterative closest point (ICP). On the other hand, LGSP refines the rational function models (RFM) of the satellite stereo images using alignment parameters between laser altimetry data and DSM generated by satellite stereo images.

METHODOLOGY
The flowchart of LGSP is shown in Figure 1: (1) the photon outlier removal method removes inaccurate laser altimetry data to obtain high-accuracy photons.

Photon Outlier Removal
The mission of the ICESat-2, launched on September 15, 2018, is to measure and monitor the impacts of the changing environment (Neumann et al., 2019). The ICESat-2 observatory contains an Advanced Topographic Laser Altimeter System (ATLAS), which sends 10,000 laser light pulses per second to measure the topography, global vegetation, sea ice freeboard, and cloud and atmospheric properties. The spacing of the photons received by ICESat-2 is approximately 70 cm along the satellite's orbit. Compared with the full-waveform measurement ICESat, ICESat-2 measures photons with a higher accuracy of better than 0.1 meters.
ATL03 products distributed by ICESat-2 mainly record the acquisition time, geodetic coordinates (latitude, longitude, and height above the WGS84 ellipsoid), and measurement confidence of each photon (Neumann et al., 2020). The measurement confidence levels of each photon can be classified as 0 (noise), 1 (background), 2 (low confidence signal), 3 (medium confidence signal), or 4 (high confidence signal). The classification method is performed by first calculating the height histogram of photons in a certain range and calculating their signal-to-noise ratio. Photons with a signal-to-noise ratio greater than a threshold are classified as signals, and other photons are classified as background noise. Therefore, most of the outliers in the photons are filtered out during the preprocessing process ( Figure 2). However, as shown in Figure 2 (b), inaccurate photons remain in the red ellipse.
In addition to the aforementioned confidence selection method, TanDEM-X DSM has been used to remove the outliers of photons. TanDEM-X, a global DSM produced by the German Aerospace Center, was produced in September 2016 (Grohmann, 2018;Krieger et al., 2007). Considering that the acquisition time of TanDEM-X and laser altimetry data are close, the height values of photons and TanDEM-X at the same location should also be close. Thus, based on height consistency constraints , photons with large height offsets can be removed.

Variant of tSGM for Regions of Interest
The tSGM algorithm performs semi-global matching (Hirschmuller, 2007) layer by layer on the pyramid image. The disparity search range of each pixel in the next layer is limited by the minimum and maximum values of the valid disparities in the rectangular window around the same position of the current layer (Rothermel et al., 2012). Therefore, the outliers of tSGM's results are lower than that of the classic SGM's results.
Considering that the photons are not uniformly distributed on the image but in strips, only the DSM around high-accuracy photons needs to be generated. When high-accuracy photons are projected onto the epipolar images, the region of interest is generated by the projected points and surrounding pixels. The window size of the surrounding pixels is determined by the direct positioning accuracy of the satellite image. Except for the original stereo images, the pyramid images are matched via tSGM to obtain an accurate disparity map. We propose a variant of tSGM for regions of interest, which comprises two cases. If the pixel of the original image is successfully matched in the upper pyramid and the pixel is outside the region of interest, the disparity ranges are determined based on the disparity pyr d of the pixel in the upper pyramid disparity map and extend to a small range γ .
In other cases, tSGM is performed to ensure matching accuracy. Finally, the matched points for regions of interest are output as the DSM point clouds as shown in Figure 3 (b). Unlike DSM for the Whole Region of traditional tSGM as shown in Figure 3 (a), the variant of tSGM for regions of interest only produces point clouds around laser altimetry data. Thus, the variant of tSGM for regions of interest reduces not only the running time of dense matching but also memory.

Fast Point Cloud Alignment and DSM Refining
When both laser altimetry data and initial DSM have high initial positioning accuracy, only a fine alignment is required. Otherwise, the coarse alignment will be performed by interactively selecting three common points on both surfaces before fine alignment. ICP solves rigid transformation parameters iteratively by least squares estimation, which is simple to implement and has high alignment accuracy (Besl and McKay, 1992). At each iteration, the points of the DSM point clouds that are too far from the laser altimetry data are removed. This farthest points removal alleviates local inconsistencies between laser altimetry data and DSM. Finally, based on the transformation matrix containing three translation parameters and three rotation parameters, all points of the DSM point clouds are transformed to generate the refining DSM.

Refinement of High-resolution Stereo Models Based on Alignment Parameters
The RFM has become the mainstream positioning model for high-resolution satellite images because of its parameter confidentiality, generality, and high fitting accuracy (Tao and Hu, 2001). The RFM establishes the relationship between the image point ( c , r ) and the ground point ( P , L , H ) by using ratios of polynomials: where 0 c , 0 r , 0 P , 0 L , and 0 H are the offset parameters; and s c , s r , s P , s L , and s H are the scale parameters.
The solution of rational polynomial coefficients (RPCs) can be divided into terrain-dependent and terrain-independent strategies. The terrain-dependent strategy uses many control points to solve the RPCs, and the terrain-independent strategy uses a physical sensor model to solve the RPCs. With the physical sensor model available, the RFM can be solved using the object grids with their grid-point coordinates determined using a physical sensor model (Tao and Hu, 2001).
Inspired by the terrain-independent strategy, we propose a method to optimize the high-resolution stereo model based on the alignment parameters. First, many virtual image grids ( c , r ) with uniform distribution in image space are established. Second, several object coordinates ( P , L ) of a virtual grid are calculated by the RFM of the left and right images onto several elevation layers. Third, Applying the alignment parameters transforms the object coordinates ( P , L , H ) of all virtual grids into the aligned object coordinates ( A P , A L , A H ): where R is a rotation matrix, and T is a translation vector. Finally, the RPCs of the left and right images are refined with these aligned object point coordinates ( A P , A L , A H ) and its image point coordinates ( c , r ) by using an iterative least squares solution. As for the iterative least-squares solution, the initial value of the coefficients can be first set using the initial RPCs to speed up the convergence.
With the initial orientation parameters and the alignment parameters available, LGSP can be used to refine the RPCs of stereo images and provide the estimated accuracy for the RFM fitting solutions. The fitting accuracy of RFM is calculated by the remaining grid point, which are not involved in refining the RFM. Further, the feature matching points of the stereo images  are used to evaluate the alignment substitution accuracy of LGSP. Section 3.2 will discuss the definition, calculation procedure and results of both the fitting accuracy and the alignment substitution accuracy.
In addition, we compared the performance of LGSP with bundle adjustment in which four GCPs were deployed around the four corners of the stereo images. Bias compensation in bundle adjustment is the affine compensation in the image space to compensate the yaw angle and principal distance errors (Fraser and Hanley, 2005;Pan et al., 2018).

EXPERIMENTS
The inputs of LGSP are satellite stereo images and laser altimetry data. To test the performance of LGSP, the Ziyuan-3 (ZY-3) and Gaofen-7 (GF-7) stereo datasets are used as the experimental data. The ZY-3 satellite is China's first civilian high-resolution stereo mapping satellite, and the GF-7 satellite, which was newly launched on November 3, 2019, is China's first submeter high-resolution optical stereo-mapping satellite. The ZY-3 stereo datasets are located in a hilly area and contain forward and backward images with a 2.7 m ground sample distance (GSD) and the corresponding initial orientation parameters. The inclination angle of both forward and backward cameras is 22 degrees. The GF-7 stereo datasets are located in a mountainous area and contain forward and backward images with GSDs of 0.8 m and 0.65 m and the corresponding initial orientation parameters. The inclination angle of the forward and backward camera is 26 degrees and 5 degrees, respectively. The ZY-3 and GF-7 stereo datasets have 17 and 62 evenly distributed GCPs, respectively. These points were measured by GPS with an accuracy of approximately 0.1 m in both horizon and height and were used as either control points or checkpoints (CKPs) in the experiments. Laser altimetry data is the ICESat-2 ATL03 products close to the image capturing time. As shown in Figure 4, the ZY-3 and GF-7 stereo datasets have 1600460 and 307472 evenly distributed photons, respectively, which are superimposed onto the DSM generated by satellite stereo images. The above DSMs are drawn over the entire area, The above DSMs are plotted over the entire area for topographic display. In the actual processing, only the DSMs around the laser altimetry data are generated using the variant of tSGM.

Alignment of Laser Altimetry Data and DSM
After the laser altimetry data and the pre-alignment DSM point cloud were quickly aligned, the alignment transformation parameters were obtained using the ZY-3 stereo dataset. The alignment transformation parameters were applied to the DSM point cloud to obtain an aligned DSM based on the alignment parameters.
We first qualitatively compared the height changes in ICESat-2, the pre-alignment DSM, and the aligned DSM. Figure 5 shows the height profiles along the orbital direction of ICESat-2. Before alignment, ICESat-2 and the pre-alignment DSM have a remarkable offset in the height direction. After alignment, the profiles of ICESat-2 and the aligned DSM are matched, indicating that the alignment of ICESat-2 and the pre-alignment DSM in the height direction is necessary and valuable. However, there are many local inconsistencies between the aligned DSM and ICESat-2, mainly due to the limitation of image ground sampling distance, matching errors, and height changes of the surface. Additionally, many cusps appear in ICESat-2 profiles, which are mainly caused by the multiple echoes of the laser on the surface. The aforementioned inconsistencies are less in proportion and have been well removed by introducing farthest points removal in the alignment process.
We then quantitatively compared the height accuracy of the prealignment and aligned DSM. When the DSM was not aligned with the laser altimetry data, the root-mean-square error (RMSE) of the CKPs was 5.79 m in height. The RMSE was reduced to 0.85 m in height when DSM was aligned and transformed with the laser altimetry data.

Accuracy Verification of RFM Refinement Based on Alignment Parameters
After the alignment of the laser altimetry data with the DSM, the RFMs of the satellite stereo images were refined by the alignment parameters. To verify the effectiveness of the RFM refinement method based on alignment parameters, we used fitting accuracy and alignment substitution accuracy to evaluate the RFM of single images and stereo models composed of stereo images, respectively. The fitting accuracy was calculated from the grid points without refining the RFM, and the alignment substitution accuracy was calculated from the feature matching points.
In the process of recalculating the RPCs using virtual grids, grids checkpoints of the same order of magnitude were generated simultaneously and at equal intervals to verify the fitting accuracy of the RFM refinement. The 3D coordinates of the grid checkpoints were compared with the image square coordinates of the grid checkpoints by calculating the image square coordinates according to the updated RPCs, and the statistical information in Table 1. According to the statistical results in Table 1, the minimum value, maximum value, and RMSE of the difference are better than 1% of the image elements, which verifies that the refined RFM has high fitting accuracy.
For verifying the accuracy of the refined RFMs, the accuracy can be evaluated based on the 3D coordinates obtained from the tie points. The specific process is as follows: (1) the 3D ground coordinates of the tie points are obtained from the initial RFMs using triangulation, and then the new 3D coordinates are obtained by adding the alignment parameters; (2) the 3D ground coordinates of the tie points are obtained from the refined RFMs; and (3) the difference between the 3D ground coordinates of (1) and (2) is used, and the alignment substitution accuracy is shown in Table 1. The minimum value, maximum value, and RMSE of planemetry and height are below the centimeter level, which verifies the reliability of optimizing the high-resolution stereo model based on the alignment parameters.

Positioning Accuracy of High-resolution Satellite Stereo Images
In this experiment, the ZY-3 and GF-7 stereo datasets were carried out by LGSP. Thus, we also obtained the refined RFMs in addition to the initial RFMs. The height errors of the entire image region are illustrated in Figure 6. For the initial RFMs of the ZY-3 and GF-7 stereo datasets, there are large height errors at most CKPs. The height errors of all CKPs are almost equal for the GF-7 stereo datasets, but there are different error characteristics for the ZY-3 stereo datasets. After the RFMs were refined by LGSP, height errors of most CKPs are less than one meters for both the ZY-3 and GF-7 stereo datasets. Compared with ZY-3 stereo datasets, the height errors of GF-7 are smaller due to the difference in resolution resulting in a different level of detail in the DSM.
We then evaluated the bias and RMSE of the initial RFM, Bundle adjustment method with four GCPs and LGSP without GCPs. The bias and RMSE are calculated from CKPs that are not involved in the bundle adjustment. Table 2 and Table 3 show that the RFM of both the bundle adjustment method and LGSP that use the ZY-3 and GF-7 stereo datasets are superior to the initial RFM in both horizon and height. Numerically, On GF-7 stereo datasets, the bias of horizontal direction was reduced from 14.8 m to 1.24 m, and the RMSE was reduced from 14.81 m to 1.32 m. The RFM of both the bundle adjustment method and LGSP has higher positioning accuracy on the GF-7 stereo dataset than the ZY-3 stereo dataset, due to the GSD of GF-7 stereo images being smaller than that of the ZY-3 stereo images. More notably, LGSP reduced the height errors from the original 5.99 m to 1.00 m for the Ziyuan-3 stereo dataset and from the original 20.39 m to 0.83 m for the Gaofen-7 stereo dataset, and slightly outperformed bundle adjustment in which four ground control points were deployed around the four corners of the stereo images. Considering that LGSP does not require GCPs compared with the bundle adjustment method, LGSP is low-cost and efficient.

CONCLUSION
In this paper, we proposed a high-accuracy positioning method for high-resolution stereo images based on laser altimetry data, which consists of photon outlier removal, DSM generation by the variant of tSGM, fast point cloud alignment, DSM refining, and RFM refining. Experiments on two satellite stereo datasets with different loads, resolutions, and topographies verified that LGSP significantly improved the positioning accuracy of satellite images. The main work of this paper can be concluded as follows: 1. We proposed a high-resolution stereo model refinement method based on alignment parameters. This method has high fitting accuracy verified by grid checkpoints and high alignment substitution accuracy verified by feature matching. The experiments show that this method is essentially equivalent to the direct transformation of object coordinates based on alignment parameters.

Experiments on ZY-3 and GF-7 stereo datasets show that
LGSP significantly improved the height positioning accuracy of satellite images, and slightly outperformed bundle adjustment in which four GCPs were deployed around the four corners of the stereo images. When the plane error of the initial orientation parameters is large, LGSP also significantly improved the plane accuracy of satellite images.
This study confirms that LGSP is suitable for improving the positioning accuracy of high-resolution stereo images. With the increased availability of satellite stereo images, we plan to extend our method in further research to more types of satellite images, for example Pléiades and WorldView-3.