Registration and Combined Adjustment for the Laser Altimetry Data and High-Resolution Optical Stereo Images of the GF-7 Satellite

Abstract

The GF-7 satellite is China’s first civil sub-meter resolution stereo mapping satellite, aiming at 1:10,000-scale mapping. To achieve this goal, apart from the stereo optical cameras that reach sub-meter resolution, the GF-7 satellite is equipped with a laser altimetry system capable of obtaining three-dimensional laser points (LPs) with high elevation accuracy. However, the combination of laser altimetry data and optical stereo images has not been thoroughly studied. In this paper, we exploit the images recorded by the highly integrated laser footprint cameras and propose a hierarchical phase correlation method based on a geographic pyramid for the registration of laser altimetry data and high-resolution optical stereo images, which lays a solid foundation for the following combined adjustment. Experiments show that the proposed registration method can automatically locate the LPs on high-resolution stereo images and meet the requirements of bundle adjustment. A series of bundle adjustment experiments were carried out, showing that laser altimetry data can significantly enhance the vertical accuracy of optical image stereo mapping and that elevation accuracy can reach roughly 1.0 m (RSME) without ground control points. Therefore, this study could be a good guide for global high-precision DSM acquisition with the GF-7 satellite.

1. Introduction

The GF-7 satellite is China’s first civilian optical stereo mapping satellite with sub-meter resolution [1]. Apart from two line-array stereo optical imaging cameras, the GF-7 satellite carries a laser altimetry system, which is China’s first full-waveform laser altimeter for Earth observation. The principal goal of the altimetry system is to obtain sparse elevation control points and ultimately improve the stereo mapping accuracy of high-resolution optical images without ground control points. To achieve this goal, the raw data from the GF-7 satellite laser altimetry system is processed into standard data products [2] that are as accurate as possible, so that they are convenient for application in stereo mapping as well as other fields.

In previous research practices, it has been proved that the vertical accuracy of stereo images can be improved with elevation control points obtained by satellite laser altimetry technology. High elevation accuracy of optical satellite stereo mapping requires demanding attitude measurements [3], so taking advantage of laser altimetry should be a good alternative. The idea of laser ranging-assisted stereo mapping has been proposed, and simulations have been carried out for validation [3,4,5] which lay the theoretical foundation for laser altimetry technology-assisted optical image stereo mapping. The ICESat (Ice, Clouds, and Land Elevation Satellite) [6], launched by NASA in 2003, provided the data to validate these theories. ICESat/GLAS laser altimetry data are used as a reference for the elevation optimization of DSM produced by ALOS PRISM stereo imagery [7]. High-accuracy ICESat laser altimetry data are extracted and introduced as elevation control points into a bundle adjustment of optical stereo images from the ZY-3 satellite or Mapping Satellite-1, which considerably improved the stereo mapping accuracy [8,9,10]. These works have proved the effectiveness and feasibility of satellite laser altimetry-assisted stereo mapping of optical remote sensing images. In May of 2016, an experimental laser altimeter for the purpose of elevation control was launched on-board the ZY3-02 satellite [11]. Without waveform-collecting ability, the data process was relatively simple and could only provide three-dimensional points as products. Nevertheless, it was enough to demonstrate the fact that laser altimetry can help improve the elevation accuracy of stereo mapping [12,13,14].

The integration of laser altimetry data and optical images is also used in planetary remote sensing. Surface matching with a 3D rigid transformation model is used to make the laser altimetry data and DEM generated by stereo images consistent [15], which is only valid when the surveying area is large enough and when there is plenty of laser altimetry data. A combined block adjustment is proposed to integrate multiple strips of the Chinese Chang’E-2 imagery and NASA’s Lunar Orbiter Laser Altimeter (LOLA) data for topographic modeling [16]. A local surface constraint is developed to incorporate the LOLA points into the block adjustment since the exact image coordinates of the Lunar Orbiter Laser Altimeter (LOLA) points cannot be determined in the stereo images. However, the terrain features on Earth are far more complex than those on the Moon. As a result, this method is inapplicable to earth observation.

