Integrating Stereo Images and Laser Altimeter Data of the ZY3-02 Satellite for Improved Earth Topographic Modeling

Abstract

The positioning accuracy is critical for satellite-based topographic modeling in cases of exterior orientation parameters with high uncertainty and scarce ground control data. The integration of multi-sensor data can help to ensure precision topographical modeling in such situations. Presently, research on the combined processing of optical camera images and laser altimeter data has focused on planetary observations, especially on the Moon and Mars. This study presents an endeavor to establish a combined adjustment model with one constraint in image space for integration of ZY3-02 stereo images and laser altimeter data for improved Earth topographic modeling. The geometric models for stereo images and laser altimeter data were built first, and then, the laser ranging information was introduced to construct a combined adjustment model on the basis of the block adjustment model. One constraint that minimized the back-projection discrepancies in image space was incorporated into the combined adjustment. Datasets in several areas were collected as experimental data for the validation work. Experimental results demonstrated that the inconsistencies between stereo images and laser altimeter data for the ZY3-02 satellite can be reduced, and the elevation accuracy of stereo images can be significantly improved after applying the proposed combined adjustment. Experiments further proved that the improved height accuracy is insensitive to the number and relative position of laser altimeter points (LAPs) in stereo images. Moreover, additional plane control points (PCPs) were incorporated to achieve better planimetric accuracy. Experimental results in the Dengfeng area showed that the adjustment results derived by using LAPs and additional four PCPs were only slightly lower than those for the block adjustment with four ground control points (GCPs). Generally, the proposed approach can effectively improve the quality of Earth topographic model.

1. Introduction

Surveying and mapping play a very important role in various socioeconomic and national defense related construction projects. Especially, in regard to the overall issues and strategic decision making, most of these projects rely on the geospatial information obtained by surveying and mapping. Since space-based photogrammetry technology allows for fast updates and low-cost data acquisition that is not hindered by regional and national boundaries, countries around the world have been eager to develop their own surveying and mapping satellite programs. Geometric positioning accuracy is a key index to measure the mapping satellite performance and is a critical factor that constrains the accuracy of satellite image mapping data.

The geometrical positioning accuracy of surveying and mapping satellite is mainly limited by the attitude performance of the satellite platform which affects the precision of orientation parameters [1]. Therefore, it is possible to achieve high geometric positioning accuracy by acquiring high-precision orientation parameters via the high-accuracy attitude determinations [2] or by using a large number of ground control points to correct positioning errors [3]. However, it is particularly paramount to ensure the positioning accuracy in cases of exterior orientation parameters with high uncertainty and scarce ground control data.

As a high-precision ranging instrument, laser altimeters are widely used in the field of aerospace photogrammetry to obtain high-precision vertical information because of their good directivity and high ranging accuracy. Space-borne laser altimetry systems such as the Geoscience Laser Altimetry System (GLAS) system [4,5], the Mars Orbiter Laser Altimeter (MOLA) system for the observation of Mars [6,7], the NEAR Laser Rangefinder (NLR) system for the observation of space asteroids [8], the Clementine system [9], as well as China’s ChangE-1 laser altimeter system [10,11] for lunar observations have been found to be useful in the exploration of celestial bodies such as the Earth, Moon, and Mars; these systems are the ones most commonly used to generate topographic models of celestial bodies.

In previous research involving optical stereoscopic mapping, compared with the horizontal accuracy, the vertical accuracy was found to be more difficult to achieve because of the influence of the base-to-height ratio, platform stability, and other factors [12]. Consequently, integrating the stereo imagery and laser altimeter data has the potential to generate better accuracy, thereby making up for the defects such as poor attitude measurements and difficulty in obtaining ground control data.

At present, much of the work on the combined processing of optical camera images and laser altimeter data so far has focused on planetary observations, especially on the Moon and Mars.

In applications of Mars topographic mapping, the combined processing of optical imagery and laser altimeter data has been proposed. Anderson and Parker [13] examined the precision registration between Mars Orbiter Camera (MOC) imagery and MOLA data at selected candidate landing sites. The integration of MOC imagery and MOLA data was further studied by Yoon and Shan [14]. They introduced combined processing between the Mars Orbiter Camera (MOC) imagery and Mars Orbiter Laser Altimeter (MOLA) data using an adjustment method to correct the mis-registration and accurately determine ground positions. The results indicated that the large mis-registration between the two datasets can be corrected to a certain extent by the combined adjustment. Spiegel [15] and Ebner et al. [16] developed a high-resolution stereo camera (HRSC) imagery bundle adjustment technique, in which a sparse stereo point cloud was adjusted to optimize its fit to a surface interpolated from the MOLA data. The results showed the potential of the image matching and bundle adjustment approaches for achieving improved exterior and interior orientations with the MOLA digital terrain model (DTM) as control information. However, these methods failed to take into account the errors induced by MOLA data.

The integration of optical imagery and laser altimeter data also has been investigated in the applications of lunar topographic mapping. Rosiek et al. [17] presented an endeavor to integrate the Clementine images with the Clementine laser altimeter data, in which the Clementine global mosaic was used to establish horizontal control and Clementine laser altimeter points (LAPs) were used for vertical control. Di et al. [18] presented an integration method involving a lunar digital elevation model (DEM) derived from the Chang’E-1 stereo images and the laser altimeter data to reduce the inconsistencies between them. In their method, a DEM was generated first from the stereo images, and then the DEM was registered to the laser altimeter data through surface matching using an ICP (iterative closest point) algorithm [19], after which the exterior orientation (EO) parameters of the images were adjusted so that the inconsistency between the charged-couple device (CCD) images and the laser altimeter data was significantly reduced by this co-registration. In order to reduce the inconsistencies between adjacent orbits, Di et al. [20] further improved the co-registration of Chang’E-1 stereo images and Chang’E-1 laser altimeter data by incorporating a crossover adjustment of the laser altimeter data and refinement of the CCD image sensor model. In their method, crossover adjustment was employed to reduce the inconsistencies of different laser altimeter profiles, and this yielded more accurate and consistent control data; refinement of the image sensor model was realized by adding attitude angle bias corrections through a least-squares adjustment, from which consistency between the refined DEM from stereo imagery and LAPs is improved. However, both of the two methods were fixed in terms of the registration process and failed to consider the errors induced by laser altimeter data, which bring about the accuracy of the final generated topographic models; these models are totally dependent on the accuracy of the laser altimeter data.

In addition to co-registration, Wu et al. [21] presented an endeavor to integrate the Chang’E-1 imagery and laser altimeter data for consistent and precision lunar topographic modeling through a combined adjustment with a local surface constraint, the LAPs, image EO parameters, and tie points collected from the stereo images were the participants, and the output included the refined image EO parameters and laser ground points. The proposed combined adjustment approach can reduce the mis-registrations between the imagery and the laser altimeter data by a maximum of 1–18 pixels in image space. However, the method introduced the local surface constraint in which the stereo images were adjusted to optimize the fit to a surface interpolated from the laser altimeter data and it failed to consider the orientation errors of the laser altimeter data. This will bring about error when the local surface interpolated from the laser altimeter data cannot describe the real surface accurately. Moreover, the method works only under the condition that there are enough laser ground points.

To support Earth observation applications, in 2003, the United States launched an Earth observation laser altimetry satellite, called the Ice, Cloud and land Elevation Satellite (ICESat), [22], which is only equipped with GLAS and not with optical cameras. Most research uses ICESat altimeter data, as an elevation control in the block adjustment as Teo et al. [23] and Yamanokuchi et al. [24]. In this research, additional datasets (e.g., data from DEMs or laser altimeter data) were used as an absolute control and the adjustments were only made to the image data to adjust the images to fit the additional datasets. However, the accuracy of the final topographic models is dependent on the accuracy of the additional datasets. Thus, when there are only a few laser points, the method will not work well.