For more demanding mapping tasks, the accurate registration between laser altimetry data and optical stereo images should be taken into consideration. To further realize the potential of laser altimetry technology in assisting optical stereo mapping, the GF-7 satellite, launched in November 2019, carried a more powerful laser altimetry system [1]. It has two laser beams, each integrated with a laser footprint camera. The laser footprint cameras share the optical telescope and work closely with the laser instruments. Therefore, it gives a chance for strict registration between laser altimetry data and optical remote sensing images, which is impossible for the former ICESat/GLAS or ZY3-02 satellite. There is, however, limited commercial software that directly supports this kind of data; additional specific developments should be completed for technology adaptation.

This study proposes a strict registration method for laser altimetry data and optical stereo images according to the characteristics of the GF-7 satellite, and assesses the combined stereo mapping of optical images and laser altimetry data. In Section 2, we first elaborate on the laser altimetry system of the GF-7 satellite and its data products. Then, we propose a hierarchical phase correlation method based on geographic pyramids for the registration of laser altimetry data and optical stereo images and design the experimental scheme. In Section 3, the experiments are presented with the results, which we discuss in Section 4. Finally, conclusions and prospects are provided.

2. Materials and Methods

2.1. The GF-7 Satellite Data

The GF-7 satellite laser altimetry system is comprised of two laser instruments, each equipped with a footprint camera. The lasers fire three times per second at an altitude of about 500 km. As a result, the distance between adjacent ground laser spots within a beam is about 2.4 km. Each laser instrument points at 0.7 degrees across the direction of flight, and the ground tracks are about 12.25 km apart. The diameters of the ground laser footprint are both about 20 m. The laser footprint camera, unlike the linear array high-resolution optical camera, employs a CMOS (Complementary Metal Oxide Semiconductor) area-array imaging sensor to capture a two-dimensional image array for each exposure. The imaging spectrum of the footprint camera ranges from 500 nm to 720 nm and also 1064 nm, which is specified for recording the output laser energy. Part of the output laser pulse is introduced into the laser footprint camera to form a laser spot on the laser footprint image. The pixel size of the CMOS sensor is 5.5 µm. However, a 3 × 3 pixel-merging is conducted during imaging; thus, the equivalent pixel size is 16.5 µm. As a result, the ground resolution of the footprint image is reduced to 3.2 m. Only a 550 × 550-pixel window of the entire sensor array is preserved, resulting in a ground range of 1.76 km × 1.76 km for the laser footprint images. Figure 1 shows the working scheme and the simplified geometry relationship for the GF-7 laser altimetry system [17].

After the raw laser altimetry data is downlinked, a series of steps are performed to convert it into standard data products, including waveform decomposition and feature extraction, the rigorous geolocation of laser footprint images, ranging correction for atmospheric delay, height corrections for tides, elevation control flag determination, and so on [2]. Apart from the 3D coordinates of laser ground spots, which are referred to as “laser points” (LPs), the standard data products provide waveform data, laser footprint images, and various feature parameters or corrections for accuracy refinement [18]. On-orbit calibrations are carried out for both the laser altimeters and the laser footprint cameras [17,19]. The planimetric geolocation accuracy of the ground laser points is within 5.0 m, and the elevation accuracy is better than 0.1 m in the calibration region [2]. According to the waveform features and other prior information such as global land cover, each laser point is evaluated to determine whether it can be used for elevation control [9]. The rigorous geometric model of the laser footprint cameras is built [17], based on which the Rational Polynomial Coefficients (RPCs) [20] for laser footprint images are calculated and provided in the standard data products. The laser instrument and the corresponding footprint camera are highly integrated and have a stable geometrical relationship. In addition, the output laser spot on the laser footprint image can be used to detect whether the relationship has changed [21]. Therefore, the position of the laser point in the corresponding footprint image can be determined directly.

The dual-line-array stereo optical cameras onboard the GF-7 satellite acquire sub-meter resolution stereo images, including forward and backward images (FWD and BWD). The forward and backward view angles are 21° and −5°, respectively. The ground resolution of the stereo images is 0.64 m for BWD and 0.8 m for FWD. There is also a multispectral camera that acquires 4-band multispectral images at a resolution of 2.56 m. The laser altimeter and stereo cameras can be operated simultaneously or individually. The simultaneously obtained data distribution is shown in Figure 2. In the range of a single high-resolution image scene, there are 2 columns, each containing up to 11 LPs. Satellite laser altimetry and optical remote sensing technology have their own advantages and are complementary to each other. The integration of this data is essential for high-accuracy stereo mapping.

2.2. Registration for Laser Altimetry Data and Stereo Images

To introduce laser altimetry data into stereo mapping, the LPs must be automatically located in the stereo images. However, directly projecting the laser spot onto stereo images may result in significant errors, because the projected positions on different stereo images are not homologous points due to the geolocation bias of the stereo image. As previously mentioned, the geometry relationship between the laser footprint camera and the laser instrument is robust, allowing the laser altimetry data to be rigorously aligned with the stereo images. This can be accomplished by matching a fixed point in the footprint image that corresponds to the laser ground spot to the high-resolution optical stereo images. However, registration between the footprint and stereo images is no easy task. Apart from the extreme difference in resolution, there are also obvious radiant differences between laser footprint images and stereo images from different types of imaging sensors. Such differences make registration unachievable for matching methods utilizing intensity- or gradient-based image features, such as SIFT [22], SURF [23], etc.

In this paper, we propose a hierarchy phase correlation method for the registration of altimetry data and optical stereo images. Taking the elevation of the LP as reference, the footprint image and the stereo image are both rectified based on their corresponding RPCs. As a result, the rotation and scaling of the original footprint image and stereo images are eliminated. Then, a sub-pixel registration based on a Fourier transform, i.e., phase correlation, is applied to the rectified image chips to determine the pixel shift.

Phase correlation has been described in previous work [24], but will also be briefly summarized here. Suppose there are two images, i.e., $f_a$, $f_b$, with a shift of $\left(x_0,y_0\right)$; their relationship can be represented as:

$$f_a\left(x,y\right)=f_b\left(x-x_0,y-y_0\right)$$

(1)

According to the shift theorem, taking the Fourier transform of both images, the relationship is rewritten as:

$$F_b\left(u,v\right)=F_a\left(u,v\right)e^{-j\left(u{x}_0+v{y}_0\right)}$$

(2)

where $F_a$ and $F_b$ are the Fourier transform of $f_a$,$f_b$, respectively. Then, the cross-power spectrum of two images is:

$$\frac{F_a\left(u,v\right)F_b{}^{\ast }\left(u,v\right)}{\left|F_a\left(u,v\right)F_b{}^{\ast }\left(u,v\right)\right|}=e^{-j\left(u{x}_0+v{y}_0\right)}$$

(3)

where $F^{\ast }$ is a complex conjugate of $F$. By inverse Fourier transform, we can obtain a pulse function whose peak coordinates are located on the offset $\left(x_0,y_0\right)$. In practice, the amount of translation is not necessarily an integer. Here, we take the advantage of sub-pixel peak location, which is approximated by a simple formula involving the peak pixel value and the values of its nearest neighbors [25].