With the continuous improvement of China’s earth observation technology and the demand for large-scale stereo mapping, the idea of introducing geometric constraints provided by laser altimeter data to improve positioning accuracy has been put into practice. The recent Chinese ZY3-02 mission successfully produced 44 orbiter laser altimeter datasets. Subsequently, Tang et al. [25] and Xie et al. [26] built a rigorous geometric calibration model with pointing and ranging for correcting system errors of the laser altimeter onboard the ZY3-02 satellite, and this work realized high-precision in-orbit geometry calibration. It has been verified by Xie et al. [26] that the laser precision is about 2–3 m in areas with a slope less than 2°, and the absolute accuracy is better than 1 m in flat areas after calibration. Li et al. [27] performed experimental investigations of the integration of ZY3-02 satellite laser altimetry data and optical stereo images using combined adjustment by RFM with laser elevation constraint (RFM_EC) or RSM with laser ranging constraint (RSM_RC) with all result better than 3 m. However, there are two defects in their methods. On the one hand, the result obtained by RFM_EC is dependent on the laser altimeter data. On the other side, although the work also proposed the ranging constraint when using RSM, it failed to consider there exist attitude error for the altimeter because altimeter and stereo cameras are onboard the satellite.

This paper presents an investigation into the integration of stereo imagery and laser altimeter data, which is achieved through combined adjustment with one constraint in image space so as to improve Earth topographic modeling. Compared with previous work, the proposed method highlights the following two aspects: (1) Laser ranging information is used as to improve the elevation accuracy of stereo images; (2) a constraint to minimize the back-projection discrepancies in image space is incorporated into the combined adjustment. Stereo images and laser altimeter data covering research areas in Songshan, Tianjin, Dengfeng, and Taihang were collected to conduct experiments, and a total of seven kinds of adjustment schemes were designed.

2. Materials and Methods

2.1. Specifications and Geometric Modeling of ZY3-02 Stereo Images and Laser Altimeter Data

2.1.1. Specifications of ZY3-02

Being the first operational surveying and mapping satellite for the development of China’s civil space infrastructure (CCSI) and the second satellite in the ZiYuan3 series [28], the ZY3-02 satellite was launched successfully on May 30, 2016 [29]. Phased 183° from each other and operating as a true constellation on the same orbit, the ZY3-01 and ZY3-02 constellation has enabled the revisit cycle to be shortened from 5 days to 3 days, which has increasingly improved the country’s ability to anticipate risks, managing crises effectively, and speed up the collection of data over large areas.

ZY3-02 is the first satellite equipped with an experimental laser altimeter for Earth observations in China. This configuration is mainly used to test the function and performance of the laser altimeter; it can be also used to explore the feasibility of obtaining high-precision elevation control data on Earth surface and the potential of using laser ranging data for increasing the precision of stereo mapping without using GCPs. In addition to the laser altimeter, other core payloads onboard the ZY3-02 include triple linear array cameras (TLCs), a multi-spectral camera, a dual-frequency Global Positioning System (GPS) receiver, star sensor, etc. For the ZY3-02, the orbit accuracy reaches the centimeter level and the attitude stability is approximately 5 × 10⁻⁴ °/s. A schematic diagram showing the installation of the laser altimeter and the TLCs is illustrated in Figure 1, wherein the pointing of the laser altimeter is substantially parallel to the visual axis of the nadir camera, with both pointing to the center of the Earth.

As illustrated in Figure 2, the TLCs onboard ZY3-02 are comprised of nadir-view (NAD), forward-view (FWD), and backward-view (BWD) linear CCD cameras, where the viewing angles of the FWD and BWD cameras are installed at an inclination of ±22° from the NAD to allow for the acquisition of stereo pairs along the flight direction. The number of CCDs in the FWD and BWD cameras focal plane has been decreased from four to three and the size has been changed from 10 μm to 7 μm with respect to the ZY3-01 satellite to enable ground resolution increases from 3.5 m to 2.7 m. The ground sample distance (GSD) for NAD imagery is still 2.1 m. For each camera among the TLCs, there are three CCDs with 8192 pixels per CCD and approximately 30 pixels between the overlapping areas of adjacent CCDs.

The laser altimeter onboard ZY3-02 is a test load and its main design parameters are shown in

Table 1. Design parameters for the laser altimeter onboard ZY3-02.

Parameters	Designed Value
Laser pulse width	7 ns
Repetitive frequency	2 Hz
power consumption	45 W
Weight	40 kg
Single pulse energy	200 mJ
Effective distance	480~520 km
Laser footprint size	>75 m
center wavelength	1064 nm
Instantaneous emission field of view	0.1 mrad
Instantaneous receiving field of view	0.5 mrad
Detection probability	≥95% (flat area @ 500 km)
False alarm probability	≤5 × 10−4 (flat area @ 500 km)
Ranging accuracy	1 m (3σ, slope < 2°)

[26,30]. Due to the limited mounting space, its volume is small and the function is relatively simple. Besides ranging, there is no echo waveform recording system and no laser footprint camera for the laser altimeter. In comparisons with ICESat/GLAS [4,5], it was found that the single-pulse laser energy emitted by the laser altimeter onboard ZY3-02 is about 200 mJ, which is higher than the 75 mJ of GLAS. Additionally, the ZY3-02 satellite is lower than GLAS in terms of the repetition frequency, while the designed ground size of the laser footprint spot is equivalent to that of GLAS. When the slope is less than 2°, the atmospheric transmittance is good, and when the surface reflectivity is high, the accuracy of laser ranging can be better than 1 m based on the geometric calibration and validation for the laser altimeter onboard ZY3-02 [26,30]. Although there is still a certain gap in comparison to the 0.1 m ranging accuracy level of GLAS, it appears promising to use these data to achieve ZY3-02 1:50,000 mapping without ground control points.

2.1.2. Geometric Model of TLCs

The construction of the laser altimeter and optical geometric imaging models formed the basis for the geometric positioning method and establishment of combined adjustment equations. Therefore, this section introduces the laser geometric imaging model, the optical geometric imaging model, and its inverse transformation model, which were used to construct the basic model needed for this research.

A typical rigorous geometric model for an optical linear camera is collinear equation [31,32], and hence, the geometric model for the TLCs hence was established as follows:

$${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{WGS84}={\left[\begin{array}{c}{X}_{GPS}(t)\\ Y_{GPS}(t)\\ Z_{GPS}(t)\end{array}\right]}_{WGS84}+m\cdot R_{Body}^{WGS84}(t)\cdot {\left[\begin{array}{c}x-x_0\\ y-y_0\\ f\end{array}\right]}_I$$

(1)

where ${\left[\begin{array}{ccc}{X}_{GPS}(t)& Y_{GPS}(t)& Z_{GPS}(t)\end{array}\right]}_{WGS84}^T$ are the coordinates of satellite in the World Geodetic System 1984 (WGS84), which are determined by precise orbit determination of ZY3-02 satellite [33]. $R_{body}^{WGS84}$ is the rotation matrix from the satellite body fixed coordinated system to WGS84. ${\left[\begin{array}{ccc}{X}_{GPS}(t)& Y_{GPS}(t)& Z_{GPS}(t)\end{array}\right]}_{WGS84}^T$ and $R_{Body}^{WGS84}(t)$ together form the EO elements and both vary with scanning time t. ${\left[\begin{array}{ccc}x-x_0& y-y_0& f\end{array}\right]}_I^T$ represent the direction of the light ray of a pixel in the satellite body fixed coordinate system, where I stands for the NAD/BWD/FWD camera. $\left[\begin{array}{cc}{x}_0& y_0\end{array}\right]$ are the image coordinates of the principal point and f is the focal length; these form the interior orientation (IO) parameters. ${\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_{WGS84}^T$ are the coordinates in WGS84 corresponding to the pixel, and m is a scale factor and an unknown.

It can be seen from the above Equation (1) above that the factors affecting the positioning accuracy are mainly the IO and EO parameters. The errors of IO parameters belong to the static error, and its errors can be accurately compensated after geometric calibration. In contrast, the errors of EO parameters including the measurement errors of the attitude and trajectory, will vary with time. To eliminate the errors caused by EO parameters, the offset matrix $R_u$ is used, in which the deviation between the real light ray and the uncorrected light ray is corrected [34,35]. Based on Equation (1), a corrected geometric model then can be reconstructed for the TLC of ZY3-02 as expressed in Equation (2):

$${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{WGS84}={\left[\begin{array}{c}{X}_{GPS}(t)\\ Y_{GPS}(t)\\ Z_{GPS}(t)\end{array}\right]}_{WGS84}+m\cdot R_{body}^{WGS84}(t)\cdot R_{u\_camI}{\left[\begin{array}{c}x-x_0\\ y-y_0\\ f\end{array}\right]}_I$$

(2)

$R_{u\_camI}$ is the offset matrix for corresponding camera I, and $R_{body}^{WGS84}(t)\cdot R_{u\_camI}$ is defined as follows:

$$R_{body}^{WGS84}(t)\cdot R_{u\_camI}=\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]$$

(3)

The above equation describes the calculation from image space to object space. For linear array camera, the inverse calculation from object space to image space can also be described as shown in Equation (4):

$$\begin{array}{l}x=x_0+f\frac{a_1(X_{GPS}(t)-X)+b_1(Y_{GPS}(t)-Y)+c_1(Z_{GPS}(t)-Z)}{a_3(X_{GPS}(t)-X)+b_3(Y_{GPS}(t)-Y)+c_3(Z_{GPS}(t)-Z)}\\ y=y_0+f\frac{a_2(X_{GPS}(t)-X)+b_2(Y_{GPS}(t)-Y)+c_2(Z_{GPS}(t)-Z)}{a_3(X_{GPS}(t)-X)+b_3(Y_{GPS}(t)-Y)+c_3(Z_{GPS}(t)-Z)}\end{array}$$

(4)

For convenience in the following derivations, the Equation (4) can be abbreviated as follows:

$$\begin{array}{l}x=f_{camI\_x}(R_{u\_CamI},X,Y,Z,I{O}_{camI},E{O}_{camI})\\ y=f_{camI\_y}(R_{u\_CamI},X,Y,Z,I{O}_{camI},E{O}_{camI})\end{array}$$

(5)

2.1.3. Geometric Model for the Laser Altimeter

A typical geometric model for laser altimeter is the range equation, which can be described as follows:

$${\left[\begin{array}{c}{X}_{lidar}\\ Y_{lidar}\\ Z_{lidar}\end{array}\right]}_{WGS84}={\left[\begin{array}{c}{X}_{GPS}(t)\\ Y_{GPS}(t)\\ Z_{GPS}(t)\end{array}\right]}_{WGS84}+\rho \cdot R_{body}^{WGS84}(t)\cdot \left[\begin{array}{c}{l}_x\\ l_y\\ l_z\end{array}\right]$$

(6)

where ${\left[\begin{array}{ccc}{X}_{GPS}(t)& Y_{GPS}(t)& Z_{GPS}(t)\end{array}\right]}_{WGS84}^T$ are the coordinates of satellite in WGS84 and $R_{body}^{WGS84}$ is the rotation matrix from the satellite body fixed coordinate system to WGS84. For the laser altimeter, ${\left[\begin{array}{ccc}{X}_{GPS}(t)& Y_{GPS}(t)& Z_{GPS}(t)\end{array}\right]}_{WGS84}^T$ and $R_{body}^{WGS84}$ are EO parameters and both vary with time; these values can be obtained by interpolation at any time t. ${\left[\begin{array}{ccc}{l}_x& l_y& l_z\end{array}\right]}^T$ represent unit direction vector of the laser beam in the satellite body fixed coordinate system, where $l_x^2+l_y^2+l_z^2=1$. $\rho$ is the range measured by the laser altimeter. ${\left[\begin{array}{ccc}{X}_{lidar}& Y_{lidar}& Z_{lidar}\end{array}\right]}_{WGS84}^T$ are the coordinates of the laser ground point in WGS84.

After the geometric calibration of laser altimeter, the main factors affecting the positioning accuracy are the errors of EO parameters. Similarly, the offset matrix is applied to compensate the errors of EO parameters. Then, Equation (6) can be reconstructed as follows:

$${\left[\begin{array}{c}{X}_{lidar}\\ Y_{lidar}\\ Z_{lidar}\end{array}\right]}_{WGS84}={\left[\begin{array}{c}{X}_{GPS}(t)\\ Y_{GPS}(t)\\ Z_{GPS}(t)\end{array}\right]}_{WGS84}+\rho \cdot R_{body}^{WGS84}(t)\cdot R_{u\_lidar}\cdot \left[\begin{array}{c}{l}_x\\ l_y\\ l_z\end{array}\right]$$

(7)

Equation (7) can be furtherly abbreviated as follows:

$$\begin{array}{l}{X}_{lidar}=F_{lidarX}(R_{u\_lidar},\rho ,E{O}_{lidar})\\ Y_{lidar}=F_{lidarY}(R_{u\_lidar},\rho ,E{O}_{lidar})\\ Z_{lidar}=F_{lidarZ}(R_{u\_lidar},\rho ,E{O}_{lidar})\end{array}$$

(8)

where $R_u{\underset{\_}{}}_{lidar}$ strands for the offset matrix to compensate the errors for the errors of EO parameters for the laser altimeter. Equation (7) also can be modified as follows:

$$\rho =\left|R^{-1}{}_{u\_lidar}\cdot R_{body}^{WGS84}{(t)}^{-1}\left(\left[\begin{array}{c}{X}_{lidar}\\ Y_{lidar}\\ Z_{lidar}\end{array}\right]-\left[\begin{array}{c}{X}_{GPS}(t)\\ Y_{GPS}(t)\\ Z_{GPS}(t)\end{array}\right]\right)\right|$$

(9)

The Equation (9) is the constraint of laser ranging parameter. Similarly, for convenience in following derivations, the Equation (9) is abbreviated as:

$$\rho =\left|f_{lidar}(R_{u\_lidar},X_{lidar},Y_{lidar},Z_{lidar},E{O}_{lidar})\right|$$

(10)

2.2. Integration of Stereo Images and Laser Altimeter Data for ZY3-02 Satellite

2.2.1. Overview of the Approach

A combined adjustment model with one constraint in image space is proposed for integrating the stereo images and laser ranging data of the ZY3-02 satellite. The block adjustment is constructed based on the principle of the intersection of corresponding light rays in the overlap region between stereo images. On the basis of the block adjustment model, the laser ranging information is introduced to form a combined adjustment model.

Given that the ZY3-02 experiment dataset used in this paper was level-0 data, the rigorous geometric positioning model for stereo cameras and laser altimeter data was built based on the Section 2.1 equations with initial IO and EO parameters. When the combined adjustment is employed, the initial EO parameters of the stereo images and laser altimeter data, the image coordinates, and the ground coordinates of stereo images’ tie points represent the main adjustment observations. The initial EO parameters can be obtained directly from the ZY3-02 level-0 data. The image coordinate observations are divided into the following two types: the tie points of stereo images and the image point coordinates of laser ground points. The stereo image tie points are the conjugate points of stereo images obtained by image registration. By applying the geometric model established by corresponding IO parameters and initial EO parameters, the image coordinates of laser ground points onto stereo images can be derived through the back-projection. There is an image space constraint when calculating the adjustment parameters in the combined adjustment for confining the process of adjustment. With the adjustment parameters, which mainly include the improved EO parameters, i.e., the offset matrix $R_u$, and with the improved ground coordinates of laser points as final outputs, the improved geo-positioning accuracy of Earth topographic models can be achieved. The proposed approach can be summarized as follows:

(1)

Based on the initial orientation parameters, the geometric models of stereo images and laser altimeter data are built respectively;

(2)

Obtain the tie points of stereo images by image registration and their ground coordinates through space intersection;

(3)

Calculate the ground coordinates of LAPs and obtain the corresponding image coordinates in stereo images by back-projection; by using the laser points projected on the NAD camera image as references, obtain the LAPs on the BWD and FWD camera images by registration; the discrepancies between these two image coordinates are defined as the constraint in the combined adjustment;

(4)

Establish the observation equations for the combined adjustment model with the constraint in image space and solve the adjustment parameters;

(5)

Iterate steps (2)–(4) until the residuals of the adjustment parameters meet the threshold;

(6)

Generate the improved topographical model with the refined image orientation parameters and laser ground points.

2.2.2. Constraints in the Combined Adjustment Model

During the process of combing observation collected by stereo cameras and laser altimeter, because of the existence of various errors, e.g., attitude error, hardware errors, time synchronization errors, errors caused by changes in the space thermal environment, there might be discrepancies between the calculated and true values. In particular, the calculated laser ground points may differ from their real locations. Also, the calculated ground coordinates of tie points through space intersection may deviate from their true locations. Moreover, the back-projected laser ground points on stereo images may not be conjugate points and represent the different textural feature. These deviations reflect the inconsistency between stereo images and laser altimeter data by which we can establish the constraint equations, thus achieving the combined adjustment for ZY3-02 stereo images and laser altimeter data.

The key of combined adjustment is to take advantage of the high precision of laser ranging data, thereby confining the accuracy along the vertical direction in stereoscopic mapping. Figure 3 presents a schematic diagram showing the geometric relationship before and after the combined adjustment. In Figure 3, the light rays emitted by the FWD and BWD cameras are black dotted lines, and the laser light is a blue solid line. With respect to the laser altimeter, FWD and BWD cameras shoot the same ground point sequentially. It can be seen from Figure 3a that for the same ground point (black grid point in the figure), the elevation calculated by stereo images and the elevation measured by the laser altimeter are obviously inconsistent. Thanks to the high precision of the orbit and laser range, this combination provides a reliable measurement in the reference coordinate system along the vertical direction. Subsequently, the inconsistency between stereo images and laser altimeter data can be used to establish the range constraint of the combined adjustment, which is realized by minimizing discrepancies in image space, and the corrected intersection point is obtained. As shown in Figure 3b, the elevation calculated by stereo images’ space intersection is consistent to that of the LAP after combined adjustment, which improves the elevation accuracy of the stereo mapping and lays a foundation for high-precision stereo mapping.

The above range constraint is from the perspective in object space. Reversely, the range constraint can be analyzed in image space. Due to inaccurate EO parameters and inconsistencies between stereo images and laser altimeter data, the back-projected image coordinates on stereo images for the same laser ground point are not necessarily the true conjugate points. The issue is just to be solved in the combined adjustment.

In the case of a laser footprint camera onboard the satellite, the image point of laser footprint camera corresponding to the laser ground point can be directly extracted from the laser footprint image, and then, the accurate homologous points of the stereo camera for the laser ground point can be further obtained by matching the laser footprint image with stereo images using a high-precision matching method; however, this is not feasible for ZY3-02 because it is not equipped with a footprint camera. Nevertheless, according to the ZY3-02 satellite designer, the NAD camera and the laser altimeter can simultaneously work together. Moreover, the installation angles of the them are approximately identical with a deviation angle 0.92° [25,36]. Therefore, the image coordinates for the NAD camera derived by back-projecting the laser ground point can be seen as the intended position of the laser ground point to a large extent. As shown in Figure 4, the red image points on stereo images were obtained by back-projecting the blue laser ground point, and the green image points for the FWD and BWD cameras were obtained by high-precision matching with respect to the back-projected image point for the NAD camera. Theoretically, the discrepancies between the back-projected image points’ coordinates and matched image points’ coordinates should be zero for the FWD and BWD cameras, which is the basic idea of image space constraints that minimize the discrepancies in image space.

2.2.3. Combined Adjustment Model

The observation equation system for the combined adjustment of the ZY3-02 stereo images and laser altimeter data can be described in matrix form as follows:

$$\begin{array}{ll}{v}_1=A{x}_1+B{x}_2-l_1& P_1\\ v_2=A{x}_1+B{x}_3-l_2& P_2\\ v_3=C{x}_3+D{x}_4-l_3& P_3\\ v_4=E{x}_5-l_4& P_4\end{array}$$

(11)

There are four types of observation equations in the above equation. The first and second observation equations are related to the image tie points and the back-projected image coordinates of the laser ground point on the stereo images. The third is the laser range constraint observation equation, and the fourth is the additional constraint equation in the image plane. Here, x₁ is the vector of the stereo cameras’ offset matrixes’ corrections, x₂ is the vector of corrections to the ground coordinates of the tie points, x₃ is the vector of corrections to the laser ground point, x₄ is the vector of corrections to the laser altimeter’s offset matrix, and x₅ actually includes x₁, x₃, and x₄. A, B, C, D, and E are the partial derivations corresponding to the unknown parameters, and l is the constant term of the error equation. The weight P_i (i = 1, 2, 3, 4) in Equation (11) represents the contribution of each type of observation and can weaken the influence of the gross error in the adjustment system to a certain extent. Normally, weight parameter P_i is used with 10 times of a priori standard deviation ${\sigma }_i$ for each type of observation, and then, recalculate the weights P_i of each observations after each iterative adjustment.

The first kind of observation equation is actually the classic block adjustment model for tie points identified from stereo images. After linearization of Equation (5), which describes the relationship between the image point and ground point, the concrete expression for the first observation equation is established as follows:

$$\begin{array}{l}{v}_{x\_I}=\frac{\partial f_{camI\_x}}{\partial R_{u\_camI}}d{R}_{u\_camI}+\frac{\partial f_{camI\_x}}{\partial (X,Y,Z)}d(X,Y,Z)\\ +f_{camI\_x}(R_{u\_camI0},X_0,Y_0,Z_0,I{O}_{camI},E{O}_{camI})-x_I\\ v_{y\_I}=\frac{\partial f_{camI\_y}}{\partial R_{u\_camI}}d{R}_{u\_camI}+\frac{\partial f_{camI\_y}}{\partial (X,Y,Z)}d(X,Y,Z)\\ +f_{camI\_y}(R_{u\_camI0},X_0,Y_0,Z_0,I{O}_{camI},E{O}_{camI})-y_I\end{array}$$

(12)

where the first and second coefficients represent the partial derivative of Equation (5) with respect to the offset matrix $R_{u\_camI}$ and the ground point coordinates $(X,Y,Z)$ corresponding to A and B in Equation (11), respectively. $d{R}_{u\_camI}$ and $d(X,Y,Z)$ are corresponding to x₁ and x₂ in Equation (11) respectively. $f_{camI\_x}(R_{u\_camI0},X_0,Y_0,Z_0,I{O}_{camI},E{O}_{camI})-x_I$ and $f_{camI\_y}(R_{u\_camI0},X_0,Y_0,Z_0,I{O}_{camI},E{O}_{camI})-y_I$ are the constant terms corresponding to the term l₁. $(X_0,Y_0,Z_0)$ represent the initial ground point coordinates of the image tie point determined by space intersection based on the initial values of the offset matrix $R_{u\_camI0}$, initial stereo images’ IO and EO parameters.

The second kind of observation equation is very similar with the first, but here the image points are derived by the back-projection of laser ground points on stereo images and they are possibly not conjugate points because of existing various errors, which can be corrected by the process of combined adjustment. Similarly, the concrete expression of the second observation equation is established as follows:

$$\begin{array}{l}{v}_{x\_I}=\frac{\partial f_{camI\_x}}{\partial R_{u\_camI}}d{R}_{u\_camI}+\frac{\partial f_{camI\_x}}{\partial (X_{lidar},Y_{lidar},Z_{lidar})}d(X_{lidar},Y_{lidar},Z_{lidar})\\ +f_{camI\_x}(R_{camI0},X_{lidar0},Y_{lidar0},Z_{lidar0},I{O}_{camI},E{O}_{camI})-x_I\\ v_{y\_I}=\frac{\partial f_{camI\_y}}{\partial R_{u\_camI}}d{R}_{u\_camI}+\frac{\partial f_{camI\_y}}{\partial (X_{lidar},Y_{lidar},Z_{lidar})}d(X_{lidar},Y_{lidar},Z_{lidar})\\ +f_{camI\_y}(R_{camI0},X_{lidar0},Y_{lidar0},Z_{lidar0},I{O}_{camI},E{O}_{camI})-y_I\end{array}$$

(13)

where the definitions of the first and second coefficients in Equation (13) are same as those in Equation (12). $d{R}_{u\_camI}$ and $d(X_{lidar},Y_{lidar},Z_{lidar})$ corresponding to x₁ and x₃ in Equation (11), respectively. $f_{camI\_x}(R_{camI0},X_{lidar0},Y_{lidar0},Z_{lidar0},I{O}_{camI},E{O}_{camI})-x_I$ and $f_{camI\_y}(R_{camI0},X_{lidar0},Y_{lidar0},Z_{lidar0},I{O}_{camI},E{O}_{camI})-y_I$ are the constant terms corresponding to the term l₂. $(X_{lidar\_0},Y_{lidar\_0},Z_{lidar\_0})$ represent the initial coordinates of laser ground point which are obtained based on the geometric model for the laser altimeter in Equation (10) by using the range, initial offset matrix $R_{u\_lidar0}$ and EO parameters.

The third kind of observation equation is a modification of Equation (10) constrained by the laser range as follows:

$$F_{d\rho }=\rho -\left|f_{lidar}(R_{u\_lidar},X_{lidar},Y_{lidar},Z_{lidar},E{O}_{lidar})\right|=0$$

(14)

After linearization, the concrete expression of the second observation equation is established as follows:

$$\begin{array}{l}{V}_{F_{d\rho }}=\frac{\partial F}{\partial R_{u\_lidar}}d{R}_{u\_lidar}+\frac{\partial F}{\partial (X_{lidar},Y_{lidar},Z_{lidar})}d(X_{lidar},Y_{lidar},Z_{lidar})\\ +\rho -f_{lidar}(R_{u\_lidar0},X_{lidar0},Y_{lidar0},Z_{lidar0},E{O}_{lidar})\end{array}$$

(15)

where the first and second coefficients represent the partial derivative of Equation (14) with respect to offset matrix $R_{u\_lidar}$ and the laser ground point coordinates $(X_{lidar},Y_{lidar},Z_{lidar})$ corresponding to C and D in Equation (11) respectively. $d{R}_{u\_lidar}$ and $d(X_{lidar},Y_{lidar},Z_{lidar})$ correspond to x₃ and x₄ in Equation (11) respectively.

$\rho -f_{lidar}(R_{u\_lidar0},X_{lidar0},Y_{lidar0},Z_{lidar0},E{O}_{lidar})$ are the constant terms corresponding to the term l₃ and $R_{u\_lidar0}$ is the initial values of the offset matrix for the laser altimeter.

The fourth kind of observation equation is the image space constraints which can be represented as follows:

$$\begin{array}{l}{F}_{dx}=x_{backprojected\_FWD/BWD}-x_{match\_FWD/BWD}=0\\ F_{dy}=y_{backprojected\_FWD/BWD}-y_{match\_FWD/BWD}=0\end{array}$$

(16)

where $x_{backprojected\_FWD/BWD}$ and $y_{backprojected\_FWD/BWD}$ are the image points’ coordinates of the back-projected laser ground points on FWD and BWD cameras and are used as adjustment parameters. $x_{match\_FWD/BWD}$ and $y_{match\_FWD/BWD}$ are the matched image points’ coordinates with respect to the back-projected image point on NAD camera and are used as constants. Based on Equations (5) and (8), the concrete expression of the fourth kind of observation equation is established as follows:

$$\begin{array}{ll}{V}_{F_{dx}}& =\frac{\partial f_{FWD/BWD\_x}(R_{u\_FWD/BWD},X_{lidar},Y_{lidar},Z_{lidar})}{\partial R_{u\_FWD/BWD}}d{R}_{u\_FWD/BWD}\\ & +\frac{\partial f_{FWD/BWD\_x}(R_{u\_FWD/BWD},X_{lidar},Y_{lidar},Z_{lidar})}{\partial X_{lidar}}\cdot \frac{\partial R_{F_{lidarX}}}{\partial R_{u\_lidar}}d{R}_{u\_lidar}\\ & +\frac{\partial f_{FWD/BWD\_x}(R_{u\_FWD/BWD},X_{lidar},Y_{lidar},Z_{lidar})}{\partial Y_{lidar}}\cdot \frac{\partial R_{F_{lidarY}}}{\partial R_{u\_lidar}}d{R}_{u\_lidar}\\ & +\frac{\partial f_{FWD/BWD\_x}(R_{u\_FWD/BWD},X_{lidar},Y_{lidar},Z_{lidar})}{\partial Z_{lidar}}\cdot \frac{\partial R_{F_{lidarZ}}}{\partial R_{u\_lidar}}d{R}_{u\_lidar}\\ & +f_{FWD/BWD\_x}(R_{u\_FWD/BWD0},X_{lidar0},Y_{lidar0},Z_{lidar0})-x_{match\_FWD/BWD}\\ V_{F_{dy}}& =\frac{\partial f_{FWD/BWD\_y}(R_{u\_FWD/BWD},X_{lidar},Y_{lidar},Z_{lidar})}{\partial R_{u\_FWD/BWD}}d{R}_{u\_FWD/BWD}\\ & +\frac{\partial f_{FWD/BWD\_y}(R_{u\_FWD/BWD},X_{lidar},Y_{lidar},Z_{lidar})}{\partial X_{lidar}}\cdot \frac{\partial R_{F_{lidarX}}}{\partial R_{u\_lidar}}d{R}_{u\_lidar}\\ & +\frac{\partial f_{FWD/BWD\_y}(R_{u\_FWD/BWD},X_{lidar},Y_{lidar},Z_{lidar})}{\partial Y_{lidar}}\cdot \frac{\partial R_{F_{lidarY}}}{\partial R_{u\_lidar}}d{R}_{u\_lidar}\\ & +\frac{\partial f_{FWD/BWD\_y}(R_{u\_FWD/BWD},X_{lidar},Y_{lidar},Z_{lidar})}{\partial Z_{lidar}}\cdot \frac{\partial R_{F_{lidarZ}}}{\partial R_{u\_lidar}}d{R}_{u\_lidar}\\ & +f_{FWD/BWD\_y}(R_{u\_FWD/BWD0},X_{lidar0},Y_{lidar0},Z_{lidar0})-x_{match\_FWD/BWD}\end{array}$$

(17)

where the first four coefficients correspond to x₅ in Equation (11) which is the partial derivative of Equation (5) with respect to the offset matrix and of Equation (8) with respect to the laser ground point. It should be noted that the $R_{u\_lidar}$ is same with $R_{u\_NAD}$ in Equations (12) and (13).

2.3. Introduction to the Experimental Area and Data Sources

To sufficiently validate the proposed integration model, ZY3-02 stereo images and laser altimeter data covering the Songshan, Taihang, Dengfeng, and Tianjin areas in China were collected. The Songshan, Dengfeng, and Tianjing datasets consisted of only one standard scene and associated laser altimeter data, where each standard scene included NAD, BWD, and FWD stereo images. The Songshan and Dengfeng research areas are mountainous areas with complex topography, while the Tianjin research area is located in the plain area. As shown in Figure 5a, the Taihang dataset consisted of one strip of scenes (12 standard scenes) and one strip of laser altimeter data (orbit 944). This dataset crossed the Taihang Mountains and covered a total area of around 55 km × 550 km; the height range varied from 75 to 2,750 m. The topography of the Taihang dataset was varied and complex, encompassing mountains, basins, hills, and plains. General information about the experimental data is listed in

Table 2. General information for the experimental data.

Area	Stereo Images		Laser Altimeter Data
	Acquisition Date	Elevation Range (m)	Acquisition Date	Track ID	* Location
Songshan	26 July 2016	(123, 1491)	31 July 2016	Orbit 944	Center-right
Taihang	26 July 2016	(75, 2750)			Center-right
Dengfeng	3 October 2016	(−1, 778)			Left
Tianjin	30 August 2016	(−7, 0)	4 September 2016	Orbit 1476	Center-right

* Location represents the relative position overlaid on the stereo images.

, where the location stands for the location of the laser altimeter in relation to stereo images. There were two location groups to verify the influence of the location of the laser altimeter in relation to stereo images on the results of the combined adjustment. It should be noted that the stereo images and laser altimeter data were not acquired simultaneously. Ideally, the acquisition time should be the same so that the image coordinates obtained by back projecting the LAPs on NAD camera imagery can be more consistent with the true value. However, because of the acquisition time interval between the stereo images and laser altimeter data was relatively short and the pointing angle between them was consistent, the error induced by the non-simultaneous acquisition could be neglected.

With the final output based on the improved EO parameters and improved ground coordinates of laser points, improved geo-positioning accuracy of Earth topographic models can be achieved. Therefore, the experiments were mainly used to verify the more intuitive geo-positioning accuracy by using check points (CPs). The CPs for the Songshan, Taihang, and Dengfeng research areas were manually extracted from a 1:5000 digital orthophoto map/digital elevation model (DOM/DEM) of Songshan, 1:2000 DOM/DEM of Dengfeng, and 1:2000 DOM/DEM of Tianjin, for which specific information is given in

Table 3. Specific information about the reference data.

Area.	* GSD of DOM (m)	Planimetric Accuracy of the DOM RMS (m)	Height Accuracy of the DEM RMS (m)	Center Latitude and Longitude
Tianjin	0.5	1	1.5	38.00° N, 112.52° E
Songshan	0.5	1	1.5	34.65° N, 113.55° E
Dengfeng	0.2	0.4	0.7	34.45° N, 113.07° E

* GSD denotes the ground sample distance.

. For Taihang, there were 621 GCPs available to check the performance of the combined processing technique. The WGS84 coordinates of the GCPs were determined by static GPS processing and with accuracies exceeding 0.1 m; most of them were located in road intersections and manually picked its’ corresponding image point coordinates (x,y). Besides, four plane control points (PCPs), i.e., ground points with only known plane coordinates and unknown heights, were used. Figure 5b presents a detailed view of a GCP in the field.

It should be noted, that invalid laser altimeter data existed because of the influence of the atmosphere, vegetation, complex terrain, and other factors, so not all of the original laser altimeter data were used in the combined adjustment. For the laser altimeter of the ZY3-02 satellite, because it has no full waveform recording function, only ranging information and relatively simple ranging attributes such as echo detection anomalies, transmission anomalies, and other basic parameters, invalid LAPs could only be filtered out according to the basic ranging parameters. Then, ground coordinates of laser footprint points were calculated and invalid laser points within the existing terrain reference data were eliminated, such as those within the Shuttle Radar Topography Mission (SRTM) data. Invalid laser points were eliminated by following steps:

1.

Calculate the ground coordinates $(X_{lidar},Y_{lidar},Z_{lidar})$ in the geocentric Cartesian coordinate system for laser footprint based on the geometric model of laser altimeter, and convert the coordinates $(X_{lidar},Y_{lidar},Z_{lidar})$ in the geocentric Cartesian coordinate system to the coordinates $(La{t}_{lidar},Lo{n}_{lidar},H_{lidar})$ in the geographical coordinate system;

2.

Second, extract the elevation value $H_{SRTM}$ that can be directly interpolated in SRTM data according to the plane coordinates $(La{t}_{lidar},Lo{n}_{lidar})$ of the laser footprint. It is noted that the $H_{SRTM}$ should be converted to the same height Datum with the elevation $H_{lidar}$;

3.

Compare elevation of the laser footprint $H_{lidar}$ with the elevation value of SRTM $H_{lidar}$ to calculate the ${\delta }_h$, and eliminate the laser points satisfying the condition ${\delta }_h>3\sigma$, $\sigma$ is the absolute elevation accuracy of SRTM data.

Besides, some plane control points (PCPs) in Dengfeng area were manually extracted from YG-13 SAR image. For the PCP, its plane coordinates in object space are known while the height coordinate is unknown. The plane accuracy of PCPs was within 1.5 m [37] and Figure 6 illustrates a detail view of PCP.

3. Results

3.1. Inconsistency between ZY3-02 Stereo Camera Images and Laser Ranging Data

Using the initial orientation parameters of ZY3-02 stereo camera images and laser altimeter, the geometric model of them can be established. Then, the ZY3-02 LAP can be overlaid on the stereo images by back-projecting the LAPs directly. Figure 7 shows the ZY3-02 stereo images and corresponding laser altimeter data at the Taihang research area. Figure 7a shows the planar view of the laser altimeter data overlaid on the stereo images, where the left image shows the NAD, the center one shows BWD and the right shows the FWD image. The distance between adjacent laser footprints in NAD camera image is almost 1706 pixels, whose corresponding distance of ground laser spot is about 3.4 km. Due to the existence of invalid laser data which has been eliminated between the laser footprints, the distance between some of the individual laser footprint is several times of 1706 pixels. For example, the distance between laser points 1 and 2 on the NAD image is about 5118 pixels, with two invalid data in the middle.

For the same LAPs overlaid on the NAD, BWD, FWD images, the image points on three-view images are not the homologous points which represent the same textural feature as shown in Figure 7b. The range calculated by stereo images’ space intersection is not conform with the laser range. Using the laser points projected on the NAD camera image as references, the LAPs on the BWD and FWD camera images can be obtained by registration which named matched image points. After making the statistics of discrepancies between the back-projected image points’ coordinates and matched image points’ coordinates for BWD and FWD images respectively, it can be found that for BWD image, the Root Mean Square Error (RMSE) are 2.7 pixels across track and 2.0 pixels along track respectively, for FWD image, the RMSE are 3.1 pixels across track and 2.9 pixels along track respectively. This shows that inconsistency between stereo images and laser altimeter data exists. The main error sources are the orientation parameters. The proposed combined adjustment is to utilize the laser range with high accuracy for correcting the orientation parameters and elimination of these inconsistencies.

3.2. Validation of the Proposed Combined Adjustment Model

3.2.1. Validation Using the Standard Scene

In the Songshan, Dengfeng, and Tianjin research areas, there were evenly distributed GCPs that were manually extracted from the DOMs and DEMs. These GCPs served as CPs to validate the accuracy of the combined adjustment. Considering that the picking accuracy was limited by the blur effect [38], the final accuracy of CPs was approximately within 1 m. As mentioned above, not all of the laser altimeter data were used in the combined adjustment because of the existence of invalid data. Given that the ranging accuracy of the laser altimeter of ZY3-02 can be better than 1.0 m when the slope is less than 2° and the laser footprint size is about 75 m, all of the laser altimeter data used in the experiment were prescreened by using DEM data and most of the data used were located in the plain area. Five types of adjustment experiments were carried out, and the specific plans of the experiments were designed as follows:

1.

Plan 1: the free-net block adjustment without GCPs;

2.

Plan 2: block adjustment with four GCPs distributed at four corners of the research area;

3.

Plan 3: proposed combined adjustment using one LAP located in the middle of the research area;

4.

Plan 4: proposed combined adjustment using two LAPs located in the middle and start or end of the research area;

5.

Plan 5: proposed combined adjustment using all of the available LAPs in the research area.

In addition to the five types of adjustment experiment, another two type were designed for attempting to improve the plane accuracy in the combined adjustment.

6.

Plan 6: proposed combined adjustment using the additional PCPs in the research area;

7.

Plan 7: proposed combined adjustment using all of the available LAPs and the additional PCPs simultaneously in the research area.

Specifically, for the image tie points with known PCPs, its plane coordinates (X,Y) in object space did not participate in the adjustment and were deemed constants $(X_{PCP},Y_{PCP})$. Only the height value of the tie points with known PCPs is to be solved. This can be accomplished with the following modified version of Equation (12):

$$\begin{array}{l}{v}_{x\_I}=\frac{\partial f_{camI\_x}}{\partial R_{u\_camI}}d{R}_{u\_camI}+\frac{\partial f_{camI\_x}}{\partial (Z)}d(Z)\\ +f_{camI\_x}(R_{u\_camI0},X_{PCP},Y_{PCP},Z_0,I{O}_{camI},E{O}_{camI})-x_I\\ v_{y\_I}=\frac{\partial f_{camI\_y}}{\partial R_{u\_camI}}d{R}_{u\_camI}+\frac{\partial f_{camI\_y}}{\partial (Z)}d(Z)\\ +f_{camI\_y}(R_{u\_camI0},X_{PCP},Y_{PCP},Z_0,I{O}_{camI},E{O}_{camI})-y_I\end{array}$$

(18)

Figure 8 shows the points’ (CPs, GCPs, LAPs, PCPs) distribution in the Songshan, Dengfeng, and Tianjin area.

As shown in

Table 4. Accuracy statistics for the proposed combined adjustment under different types of points’ (CPs, GCPs, LAPs) distribution in the Songshan, Dengfeng, and Tianjin areas.

Area	Plan	No. of LAPs/PPCs	No. of GCPs/CPs	RMSE of CPs (m)			* CE90 (m)	* LE90 (m)
				East	North	Height
Songshan	P1	0/0	0/40	13.84	1.96	13.93	15.80	16.73
	P2	0/0	4/36	1.61	1.46	2.37	3.36	3.95
	P3	1/0	0/40	13.30	2.02	2.33	15.03	3.79
	P4	2/0	0/44	13.31	2.05	2.35	15.31	3.76
	P5	3/0	0/26	13.30	2.02	2.34	15.31	3.77
Tianjin	P1	0/0	0/40	26.35	7.84	18.51	28.90	20.51
	P2	0/0	4/36	1.42	1.24	1.86	2.91	1.64
	P3	1/0	0/4	25.38	7.99	1.59	28.00	2.49
	P4	2/0	0/40	25.34	8.00	2.03	27.97	1.81
	P5	7/0	0/40	25.43	7.98	1.72	28.05	2.29
Dengfeng	P1	0/0	0/45	16.91	12.01	20.99	21.89	23.01
	P2	0/0	4/41	1.57	0.98	1.58	`2.78	2.01
	P3	1/0	0/45	16.93	12.04	1.60	21.91	2.65
	P4	2/0	0/45	16.93	12.05	1.66	21.91	2.55
	P5	5/0	0/45	16.93	12.04	1.63	21.92	2.60
	P6	0/4	0/45	1.97	1.48	14.45	3.47	16.42
	P7	5/4	0/45	2.01	1.42	1.61	3.45	2.63

* CE90 and LE90 are the circular and linear errors at the 90th percentile respectively.

, the height accuracy in the research areas improved significantly after integration with the laser altimeter data by using the proposed combined adjustment. In Songshan, for example, the height accuracy improved from 16.73 m to 3.79 m by using one LAP, which was equivalent to that obtained from block adjustment with four GCPs, thus demonstrating the effectiveness of the proposed combined adjustment of ZY3-02 stereo images and laser altimeter data. Besides, the improved height accuracy was found to be insensitive to the number of LAPs and the location with respect to stereo images. On the one hand, when different numbers of laser points were used in the combined adjustment, the results improved almost to the same extent and even one LAP achieved good results. However, it is worthwhile to mention that all of the available LAPs in the research area were involved in the processing of the combined adjustment for the results with high reliability. On the other hand, the location of laser altimeter data on the stereo images was on the left side for Tianjin and in the center–right side for the other areas, but the final results revealed that the relative position had little effect on the improved efficiency.

Moreover, when both PCPs and LAPs were used, it can improve the plane and elevation accuracy at the same time, which is slightly lower than the results of block adjustment with four GCPs. This is because the solutions of combined adjustment are constrained by the known PCPs. After PCPs were allowed to participate in the proposed combined adjustment, the method was found to be promising for improving the plane and elevation accuracy simultaneously, and thus, additional improvements in topographical modeling could be achieved. However, it should be noted that the PCPs were directly used as constants rather than as adjustment parameters, which means that the results for the planimetric accuracy will be decided by the accuracy of the PCPs.

To further quantitatively assess the performance of the approach, two DEMs in Dengfeng were generated based on the P6 and P7 combined adjustment results respectively. The Semi-Global Matching (SGM) method [39] was introduced to acquire the point cloud, and then, the digital surface model (DSM) at a 5 m sampling distance was directly generated from the point cloud. The corresponding reference DEM in Dengfeng was used, and detailed information about Dengfeng is listed in Table 3. The generated DEM was tested for different terrains such as mountainous, hilly and flat plain areas to comprehensively evaluate the DEM quality. As show in

Table 5. Statistics for the digital surface model (DSM) elevation accuracy in Dengfeng area over various terrain types.

Area	Terrain	Plan	Minimum (m)	Maximum (m)	MEAN (m)	RMSE (m)
Dengfeng area1	Mountainous	P6	−24.29	61.70	22.94	5.58
		P7	−54.30	40.22	3.61	4.42
Dengfeng area2	Hilly	P6	−8.48	61.20	24.64	4.98
		P7	−23.23	30.12	1.02	4.53
Dengfeng area3	Plain	P6	−0.32	51.80	22.96	3.53
		P7	−23.55	26.44	1.76	2.98

, the results indicate that the participation of laser altimeter data can bring promising height accuracy.

As it can be seen in Figure 9a, in that the DEM is generated from the P6, there obviously exists a systematic error in the elevation difference resulted from where the terrain can even be observed. Figure 9b shows the DEM elevation difference using P7. The elevations are generally in line with the reference DEM. It can be noticed that there are places like black holes, which represent the larger elevation difference whatever in Figure 9a or in Figure 9b. The reference DEM used in the experiment were obtained in 2012 year, while the generated DEM were obtained in 2016 year. These mining areas were exploited after 2012 year. When making statistics of the difference, the difference in mining areas hence influence the result and cause large difference finally. Generally, the proposed combined adjustment can effectively improve the elevation accuracy.

3.2.2. Validation at the TaiHang Area

Unlike in Songshan, Dengfeng and Tianjin, there were total 12 standard scenes covering Taihang. All of the CPs in the stereo images were manually extracted based on the GCPs in the Taihang area with distinct features and in the plain area, and the accuracy of CPs was approximately within 1 m because of the accuracy achieved with the manual selection and the blur effect. Similarly, only some screened LAPs were used in the processing of combined adjustment. The detail points’ distribution were illustrated in Figure 10.

Table 6. Accuracy statistics for the proposed combined adjustment under different adjustment schemes in the Taihang area.

Area	Plan	No. of LAPs	No. of PCPs	No. of GCPs/CPs	RMSE of CPs (m)			CE90 (m)	LE90 (m)
					East	North	Height
Taihang	P1	0	0	0/245	12.28	2.21	8.29	14.69	10.51
	P2	0	0	4/241	1.83	2.00	2.48	4.13	2.15
	P3	1	0	0/245	12.02	2.55	2.79	14.43	2.79
	P4	2	0	0/245	12.01	2.55	2.57	14.41	3.04
	P5	13	0	0/245	11.98	2.55	2.01	14.39	2.87

shows the statistics for the positioning accuracy under different adjustment schemes in the Taihang area. The plane (CE90) and elevation accuracy (LE90) were 4.13 m and 2.15 m, respectively, when using four GCPs around four corners. When LAPs were used in the proposed combined adjustment, the elevation accuracy improved significantly, with the height accuracy being improved from 10.51 m to 2.87 m when using 13 LAPs. As can be seen in Figure 11, the previous inconsistencies were removed after the combined adjustment. The projected points on the stereo images were the homologous points, and the discrepancies were reduced.

4. Discussion

As mentioned previously, the image coordinates for the NAD camera derived by back-projecting the laser ground points can be treated as the intended position of laser ground points on the basis of the NAD camera and laser altimeter onboard ZY3-02 working together simultaneously and the installation angles of them being identical. However, the NAD camera images and corresponding laser altimeter data used in the experiments were captured at different times because simultaneously captured NAD images and laser altimeter data were unavailable due to satellite shooting schedule and limited data collection. Besides, due to the impact of the satellite launch and the physical environment changes in the orbit, the NAD camera and laser altimeter installation angles may slightly change. Therefore, the difference in the capture dates of the two datasets and the slightly different viewing angles could result in accuracy loss for the proposed combined adjustment, but the errors can be analyzed quantitatively. Figure 12 shows the elevation error analysis schematic.

S is the projected center of the NAD camera. The thick solid line represents the focal plane, and SO is the focal length f. The black laser altimeter represents an ideal situation, that it simultaneously works with the NAD camera, while the blue represents the real situation due to some unfavorable factors such as changes in the thermal environment and attitude stability, in which there exists a small deviation angle $△\theta$ between the true and real pointing angles that ultimately may change and cause variations in the relative geometric relation of the NAD camera. To simplify the analysis, suppose there is $△\theta$ only along the axis y. Here, only consider the factors outside the laser altimeter, so the difference of measured range between the ideal and true is zero. Laser point A projects to the NAD image at a, and laser point A’ projects to the NAD image at a’. The view angle $\angle aSO$ of image point a in the image coordinate is a. The relative location of LAP in relation to the image can result in different view angles $\alpha$. Thus, aa’ is the projection error resulting from the capture time interval.

To calculate aa’, we draw a line over A’ parallel to aa’, which intersects LA at B, and SA at C. Then, the projection error aa’ can be calculated by:

$$\begin{array}{ll}a{a}^{\prime }& =A^{\prime }C\cdot \frac{f}{S{O}^{\prime }}=A^{\prime }C\cdot \frac{f}{H-AB}\\ & =\left(H\cdot \sin △\theta +H\cdot \tan \alpha \cdot (1-\cos △\theta )\right)\cdot \frac{f}{H\cos △\theta }\end{array}$$

(19)

Taking the derivative of Equation (19) with respect to, the derivative of aa’ is:

$$\frac{\partial a{a}^{\prime }}{\partial \alpha }=H\cdot (1-\cos △\theta )\cdot {\sec }^2\alpha \cdot \frac{f}{H\cos △\theta }$$

(20)

Since the value of $\partial a{a}^{\prime }/\partial \alpha$ is always greater than 0, aa’ is a monotonically increasing function about variable $\alpha$. Thus, the maximum value of $\alpha$ at the end pixel of NAD camera image, aa’ is the maximum projection error.

Taking the derivative of Equation (19) with respect to $\theta$, the derivative of aa’ is:

$$\begin{array}{l}\frac{\partial a{a}^{\prime }}{\partial △\theta }=\left(H\cdot \cos △\theta +H\cdot \tan \alpha \cdot (1+\sin △\theta )\right)\cdot \frac{f}{H\cos △\theta }\\ +\left(H\cdot \sin △\theta +H\cdot \tan \alpha \cdot (1-\cos △\theta )\right)\cdot \frac{f}{H}\cdot \sec △\theta \cdot \tan △\theta \end{array}$$

(21)

Since the value of $\partial a{a}^{\prime }/\partial △\theta$ is also always greater than 0, $a{a}^{\prime }$ is a monotonically increasing function about variable $△\theta$.

For the ZY3-02 satellite, f is 1.7 m, the NAD image pixel size is 7 µm, and the minimum and maximum view angles are 0° and 2.94°, respectively. Because the interval time is 5 days, suppose the maximum deviation $△\theta$ is $5^{″}$ within such a short interval. According to the calculation results, the value of aa’ is approximately 5.88 pixels in the image (12.35 m in ground space), whatever minimum or maximum view angle is used, which demonstrates that the combined adjustment results are insensitive to the location of the LAP in relation to the stereo images. Moreover, all of the LAPs engaged in the experiment were selected at the plain area with slopes of less than 2°. The 12.35 m deviation in the planimetric space will bring about a maximum 0.43 m error in the vertical direction when the slope is 2°. Overall, the error is relatively small compared with the original height accuracy acquired from stereo images of the ZY3-02 satellite. To eliminate the errors caused by capture time interval, the NAD image and laser altimeter data that are captured simultaneously remain the best choice to use for experiments.

5. Conclusions

In this study, a combined adjustment model with one constrain in image space was proposed for integration of ZY3-02 stereo images and laser altimeter data for improved Earth topographic modeling. This proposed approach takes advantage of the accurate laser range to improve the positioning accuracy of stereo images, especially the elevation accuracy. Stereo images and laser altimeter data in several areas were collected as experiment data for validation work. Several conclusions can be drawn from this research, and these are summarized below:

(1)

Inconsistencies exist between the stereo images and laser altimeter data for the ZY3-02 satellite. For the same LAP overlaid on the NAD, BWD, and FWD images, the image points on three-view images are not the homologous points, and this is mainly caused by the errors of orientation parameters. The proposed combined adjustment utilizes the laser range with high accuracy for correcting the orientation parameters and eliminating these inconsistencies.

(2)

After applying the proposed combined adjustment, a superior elevation accuracy can be achieved. Experiments demonstrated that the height accuracy can be improved significantly to the extent of that obtained from block adjustment with four GCPs. Validation work further proved that the improved height accuracy is insensitive to the number of LAPs and locations of LAPs in relation to stereo images. It is suggested that all of the available LAPs in the research area should be involved in the processing of the combined adjustment to obtain results with high reliability.

(3)

By introducing the additional PCPs engaged in the combined adjustment, the planimetric accuracy can be improved along with the elevation accuracy. In the Dengfeng area experiment, the adjustment results derived by using LAPs and additional four LAPs were slightly lower than the results from the block adjustment with four GCPs.

Despite the promising results achieved for integration of stereo images and laser altimeter data of ZY3-02, so far only the elevation accuracy can be improved. This study attempts to improve the planimetric accuracy simultaneously by integration with PCPs, which can be obtained from the results of other sensors. Nevertheless, the final results will be decided by the accuracy of PCPs. Thus, further research about integration with other sensors that can acquire reliable observations in planar space is required, rather than directly using the results. That is, integration of data from multi-sensors is an important topic for future studies.