$$\left\{\begin{array}{c}dx=\frac{v_{\left(1,0\right)}}{v_{\left(1,0\right)}\pm v_{\left(0,0\right)}}\\ dx=\frac{v_{\left(0,1\right)}}{v_{\left(0,1\right)}\pm v_{\left(0,0\right)}}\end{array}\right.$$

(4)

where $dx,dy$ are the sub-pixel shifts, $v_{\left(0,0\right)}$ represents the value at the peak of pulse function, and $v_{\left(0,0\right)}$, $v_{\left(0,0\right)}$ are the direct neighbors in the X- and Y-directions. Combining the pixel shift and the rectifying resolution, the ground bias between the footprint image and the stereo image is obtained.

$$\left\{\begin{array}{c}dLon=\left(x_0+dx\right)\ast gsd\\ dLat=\left(y_0+dy\right)\ast gsd\end{array}\right.$$

(5)

where $dLon,dLat$ are the geolocation biases of the laser footprint image and the stereo image; $gsd$ is the rectifying resolution, which can be determined according to the resolution of the laser footprint image and the stereo image. Then, the position of the LP in the stereo image can be calculated with the RPCs of the stereo image:

$$\left\{\begin{array}{c}{x}_l=f\left(lon+dlon,lat+dlat,alt\right)\\ y_l=g\left(lon+dlon,lat+dlat,alt\right)\end{array}\right.$$

(6)

where $x_l,y_l$ is the matched point in the stereo image of the LP, $lon,lat,alt$ are the 3D coordinates of the LP, which is given in standard data products of the GF-7 satellite laser altimetry system, $f$ and $g$ denote the rational polynomial functions in sample and line directions [20], respectively.

Moreover, our registration approach utilizes a coarse-to-fine pyramid strategy for optimization, with a coarse offset recognized at the top of the pyramid and a fine offset detected at the bottom. However, different from the traditional image pyramid, which changes the size of images, our new image pyramid keeps the image size constant and gradually upgrades the ground resolution. We refer to this new type of pyramid as a “geographic pyramid”, or “geo-pyramid” for short. The advantage of the geo-pyramid is that the geographic extent corresponding to each layer is reduced, and we can obtain the bias of the original images from coarse to fine. On the contrary, in traditional image pyramids, the geographic extents are the same for each layer, making it difficult to determinate the initial size of the pyramid for the reason that a large size is not conducive to detecting local deviations, while a small size would miss the corresponding area if the geolocation bias of the stereo images is large enough. Our geographic pyramid could handle this tough issue to a great extent, as our initial geographic range is a dozen times greater than that of the traditional pyramid, while the final resolutions are comparable (Figure 3). There is also another superiority of our geo-pyramid. The Fourier transform is computationally expensive. When the image sizes in a conventional pyramid increase, the computational efficiency drops quickly. In our method, however, the image sizes are maintained and can, thus, achieve a higher computational efficiency.

From the second layer (top to bottom) of the geo-pyramid, the reference point of the stereo images is updated; that is, the ground bias obtained on the lower layer is added up to the laser point. By iterative rectification and phase correlation until the bottom layer, an accurate registration of the laser footprint image and the high-resolution image is realized, which is shown in Figure 4. Therefore, the laser point is precisely located in the high-resolution optical image.

Mismatches are unavoidable during the registration process; two methods are adopted to identify them. First, the Mean Structural Similarity Index (MSSIM) [26] is calculated between image chips, and those with MSSIM values under certain thresholds are viewed as mismatches. Then, the ground coordinates are calculated through triangulation and are compared to the coordinates of the corresponding laser point. Note that the resolution and radiant quality of the stereo images is much better than the footprint images. Applying registration for a footprint image to all overlapped stereo images may cause discordance because there may be local deformations in stereo images, which are indistinguishable from the footprint images. Therefore, the laser footprint images are just used as a bridge to locate the ground laser spots in one of the stereo images. For the positions in other stereo images, registration is carried out between stereo images using the proposed method.

2.3. Bundle Adjustment for Stereo Mapping

For large-scale stereo mapping, bundle adjustment is crucial to improving the final accuracy. By bundle adjustment, the discrepancies among multiple stereo images are eliminated, and global geolocation bias can also be removed if the ground control points are provided.

The GF-7 satellite’s stereo image products come with an RPC file for georeferencing, as most imaging providers do. The RPC model [20] can establish a direct relationship between an image point and the corresponding object point independently of the sensors, as follows:

$$\left\{\begin{array}{c}l=\frac{NumL\left(X,Y,Z\right)}{DenL\left(X,Y,Z\right)}\\ s=\frac{NumS\left(X,Y,Z\right)}{DenS\left(X,Y,Z\right)}\end{array}\right.$$

(7)

where $l$ and $s$ are the line and sample in an image coordinate system, representing the rows and columns of the point; functions $NumL$, $DenL$, $NumS$, $DenS$ are the cubic polynomial functions composed by the 80 RPCs; $X,Y,Z$ are the normalized object coordinates by the other 10 RPCs.

To carry out the bundle adjustment, a polynomial function is selected for the bias compensation in image space.

$$\left\{\begin{array}{c}\Delta l=a_0+a_{1}l+a_{2}s+\dots \\ \Delta s=b_0+b_{1}l+b_{2}s+\dots \end{array}\right.$$

(8)

where $a_i$, $b_i$ are the bias-compensation parameters, and $\Delta l$, $\Delta s$ are the compensating values. Here, a second-order polynomial function is enough to compensate for the offsets, drift, or other small systematic errors [27,28,29].

The bundle adjustment is to compute the bias-compensation parameters by globally optimizing the reprojected errors in image space and also object coordinate errors if control data exists [30]. According to the type of interest points, the optimization problem can be formulated as:

$$\begin{array}{c}minE\left(x,y,Lon,Lat,Alt\right)={\sum }_i{\sum }_j\parallel RPC\left(Lo{n}_j,La{t}_j,Al{t}_j\right)-{\left(x_i,y_i\right)}^T\parallel \\ +{\sum }_{j\in GCP}\parallel {\left(Lo{n}_j,La{t}_j,Al{t}_j\right)}^T-{\left(Lo{n}_j^g,La{t}_j^g,Al{t}_j^g\right)}^T\parallel \\ +{\sum }_{j\in HCP}\parallel {\left(Lo{n}_j,La{t}_j\right)}^T-{\left(Lo{n}_j^h,La{t}_j^h\right)}^T\parallel \\ +{\sum }_{j\in VCP}\parallel Al{t}_j-Al{t}_j^v\parallel \end{array}$$

(9)

where the $\left(x_i,y_i\right)$ are images coordinates, and $\left(Lo{n}_j,La{t}_j,Al{t}_j\right)$ are the corresponding object coordinates. $\left(Lo{n}_j^g,La{t}_j^g,Al{t}_j^g\right)$ are 3D coordinates of ground control points (GCP), $\left(Lo{n}_j^h,La{t}_j^h\right)$ are the planar coordinates of the horizontal control points (HCP), and $Al{t}_j$ is the elevation of the vertical control points (VCP). Thus, the bundle adjustment is a least-squares optimization problem, which can be efficiently solved using the Ceres Solver [31].

Here, we design eight different kinds of experiments to comprehensively study the stereo mapping results, depending on whether the LPs or GCPs are included. They are as follows:

• Plan A: free network adjustment;
• Plan B: adjustment with all GCPs that are available;
• Plan C: adjustment with a modest number of GCPs;
• Plan D: adjustment with LPs as control;
• Plan E: adjustment with a few GCPs as horizontal control;
• Plan F: adjustment with a few GCPs as vertical control;
• Plan G: adjustment with all LPs as vertical control;
• Plan H: adjustment with a few GCPs as horizontal control and LPs as vertical control.

In these experiments, the LPs are located in each stereo image using the proposed method. The GCPs are manually measured in the stereo images, and those not used in the bundle adjustment are viewed as check points (CPs) to evaluate the results. It should be noted that the tie points are extracted and refined in advance and remain identical in the experiments to keep the results reliable for comparison. The overall flow chart of this study is shown in Figure 5.

3. Results

3.1. Study Area and Data

We chose two sites on which to carry out the studies. One is in China’s province of Shaanxi, and the other is on the border of Laos and Thailand. In the Shaanxi study area, we chose 4 GF-7 satellite stereo image pairs and a set of laser altimetry data, all of which were acquired in two orbits, on 26 April and 1 May 2020, respectively. Additionally, a total of 20 ground control points (GCPs) were collected. For the laser altimetry data, there are 79 LPs, of which 26 points were used in this area. In the Laos–Thailand border area, eight stereo image pairs from two orbits, on 24 February and 1 March 2020, were selected. The laser altimetry data is available only on February 24. Excluding the points with low precision, there are 53 LPs remaining for use. Due to the fact that this experimental region is located outside of China, only a limited number of GCPs are available. Figure 6 illustrates the distribution of the experimental data.

3.2. Results of Registration for Laser Altimetry Data and Stereo Images

Excluding the LPs with poor accuracy, 26 and 53 points are used for registration in the Shaanxi area and the Laos–Thailand border area, of which 24 and 49 points are successful, respectively. The success rates both exceed 92%.

An example is shown in Figure 7. Note that the footprint image has been cropped and magnified four times for demonstration purposes. The stereo images are presented at the original resolution. The laser ground spot is located in a weak-textured region. This is often the case because of the determination strategy of elevation control points [9]. As a result, matching the footprint image with the stereo images is more difficult than extracting tie points within stereo images, which is often achieved using features such as SIFT, SURF, and so on. The trail from north to south in the footprint image is undetectable because of poor resolution.

In another case, shown in Figure 8, there is a sharp radiant distinction between the laser footprint image and the stereo images. Traditional image similarity measurements, such as correlation coefficients, would produce insufficient results. In addition, ground objects or terrain would vary on images from different periods, which would also lead to poor results for classical image correlation. However, our proposed method is robust for radiant distinction and ground variation and succeeds in registration.

To assess the accuracy of matched points in stereo images, we treat them as tie points and conduct a free network bundle adjustment. The residuals are presented in

Table 1. Residuals of matched LPs in stereo images.

Experiments	Selected	Matched	Success Rate	X Res/Pixels	Y Res/Pixels
Shaanxi Area	26	24	92.3%	0.31	0.40
Laos–Thailand border Area	53	49	92.5%	0.19	0.10

. We can see that residuals in the two study areas are both less than 0.5 pixels, with the Laos–Thailand border area being slightly better than the Shanxi area. It is, therefore, clear that the results obtained with our method are resilient and sufficient for subsequent bundle adjustments.

3.3. Results of Bundle Adjustment

In the Shaanxi area, there are enough GCPs that can be used. So, when GCPs are needed in the experiment of Plan C, E, F, and H, five GCPs that are distributed at the four corners and in the center of the survey area are used, while other GCPs are evenly distributed around the area. In the Laos–Thailand border area, since the GCPs are very sparse, we select two that are as far away as possible for ground control. Though the amount and distribution of the GCPs are not ideal in this area, it still serves as a good example of the problem.

The root mean square errors (RMSE) and mean errors (ME) of all the experiments in the two study areas are presented in

Table 2. Experimental results in the Shaanxi area.

Index	GCPs	CPs	LPs	X/m		Y/m		Horizontal/m		Elevation/m
				RMSE	ME	RMSE	ME	RMSE	ME	RMSE	ME
Plan A	0	20	0	12.758	−12.703	6.393	−6.135	14.270	14.209	4.486	4.238
Plan B	20	0	0	0.827	−0.122	0.903	−0.064	1.224	1.065	0.481	0.154
Plan C	5	15	0	1.272	−0.559	1.313	−0.442	1.828	1.564	0.737	0.175
Plan D	0	20	24	5.686	−4.784	2.905	1.283	6.385	6.169	1.115	−0.388
Plan E	5H	15	0	1.452	−0.955	2.471	−1.391	2.866	2.513	5.368	5.267
Plan F	5V	15	0	12.203	−12.159	5.995	−5.866	13.596	13.557	0.781	0.217
Plan G	0	20	24V	10.675	−10.527	4.792	−4.631	11.701	11.558	1.192	−0.379
Plan H	5H	15	24V	1.586	−1.012	1.445	−0.598	2.146	1.907	0.934	−0.148

and

Table 3. Experimental results in the Laos–Thailand border area.

Index	GCPs	CPs	LPs	X/m		Y/m		Horizontal/m		Elevation/m
				RMSE	ME	RMSE	ME	RMSE	ME	RMSE	ME
Plan A	0	5	0	11.809	−11.781	5.033	−4.838	12.837	12.800	3.363	−3.294
Plan B	5	0	0	0.914	−0.486	0.835	−0.114	1.238	1.147	0.366	0.059
Plan C	2	3	0	1.630	−1.395	1.061	0.448	1.945	1.738	1.299	0.846
Plan D	0	5	49	2.475	−2.398	5.362	5.105	5.906	5.692	0.937	−0.612
Plan E	2H	3	0	2.235	−2.026	2.596	−0.843	3.405	3.174	7.273	7.197
Plan F	2V	3	0	11.024	−11.002	5.291	−5.189	12.228	12.221	1.217	0.796
Plan G	0	5	49V	6.670	−6.577	1.193	−0.704	6.776	6.659	1.037	−0.794
Plan H	2H	3	49V	1.215	−1.024	1.334	0.868	1.804	1.586	0.794	−0.020

. The accuracy assessments for check points or residuals of GCPs in the Shaanxi area are shown in Figure 9. The original horizontal accuracy of GF-7 stereo images is about 12~15 m without GCPs, while elevation accuracy ranges between 3~5 m. The original geolocation errors are to the southwest, whereas the vertical errors are all upward in the Shaanxi area, which can be seen in Figure 9a. In the Laos–Thailand border area, the horizontal errors are consistent with those in the Shaanxi area, but the vertical errors are downward. This suggests that the GF-7 stereo images contain systematic geolocation inaccuracies. The horizontal errors are comparable in both areas, while the vertical errors are diametrically opposed. When all the GCPs are included for control, the residuals are both around 1.2 m in the horizontal direction and better than 0.5 m in the vertical direction. It is obvious that, when using GCPs, the systematic errors are eliminated, seeing Figure 9b, and that the magnitude and orientation of the geolocation errors turn stochastic. The accuracy is improved remarkably with a few GCPs. The horizontal accuracy is less than 2.0 m and the elevation accuracy reaches within 1.3 m. In comparison to Plan B, it is evident that more GCPs will yield a better result, but the improvement is not that significant, as in Figure 9c.

Thanks to the laser altimetry data, without GCPs, the horizontal accuracy improved from worse than 10 m to roughly 6.0 m, and from over 3.3 m to better than 1.2 m in the vertical direction. Figure 9d shows that the systematic vertical errors have been eradicated, but the systematic horizontal errors remain because of the geolocation inaccuracy of the laser altimetry data. If the GCPs serve as horizontal controls, the horizontal accuracy is improved dramatically to about 3.0 m, while the vertical accuracy systematically declines to as low as 5~7 m, which can be seen in Figure 9e. Just the opposite happens if the GCPs are employed as vertical controls. The vertical accuracy is improved dramatically to within 1.3 m, while the horizontal accuracy remains systematic and can be as large as more than 10 m, as in Figure 9f. With all the LPs as vertical controls, the vertical accuracy is improved to about 1.0 m. This is comparable to Plan C, in which a few GCPs are used. However, the improvement of horizontal accuracy is not obvious. We can still see systematic errors in Figure 9g. In the Laos–Thailand border area, the horizontal accuracy gains some improvement, but in the Shaanxi area, it turns out to be just slightly better. For the last plan, in which the GCPs and the LPs are treated as horizontal and vertical controls, respectively, the horizontal accuracy approaches about 2.0 m, while the vertical accuracy reaches within 1.0 m, which can be comparable with the results obtained with a few GCPs, as in Figure 9h.

Apart from the limited GCPs, we also evaluate the results of the combined bundle adjustment using ICESat2 altimetry data [32,33] in the Laos–Thailand border area. The DSM is produced [34] after bundle adjustment with only LPs as control, that is, Plan D. The ATL08 products of ICESat2 within the DSM range are collected and extracted for the evaluation. The elevation accuracy of the extracted ATLAS LPs can reach 0.429 m and 0.596 m in flat and hilly areas, respectively [35]. A total number of 11,784 points, as seen in Figure 10, are available.

Table 4. Height discrepancies of DSM and extracted ATLAS LPs.

	Total	|dH| < 1.0	|dH| < 1.5	|dH| < 3.0	|dH| < 6.0
Count	11,784	5883	7948	10,721	11,587
Percentage	100%	49.9%	67.4%	91.0%	98.3%
Average/m	0.465	0.193	0.331	0.566	0.595
Standard deviation/m	2.029	0.521	0.724	1.136	1.530

shows the statistics of height discrepancies, and Figure 11 shows the discrepancy histogram, which is almost regularly distributed. Half of the extracted ATLAS LPs have a height discrepancy of less than 1.0 m and more than 90% within 3.0 m, compared to the DSM generated from stereo images and LPs without GCPs.

The extracted ATLAS LPs with different height discrepancies are shown in Figure 12. The ATLAS LPs with smaller discrepancies are dense and evenly distributed, whereas those with larger discrepancies are distributed randomly and sparsely throughout the area. The fact that the datasets were acquired at different times is the most likely cause of the large disparities. In other words, the DSM generated with stereo images and LPs from the GF7 satellite in the absence of GCPs agrees well with ICESat2 laser altimetry data. Note that, although the LPs for combined adjustment are concentrated on the western half of the survey area (Figure 6b), there is no significant difference in the height discrepancies across the survey area.

4. Discussion

To apply the high elevation accuracy LPs to the high-resolution optical image stereo mapping, it is critical to precisely determine the position of the LP on each stereo image. The highly integrated laser footprint cameras have laid a good foundation for achieving this goal. However, imaging sensors of different types and performances bring difficulties to image registration. The method proposed in this paper solves this problem well and fully brings out the function of the laser footprint camera, thereby making it possible for the automated combined stereo mapping of laser altimetry data and high-resolution optical images. This would be a good supplement to existing commercial software.

In the absence of ground control points, the horizontal positioning accuracy of GF-7 satellite stereo images is in the order of about 10~15 m, while the elevation accuracy is approximately 3~5 m. The fact that the stereo images used in this study were generated with real-time attitude and orbit data on the satellite contributes to the majority of the errors. The geolocation errors are irregular on the whole, but appear obviously systematic in local areas. The use of a modest number of ground control points can eliminate the majority of systematic bias and significantly enhance geometric positioning accuracy. The horizontal accuracy can reach within 2.0 m, and the vertical accuracy is up to 1.3 m.

After using the laser altimetry data, the accuracy of the stereo mapping has improved significantly, especially the elevation accuracy, which can reach a level of roughly 1 m. This is due to the inherent advantage of laser altimetry technology in terms of elevation accuracy. Meanwhile, the laser altimetry data is generated on the basis of post-processed precise orbit and attitude data [36,37] other than what are obtained in real time on the satellite. As a result, the accuracy of combined adjustment with high resolution stereo images and laser altimetry data is improved in both the horizontal and vertical direction. The conclusion has been validated by independent data sources, such as ground check points or ICESat2 laser altimetry data.

After incorporating laser altimetry data from the GF-7 satellite into optical image stereo mapping, horizontal positioning error is the factor that limits the accuracy of uncontrolled stereo mapping. This is because the horizontal positioning precision of the laser altimetry data is insufficient (about 6 m), which is consistent with the accuracy verification results of the laser altimetry data products. Therefore, it is recommended to use laser altimetry data only as elevation control and in conjunction with several high-precision horizontal control points. High-precision elevation points are difficult to collect and usually necessitate a field survey. In comparison, high-precision horizontal control points, which can be obtained from existing reference images or open-source maps, are relatively easy to acquire. In this way, the accuracy of the stereo mapping of the GF-7 satellite can basically reach the level of using GCPs.

5. Conclusions

Satellite laser altimetry technology has unique advantages over optical remote sensing in elevation measurements. Apart from the three-dimensional structure of forests or the detection of changes in ice elevation, satellite laser altimeters can also be used to obtain global elevation control points for optical stereo image mapping. In this paper, we studied the combined processing of GF-7 laser altimetry data and high-resolution optical stereo images. A phase correlation registration approach based on geographic pyramids is proposed, which effectively solves the key issue of the automatic positioning of LPs on stereo images and lays a solid foundation for subsequent combined bundle adjustments. Based on the registration, the LPs are then brought into bundle adjustment for stereo images. A series of experiments were carried out showing that laser altimetry data can effectively enhance the vertical accuracy of optical image stereo mapping, and elevation accuracy can reach roughly 1.0 m. This mode of combined stereo mapping with laser altimetry data and stereo images will effectively relieve the burden of field mapping, and make high-accuracy mapping possible for overseas or inaccessible regions.

The operation of the GF-7 satellite will promote 1:10,000-scale topographic mapping and realize an all-around breakthrough in the field of civil large-scale mapping satellites in China. With the help of the laser altimetry system, the stereo mapping accuracy can be significantly improved, which is especially important in overseas territories or regions that are inaccessible. However, due to the selection strategy of laser altimetry data as elevation control points, a considerable amount of laser altimetry data is left unused. How to devise a universal bundle adjustment scheme and make full use of the laser altimetry data in stereo mapping will be further studied in the future.