Combined Block Adjustment for Optical Satellite Stereo Imagery Assisted by Spaceborne SAR and Laser Altimetry Data

: To ensure the accuracy of large-scale optical stereo image bundle block adjustment, it is necessary to provide well-distributed ground control points (GCPs) with high accuracy. However, it is difﬁcult to acquire control points through ﬁeld measurements outside the country. Considering the high planimetric accuracy of spaceborne synthetic aperture radar (SAR) images and the high elevation accuracy of satellite-based laser altimetry data, this paper proposes an adjustment method that combines both as control sources, which can be independent from GCPs. Firstly, the SAR digital orthophoto map (DOM)-based planar control points (PCPs) acquisition is realized by multimodal matching, then the laser altimetry data are ﬁltered to obtain laser altimetry points (LAPs), and ﬁnally the optical stereo images’ combined adjustment is conducted. The experimental results of Ziyuan-3 (ZY-3) images prove that this method can achieve an accuracy of 7 m in plane and 3 m in elevation after adjustment without relying on GCPs, which lays the technical foundation for a global-scale satellite image process.


Introduction
Satellite image mapping plays an important role in many fields such as social economy and national defense. Remote sensing images acquired by satellites can achieve rapid update of geographic information data in large areas around the world, providing a rich source of geographic information data for decision-making and applications in related fields. In the application of satellite remote sensing images, the geometric positioning accuracy of images is the key indicator for measuring the availability of remote sensing satellite data, which has a great impact on the subsequent processing of remote sensing images and their application in large areas. To obtain high-precision satellite remote sensing images and improve the geometric positioning accuracy of images, a great deal of relevant research has been conducted in this field worldwide.
The geometric positioning accuracy of satellite images is mainly affected by the orbital attitude error on the satellite. To eliminate the influence of the system error and achieve high-accuracy geometric positioning of satellite images, the image positioning error can usually be eliminated by using a large number of ground control points (GCPs) through bundle block adjustment [1,2]. With the development of bundle block adjustment technology under sparse control conditions, the GCP-dependent adjustment method no longer requires control points for each image in the area, reducing the need for the number of control points [3]. However, this adjustment method does not completely eliminate the reliance on GCPs, and the quality and quantity of ground control points are still crucial to ensure the positioning accuracy of the images after adjustment in a large-scale area [4]. The time consumption and cost required to maintain and collect these GCPs are also important factors to be considered when using GCP-based block bundle adjustment. In addition, there will inevitably be a large number of areas with no or sparse control points such as areas with plateaus, canyons, and desert terrain, which poses a great challenge to realizing image processing on a global scale. Therefore, the processing of an adjustment method without GCPs is beginning to receive widespread attention because of its advantage of not relying on GCPs.
A satellite images process without GCPs usually refers to improving the relative consistency of the geometric positioning accuracy of satellite images without relying on GCPs measured by filed measurement, and to improving the absolute accuracy of ground positioning as much as possible. In recent years, with an increasing variety in available public geographic data and the improvement of positioning accuracy, the use of geographic data in block adjustment has become an effective method. These geographic information data include other aerospace/aerial images, digital orthophoto map (DOM), and digital elevation maps (DEMs), which have high positioning accuracy compared with the original images. At the same time, as the sources of publicly released geospatial data around the world are enriched, such as Shuttle Radar Topography Mission (SRTM) DEM, Advanced Spaceborne Thermal Emission (ASTER) DEM, Advanced Land Observation Satellite World 3D (AW3D) DEM and other data, Google Earth and other image data are gradually becoming optional ways of block adjustment without GCPs. By using these data as the control source, the image positioning accuracy can be greatly improved even without using GCPs.
For example, Wang et al. used a small number of ground control points and a variety of grid-size DEMs as control data to conduct a planar adjustment experiment with Ziyuan-3 (ZY-3) nadir-view (NAD) images in several regions, in which the accuracy of the checkpoints after adjustment was better than 7 m in an experiment with 139 images [5]. Zhang et al. investigated the process of determining ground point coordinates using DEM data as an elevation constraint, the priori weight of each observation was discussed, and a more robust adjustment strategy was proposed. The experimental results using ZY-3 images show that when the intersection angle between images is large, the accuracy of DEM and the density of tie points have obvious effects on the improvement of results [6]. In the process of block adjustment of Luojia-1 nighttime light remote sensing images, Li Xin et al. used 30 m grid SRTM DEM as auxiliary data to complete the adjustment of nighttime light images covering the whole of China, 25 control points were selected from the Google Earth map to participate in the adjustment, and the positioning accuracy of the adjusted images was better than 1.5 pixels [7]. Cao et al. proposed an adjustment method based on equivalent geometric sensor model (EGSM) that can make full use of the terrain information provided by DEM to improve the planimetric and elevation accuracy simultaneously. Based on a block adjustment test with 143 stereo images of ZY-3, the planimetric and elevation accuracies of the images after adjustment were improved from the original 17.3 m and 2.6 m to 2.5 m and 1.5 m respectively, and the results showed that the elevation accuracy achieved by AW3D as control data was basically comparable to the results using a small number of control points [8].
In addition, the laser altimetry data provided by laser rangefinders are also an effective source in adjustment without GCPs. Laser rangefinders can provide highly accurate range information by recording the interval between the transmitted and received laser beams with the laser measurement payload carried on board. At present, the major international on-board laser systems include the Near Earth Asteroid Rendezvous (NEAR) Laser Rangefinder (NLR) system [9], the measuring accuracy of which is better than 6 m; the Mars Orbiter Laser Altimeter (MOLA) [10], which is used to monitor small celestial bodies and has a measuring accuracy of approximately 37.5 cm; and the Chang'e-1 program [11], which is used by China for lunar exploration and the laser altimeter onboard has an accuracy of 1 m; and the Ice, Cloud, and Land Elevation Satellite (ICESat) with the Geoscience Laser Altimeter System (GLAS) system [12] with a measuring accuracy of 15 cm and its successor, ICESat-2 with the ATLAS system [13], which has a higher measuring accuracy better than 0.1 m.
The laser data provided by these satellite laser systems can provide detailed surface information and play a great value in remote sensing monitoring of the planets and the Earth. In addition to the fields of glacier monitoring [14], lake height variation monitoring [15], and vegetation inversion [16], laser altimetry data can be used for accuracy checking of DEM data [17,18] and assisting satellite images' precise positioning.
Wu et al. used the triple linear array cameras (TLCs) and the laser altimetry system carried by the Chang'e-1 satellite to perform the joint process of laser points and optical images by using local surface constraints to obtain a 3D digital model of the lunar surface, which provided a basis for the lunar surface landing site selection of the subsequent lunar exploration mission [19]. Li et al. used both a rational function model (RFM) and rigorous model for joint adjustment with laser points as elevation constraints to verify the feasibility of laser altimetry points (LAP) in block adjustment, and the results showed that both methods could achieve a positioning accuracy better than 3 m after adjustment [20]. Wang et al. proposed an effective laser altimetry point selection, registration, and adjustment method, and conducted experiments using mapping satellite-1 China Guangdong province data, and improved the elevation accuracy from 5.88 m to 2.51 m after using 1839 automatically extracted laser control points. It is concluded that elevation accuracy can be effectively improved by using LAPs, while plane accuracy remains unchanged [21]. Zhang et al. proposed an adjustment method based on a rigorous geometric model for using the laser system carried by the ZY3-02 satellite. The method combines the orientation parameters of the image and the ranging information of the laser data to establish a joint adjustment model, which can improve the elevation accuracy by eliminating the positioning consistency error between the image and the laser points [22].
In summary, the introduction of readily available public geographic data in the optical satellite stereo image block adjustment can improve the geometric positioning accuracy of images without GCPs, which is of great significance for realizing satellite mapping in large areas. Considering the high planimetric position accuracy of spaceborne SAR images and the high elevation accuracy of satellite-based laser altimetry data, this paper proposes a method combining the two types of data to assist the block adjustment of optical stereo images. The method first obtains the planar control points (PCP) by obtaining the tie points of optical images and SAR images through multimodal image registration, then obtains reliable LAPs by constructing a selection mechanism for a set of effective laser points, and finally introduces two kinds of control data in the block adjustment process, so as to achieve accuracy improvement in both plane and elevation directions. The method proposed in this paper is validated through experimental implementation in several regions and discussed based on the experimental results.

Materials and Methods
This paper presents a joint optical satellite stereo images block adjustment method using the DOM generated from a SAR image and the LAPs as the control source. The method first uses the registration method based on the similarity measure to achieve the tie point extraction between the optical stereo image and SAR DOM, then uses the LAPs selection method to select the LAPs and then project the LAPs into the stereo image. Finally, the observation equation of image point coordinates is established separately, the constraint of plane direction and elevation direction is realized by combining SAR DOMbased PCPs and LAPs, respectively, and finally the image orientation correction parameters are obtained by solving the equation to improve the plane and elevation geometric accuracy of the stereo image. The flow chart is shown in Figure 1.  Figure 1. Flow chart of method.

Automatic Acquisition of PCPs
Currently, the commonly used similarity metrics such as normalized cross correlation and mutual information mainly use the intensity information between images to match homonymous points, and their registration performance decreases sharply when there are significant nonlinear radiometric differences between images. However, compared with grayscale information, the geometric structural features of images are less affected by nonlinear radiometric differences, so using geometric structural features to construct similarity measures can break through the deficiency that traditional similarity measures are mainly matched by grayscale information and are more sensitive to nonlinear radiometric differences, thus enhancing the robustness of image registration and providing the basic theory and technical support. Figure 2 shows a pair of optical and SAR remote sensing images, and it can be found that although the nonlinear radiometric differences between the images are large, the geometric structural features between the images have high similarity. Considering this feature, with the idea of histogram of oriented gradient used in the field of target recognition to represent the shape and structural features of images, the magnitude and orientation of structure representation for consistency judgment are established, the feature descriptor named the histogram of orientated phase congruency (HOPC) for geometric structure similarity between images is realized on this basis, and the similarity measure for image registration is established on this basis. The phase-congruency-based registration method mainly uses the spatial structure of the same feature between heterogeneous images to construct a HOPC descriptor, which can reflect the similarity of spatial structure. The established descriptors can calculate the similarity between heterogeneous images and finally achieve accurate registration between heterogeneous images of different sources.

Automatic Acquisition of PCPs
Currently, the commonly used similarity metrics such as normalized cross correlation and mutual information mainly use the intensity information between images to match homonymous points, and their registration performance decreases sharply when there are significant nonlinear radiometric differences between images. However, compared with grayscale information, the geometric structural features of images are less affected by nonlinear radiometric differences, so using geometric structural features to construct similarity measures can break through the deficiency that traditional similarity measures are mainly matched by grayscale information and are more sensitive to nonlinear radiometric differences, thus enhancing the robustness of image registration and providing the basic theory and technical support. Figure 2 shows a pair of optical and SAR remote sensing images, and it can be found that although the nonlinear radiometric differences between the images are large, the geometric structural features between the images have high similarity. Considering this feature, with the idea of histogram of oriented gradient used in the field of target recognition to represent the shape and structural features of images, the magnitude and orientation of structure representation for consistency judgment are established, the feature descriptor named the histogram of orientated phase congruency (HOPC) for geometric structure similarity between images is realized on this basis, and the similarity measure for image registration is established on this basis.  Figure 1. Flow chart of method.

Automatic Acquisition of PCPs
Currently, the commonly used similarity metrics such as normalized cross correl tion and mutual information mainly use the intensity information between images t match homonymous points, and their registration performance decreases sharply whe there are significant nonlinear radiometric differences between images. However, com pared with grayscale information, the geometric structural features of images are les affected by nonlinear radiometric differences, so using geometric structural features t construct similarity measures can break through the deficiency that traditional similarit measures are mainly matched by grayscale information and are more sensitive to non linear radiometric differences, thus enhancing the robustness of image registration an providing the basic theory and technical support. Figure 2 shows a pair of optical and SAR remote sensing images, and it can be foun that although the nonlinear radiometric differences between the images are large, th geometric structural features between the images have high similarity. Considering th feature, with the idea of histogram of oriented gradient used in the field of target recog nition to represent the shape and structural features of images, the magnitude and or entation of structure representation for consistency judgment are established, the featur descriptor named the histogram of orientated phase congruency (HOPC) for geometr structure similarity between images is realized on this basis, and the similarity measur for image registration is established on this basis. The phase-congruency-based registration method mainly uses the spatial structur of the same feature between heterogeneous images to construct a HOPC descripto which can reflect the similarity of spatial structure. The established descriptors can ca culate the similarity between heterogeneous images and finally achieve accurate regi tration between heterogeneous images of different sources. The phase-congruency-based registration method mainly uses the spatial structure of the same feature between heterogeneous images to construct a HOPC descriptor, which can reflect the similarity of spatial structure. The established descriptors can calculate the Remote Sens. 2021, 13, 3062 5 of 19 similarity between heterogeneous images and finally achieve accurate registration between heterogeneous images of different sources.
The specific implementation usually adopts a coarse-to-fine registration scheme; firstly, a coarsely geometric correction is performed by using orbital attitude parameters or rational polynomial parameters from satellite images, which can eliminate the rotation and scale differences between images. Then, the registration metric is built based on geometric structure similarity to achieve tie point matching and mismatches are eliminated, and the obtained homonyms are used to correct the initial orbital attitude parameters or rational polynomial coefficient (RPC) parameters to finally realize the fine registration of images.
The process of constructing similarity descriptors for spatial geometric structures is shown in Figure 3. Firstly, a template window with a certain size is selected and then the phase congruency amplitude and orientation for each pixel in the template window are computed. Then, the template window is divided into overlapping blocks, which are called "cells", and each cell contains multiple units. Finally, a local histogram of phase congruency orientations is accumulated over all the pixels within the cells of each block, and collects all the blocks in the template window into a combined HOPC feature vector, a process similar to the SIFT [23]. The specific implementation usually adopts a coarse-to-fine registration scheme; firstly, a coarsely geometric correction is performed by using orbital attitude parameters or rational polynomial parameters from satellite images, which can eliminate the rotation and scale differences between images. Then, the registration metric is built based on geometric structure similarity to achieve tie point matching and mismatches are eliminated, and the obtained homonyms are used to correct the initial orbital attitude parameters or rational polynomial coefficient (RPC) parameters to finally realize the fine registration of images.
The process of constructing similarity descriptors for spatial geometric structures is shown in Figure 3. Firstly, a template window with a certain size is selected and then the phase congruency amplitude and orientation for each pixel in the template window are computed. Then, the template window is divided into overlapping blocks, which are called "cells", and each cell contains multiple units. Finally, a local histogram of phase congruency orientations is accumulated over all the pixels within the cells of each block, and collects all the blocks in the template window into a combined HOPC feature vector, a process similar to the SIFT [23].

Phase congruency magnitude and oriention
Divide the window into blocks consist of some cells

Accumulate orientation histograms for cells and blocks
Collect HOPCs for all blocks over template window cell block Figure 3. Construction process of histogram of orientated phase congruency (HOPC).
The HOPC descriptor captures the structure properties of images; the correlation coefficient between HOPC descriptors is used here as the similarity measure for image matching, and the corresponding fine registration scheme is designed based on this basis, which proceeds as follows: 1. using classical feature point detection operators (e.g., Harris, Forstner, etc.) to extract a large number of evenly distributed feature points as key points with a block strategy in the images. 2. using HOPC as a similarity measure and a template matching strategy to obtain homonymous points between images. 3. acquiring sub-pixel accuracy by a local extreme fitting technique and eliminating the mismatch points. 4. using the acquired homonymous points to correct the orbital attitude parameters or RPC parameters of the images to obtain the exact geometric transformation relationship between image coordinates and geographic coordinates, and to fine-correct the images to finally achieve exact registration between images.

Automatic Acquisition of LAPs
According to public information, the footprint diameter of a GLAS laser beam on the ground is about 70 m, and its plane accuracy is better than 15 m. If the positioning accuracy of the satellite image and the plane positioning accuracy of the laser point are the same, the relative error between the laser point and the image can be considered small, and the laser point can be introduced into the block adjustment process as an elevation constraint to improve the elevation accuracy after the block adjustment. GLAS laser altimetry data only records the coordinates of the laser point in latitude and longitude, the elevation value, and some auxiliary data. Although the nominal accuracy of GLAS data is The HOPC descriptor captures the structure properties of images; the correlation coefficient between HOPC descriptors is used here as the similarity measure for image matching, and the corresponding fine registration scheme is designed based on this basis, which proceeds as follows: 1.
using classical feature point detection operators (e.g., Harris, Forstner, etc.) to extract a large number of evenly distributed feature points as key points with a block strategy in the images. 2.
using HOPC as a similarity measure and a template matching strategy to obtain homonymous points between images. 3.
acquiring sub-pixel accuracy by a local extreme fitting technique and eliminating the mismatch points. 4.
using the acquired homonymous points to correct the orbital attitude parameters or RPC parameters of the images to obtain the exact geometric transformation relationship between image coordinates and geographic coordinates, and to fine-correct the images to finally achieve exact registration between images.

Automatic Acquisition of LAPs
According to public information, the footprint diameter of a GLAS laser beam on the ground is about 70 m, and its plane accuracy is better than 15 m. If the positioning accuracy of the satellite image and the plane positioning accuracy of the laser point are the same, the relative error between the laser point and the image can be considered small, and the laser point can be introduced into the block adjustment process as an elevation constraint to improve the elevation accuracy after the block adjustment. GLAS laser altimetry data only records the coordinates of the laser point in latitude and longitude, the elevation value, and Remote Sens. 2021, 13, 3062 6 of 19 some auxiliary data. Although the nominal accuracy of GLAS data is 15 cm, in the actual measurement process the laser will be scattered and refracted by the atmosphere after emission, and the topographic changes on the surface and forward scattering of clouds will also cause oversaturation of the laser echo signal, which will lead to the degradation of the laser ranging accuracy. Therefore, before using the laser altimetry data, it is necessary to selection the laser points first. The selection of laser points is mainly divided into two ways: one is to combine the DEM data for laser echo signal simulation, because there is a certain correlation between the echo signal and the terrain features, so by selecting the laser points corresponding to the appropriate echo signal, the accuracy of the selected laser points can be guaranteed [20]. The other is to directly use the DEM for laser point selection; by setting certain selection conditions, the laser points that do not meet the conditions are eliminated, and then the remaining laser points can ensure their positioning accuracy [21].
In this paper, we choose to use the DEM-based selection method, and considering that the diameter of laser spot is around 70 m, the data of 30 m-grid SRTM can be used to effectively extract the elevation variation within the laser spot range, and it can also be combined with the quality evaluation parameters recorded by laser data to assist the selection process. The laser point selection process based on SRTM 30 m DEM can be divided into the following steps, and Figure 4 shows a schematic diagram of the process.
(1) In the process of laser measurement, if there are thick clouds over the surface, the laser points may be reflected directly on the clouds, and at this time the laser signal recorded by the laser rangefinder is actually the elevation value of the clouds, so it is first necessary to use the DEM data to remove the laser points with large differences between the elevation value and the actual surface elevation value, and this process can be carried out by setting the threshold value of elevation difference. (2) If the slope of the terrain within the laser spot changes greatly, the increase of the slope will lead to the pulse broadening of the received laser signal. Since the laser spot diameter is 70 m and the positioning accuracy is about 15 m, if the slope within the spot range causes the elevation change of the ground surface to exceed the high accuracy of the laser spot, the error of plane positioning will cause the actual accuracy of the laser spot to decrease and affect the use of the laser spot. In order to reduce the influence of surface deformation on the laser point, the DEM can be used to estimate the change of slope within the laser spot range, and when the change of slope causes the laser height measurement error to exceed the laser point's high accuracy threshold, the laser point that exceeds the threshold will be rejected to ensure that the selected laser point's high accuracy meets the use requirements. (3) After selection of the two types of elevation anomalies of laser points using DEM, the remaining laser points can be further selected by using the quality evaluation parameters of the laser points themselves; mainly, the laser points whose atmospheric scattering gain and saturation correction parameters do not meet the requirements can be eliminated, and at this time the remaining laser points can ensure their measurement accuracy and can be extracted for LAPs.
After finishing the selection of each laser point, in order to use the laser control point in the adjustment, it is necessary to extract the coordinates of the corresponding image point of the laser control point on the optical image. Since the laser point spot diameter is about 70 m and its plane positioning accuracy is 10 ± 4.5 m, when the plane positioning accuracy of the satellite image is high (less than 15 m), it can be considered that the collocation error between the laser point and the plane position of the image in the region is small at this time, and the ground coordinates of the laser point can be used to obtain the image coordinates of the laser point on the reference image by projection directly using the RFM [24]. that the diameter of laser spot is around 70 m, the data of 30 m-grid SRTM effectively extract the elevation variation within the laser spot range, and combined with the quality evaluation parameters recorded by laser data lection process. The laser point selection process based on SRTM 30 m D vided into the following steps, and Figure 4 shows a schematic diagram o  (1) In the process of laser measurement, if there are thick clouds over t laser points may be reflected directly on the clouds, and at this time recorded by the laser rangefinder is actually the elevation value of the first necessary to use the DEM data to remove the laser points with la between the elevation value and the actual surface elevation value, a can be carried out by setting the threshold value of elevation differen (2) If the slope of the terrain within the laser spot changes greatly, the slope will lead to the pulse broadening of the received laser signal. spot diameter is 70 m and the positioning accuracy is about 15 m, if t the spot range causes the elevation change of the ground surface to e accuracy of the laser spot, the error of plane positioning will cause t racy of the laser spot to decrease and affect the use of the laser spot. duce the influence of surface deformation on the laser point, the DEM estimate the change of slope within the laser spot range, and when slope causes the laser height measurement error to exceed the las accuracy threshold, the laser point that exceeds the threshold will ensure that the selected laser point's high accuracy meets the use requ The extraction process of LAPs is mainly divided into three parts, Figure 5 shows the main process diagram.
(1) Firstly, the geographic coordinates of the LAPs can be obtained by using the latitude, longitude and elevation information of the LAPs itself, and after the conversion of the coordinate system, calculate the ground coordinates of LAPs and obtain the corresponding image coordinates in NAD images by back-projection. (2) Then, by using the laser points projected on the NAD camera image as references, obtain the LAPs on the backward-view (BWD) and forward-view (FWD) camera images by registration. (3) After obtaining the LAP's image coordinates on the NAD, BWD, and FWD images, in order to eliminate the mis-registration points, the LAPs can be regarded as the tie points to perform free network adjustment among the stereo images, and the tie points with errors larger than the threshold are eliminated until the error of all the points meets the requirements. At this time, the image coordinates and the corresponding object square coordinates of LAPs are obtained.

Combined Adjustment Model
At present, mainstream optical satellites are equipped with line array pushbroom sensors, and in order to avoid the problems of difficult to obtain the initial values of the parameters of the block adjustment solution based on the rigorous geometric model, strong parameter correlation, and unstable solution, considering the characteristics of near-parallel projection caused by the long focal length and narrow field-of-view angle of satellite photography, the RFM is generally used at present [25]. There are many error compensation models for RPC parameter systems b the image object, the most commonly used of which is the affine transformatio with six affine change parameters, which can absorb translation, rotation, and rors in the error well [5], in the following form. In order to introduce LAPs and PCPs in the adjustment process, the obs equations of the two control points need to be filled separately, and at this time t terms are calculated according to the following two categories. The RFM model establishes the transformation relationship between the ground coordinate (lat, lon, h) and the corresponding image coordinate (x, y) in the general form of [26].
   R n = P 1 (X n ,Y n ,Z n ) P 2 (X n ,Y n ,Z n ) C n = P 3 (X n ,Y n ,Z n ) Among them, (X, Y, Z) is the ground point object coordinates after regularization, (R, C) is the corresponding image point coordinates on the image after regularization, P 1 (X, Y, Z), P 2 (X, Y, Z), P 3 (X, Y, Z), P 4 (X, Y, Z) is a cubic polynomial function containing the ground point coordinates (X, Y, Z), the function X, Y, Z three coordinate components of the respective powers and their powers of the sum of the maximum does not exceed 3.
There are many error compensation models for RPC parameter systems based on the image object, the most commonly used of which is the affine transformation model with six affine change parameters, which can absorb translation, rotation, and scale errors in the error well [5], in the following form.
In Equation (3) The affine transformation model is introduced into the RFM, and the affine transformation parameters and the ground point coordinate correction parameters are solved to obtain the form of the matrix of the block adjustment [27,28].
where , l = ∆x ∆y T , P is the weight matrix.
In order to introduce LAPs and PCPs in the adjustment process, the observation equations of the two control points need to be filled separately, and at this time the filling terms are calculated according to the following two categories.
(1) Observation equation based on LAPs: Since direct free network adjustment does not significantly improve the positioning accuracy of all images after adjustment, the positioning accuracy of images after free network adjustment can be improved by introducing the ranging information of free network adjustment into the adjustment process and establishing an adjustment scheme with elevation constraints. The main difference lies in the fact that the plane coordinates of the object coordinates of the laser points are calculated by stereo intersection with the stereo image, while the elevation values of the intersection coordinates are calculated by using the elevation values of the LAPs. In this paper, the object coordinates of the laser point itself are corrected according to the adjustment results during each iteration, at which time the original error equation (Equation (3) , and its observation equation takes the following form: In the above equation, since only the elevation direction of the LAPs is utilized as the control point and its plane direction is the geographic coordinate acquired by the forward intersection, the error in its plane direction also needs to be added to the error equation.
However, in order to ensure that the elevation accuracy of the LAPs will not exceed the accuracy requirement after updating the correction parameters of the affine transformation, the elevation value h corr calculated after each iteration of the LAPs should comply with the following formula: In the above equation, h init is the initial elevation value obtained from laser altimetry data and σ Lidar is the known accuracy of the laser altimetry data used.
The adjustment iteration process is performed to correct the affine transformation parameters and update the ground point object coordinates, and then the iteration is repeated until it converges. In this process, the elevation coordinates are replaced by the elevation values of the LAPs themselves, and the plane coordinates are the plane coordinates obtained from the forward intersection. The elevation accuracy of the block adjustment based on LAPs is constrained by the elevation value of the laser points, so the elevation positioning accuracy of each image after adjustment depends mainly on the elevation accuracy of the laser control points used.
(2) Observation equation based on PCPs: Since the imaging method of satellite optical images is mainly line array pushbroom imaging, their positioning accuracy mainly depends on the measurement accuracy of on-board attitude orbit sensors. Compared with the optical image, the satellite SAR image is based on the Range-Doppler model, and the error affecting its positioning accuracy is mainly manifested as translation error on the ground, so the system error can be compensated by the on-orbit calibration method using the ground angular reflector, and the compensated image can achieve higher planar positioning accuracy [29][30][31]. Therefore, it can be considered to use the calibrated SAR image as the plane control source, and the plane direction can be constrained in the optical image block adjustment by extracting the control points between the SAR image and the optical image, so as to improve the plane positioning accuracy. Meanwhile, the object coordinates of the PCPs themselves are not modified in the adjustment process.
By using the PCPs extracted from SAR images as plane constraints, the adjustment scheme of PCPs is established. In the process of forward intersection to obtain ground point coordinates, their elevation coordinates are the elevation values obtained from direct image intersection, while their plane coordinates are replaced by the plane coordinates of the PCPs themselves. For the introduced PCPs, a suitable error equation needs to be established. Since its plane coordinates are regarded as the true value in the adjustment without correction, the difference is mainly for the partial derivative of the elevation term In addition, in order to ensure that the plane accuracy of the PCPs will not exceed the accuracy requirement after updating the affine transformation correction parameters during the adjustment solution iteration, the plane coordinates (x, y) corr should meet the following requirements in each iteration: In the above equation, (x, y) init is the planar coordinate of the PCPs itself, and σ SAR is the known accuracy of the SAR DOM data used. The method requires that the plane coordinates of the PCPs in each adjustment iteration do not exceed the accuracy range to ensure the accuracy of the adjustment results.
After listing the observation equations of the tie points, PCPs, and LAPs, the adjustment model in matrix form is obtained. In order to solve the compensation parameters of the image, the equation needs to be solved, and Equation (3) is usually transformed into the normal equation according to the least squares principle: where the P refers to weight matrix. Rewrite Equation (8) in the following form: When solving, when the number of ground points is large, the number of unknowns t will be much larger than s. Therefore, in order to reduce the number of unknowns to be solved, generally only the correction parameters of the affine coefficient of change of the image are solved, and the ground point unknowns can be solved in a modified way to eliminate the term, which can greatly improve the speed of solving, at which time we can obtain: The ground point coordinates can be updated by correcting the affine transformation coefficients after solving s by Equation (10), and then using the corrected parameters ahead of the rendezvous. When solving Equation (10), a variety of equation solving methods including the conjugate gradient method can be used for fast solution. In order to ensure the reliability of the adjustment, the adjustment results are usually optimized by multiple iterations. The results are updated by repeated iterations in the adjustment solution process until the change in two iterations is less than the set threshold, then the iteration is terminated and the adjustment results are output.

Experimental Data and Test Area
In order to verify the effectiveness of the proposed adjustment method, a total of 255 scenes of ZY-3 images of Jiangxi region, China, containing 86 stereo models, were collected from 24 December 2012 to 11 January 2014. In order to verify the positioning accuracy of the images in this region, a total of 134 ground points were collected by using Real-Time Kinematic Global Positioning System (RTK-GPS) measurements, which can ensure an accuracy better than 0.1 m and can therefore be used as independent check points (ICPs).
A total of 297 scenes of mapping satellite images of Hubei region, China within the time range of 27 April 2012 to 29 February 2018 were also collected, divided into 99 stereo models covering the whole province of Hubei, while a total of 134 ICPs were collected within the region for accuracy check, with a positioning accuracy better than 0.1 m. Figure 6 shows the image distribution in the two areas, and the detailed information of the two test areas is shown in Table 1. The topography in the Jiangxi test area is mainly plain, while hills and mountains exist in some areas with large elevation changes. The eastern part of the test area in Hubei Province is mainly plain terrain with a small amount of hills, while the western part is mainly covered with mountainous terrain, resulting in a large variation of topography in this test area.  6 shows the image distribution in the two areas, and the detailed information of th two test areas is shown in Table 1. The topography in the Jiangxi test area is mainly plain while hills and mountains exist in some areas with large elevation changes. The easter part of the test area in Hubei Province is mainly plain terrain with a small amount of hill while the western part is mainly covered with mountainous terrain, resulting in a larg variation of topography in this test area.    The SAR DOM used is a regional orthophoto of China was generated from Gaofen-3 (GF-3) image with FS2 mode, and the approximate range of this DOM is shown in Figure 7. Its ground sampling interval is about 5 m, and its positioning accuracy can be guaranteed within 7 m after verification, which can be used as control data to collect PCPs through HOPC matching and enhance the planar positioning accuracy of the image after joint adjustment. The laser data used is the L14 level product of GLAS, which mainly records more than 80 parameters such as laser footprint 3D coordinates, laser echo saturation correction, and gain value. Considering the diameter of laser spot is around 70 m, the data of 30 m-grid SRTM is used for elevation extract in the DEM-based LAP election method.

Combined Block Adjustment Test for ZY-3 Images
In this section, firstly, the registration test between ZY-3 and SAR DOM was carried out using the HOPC method. From the results, it can be seen that although the texture and radiation characteristics of SAR and optical images are different, the distribution of features actually has certain similarity in spatial scale, so that high accuracy tie points can be obtained in the case of features similar to houses as well as farmland. As shown in Figure 8, the registration between farmland and buildings can obtain good results, from which the pixel correspondence between optical images and SAR DOM can be effectively obtained, so that the latitude and longitude coordinates of the corresponding points can be read from the DOM, and finally used as PCPs to improve the planimetric positioning accuracy of the optical images after adjustment. At the same time, after selection and projection of the laser data, 9074 LAPs were obtained in the area of Jiangxi, and their distribution is shown in Figure 9a. The LAPs obtained after selection covers most of the image area and can provide enough LAPs as elevation control reference.

Combined Block Adjustment Test for ZY-3 Images
In this section, firstly, the registration test between ZY-3 and SAR DOM was carried out using the HOPC method. From the results, it can be seen that although the texture and radiation characteristics of SAR and optical images are different, the distribution of features actually has certain similarity in spatial scale, so that high accuracy tie points can be obtained in the case of features similar to houses as well as farmland. As shown in Figure 8, the registration between farmland and buildings can obtain good results, from which the pixel correspondence between optical images and SAR DOM can be effectively obtained, so that the latitude and longitude coordinates of the corresponding points can be read from the DOM, and finally used as PCPs to improve the planimetric positioning accuracy of the optical images after adjustment. At the same time, after selection and projection of the laser data, 9074 LAPs were obtained in the area of Jiangxi, and their distribution is shown in Figure 9a. The LAPs obtained after selection covers most of the image area and can provide enough LAPs as elevation control reference.
be obtained in the case of features similar to houses as well as farmland. As shown in Figure 8, the registration between farmland and buildings can obtain good results, from which the pixel correspondence between optical images and SAR DOM can be effectively obtained, so that the latitude and longitude coordinates of the corresponding points can be read from the DOM, and finally used as PCPs to improve the planimetric positioning accuracy of the optical images after adjustment. At the same time, after selection and projection of the laser data, 9074 LAPs were obtained in the area of Jiangxi, and their distribution is shown in Figure 9a. The LAPs obtained after selection covers most of the image area and can provide enough LAPs as elevation control reference.  After obtaining PCPs and LAPs as control points, in order to verify the effectiveness of adding LAPs and PCPs as control conditions to the adjustment and analyze the influence of the inconsistency of the object-image coordinates on the adjustment results, tests were conducted according to the following scheme. In addition, to study the influence of the plane position accuracy on the effect of the extracted LAPs, there are two options for acquiring the plane coordinate of the LAPs, namely, conducting adjustment with PCPs to enhance the plane positioning accuracy of the image before acquiring the LAPs and not using the PCPs before acquiring the LAPs.
Plan A: free network adjustment without PCPs or LAPs. Plan B: block adjustment with 14287 PCPs. Plan C: block adjustment with 9074 LAPs which were collected before adjustment with PCPs. After obtaining PCPs and LAPs as control points, in order to verify the effectiveness of adding LAPs and PCPs as control conditions to the adjustment and analyze the influence of the inconsistency of the object-image coordinates on the adjustment results, tests were conducted according to the following scheme. In addition, to study the influence of the plane position accuracy on the effect of the extracted LAPs, there are two options for acquiring the plane coordinate of the LAPs, namely, conducting adjustment with PCPs to enhance the plane positioning accuracy of the image before acquiring the LAPs and not using the PCPs before acquiring the LAPs.
Plan A: free network adjustment without PCPs or LAPs.  Table 2, while the distribution of the residuals of the results is shown in Figure 10. When the free network adjustment is performed without using any control points, the root mean square error (RMSE) in the east, north, plane, and elevation directions are 12.47 m, 7.64 m, 14.62 m, and 15.30 m, respectively. After adding the PCPs, the east, north, plain, and height accuracies change to 3.37 m, 6.28 m, 7.13 m, and 22.38 m, respectively. The plain accuracy is increased from 14.66 m to 7.13 m, an improvement of 51.2%, which indicates that the object coordinates of the PCPs are made to approximate their true planar positions in the adjustment iteration, and the PCPs can effectively improve the planar positioning accuracy after adjustment. Firstly, it can be seen that the elevation accuracy can be improved with the introduction of LAPs, regardless of whether the plane coordinates were corrected by adjustment with PCPs, from the original 15.30 m to 4.24 m without corrected plane accuracy and 3.09 m after correction, respectively, with 72.2% and 79.8% improvement. When LAPs collected before adjustment using PCPs are introduced as elevation control, the accuracies of the east, north, plane, and elevation directions are 7.31 m, 11.82 m, 13.90 m, and 4.24 m, respectively, while the adjustment results using LAPs collected after adjustment using PCPs are 6.85 m, 4.19 m, 9.02 m, and 3.09 m. It can also be found that, compared with the elevation accuracy after uncorrected plane coordinates, the introduction of LAPs as elevation control after corrected plane accuracy can reduce the actual elevation accuracy decrease caused by plane position error, and improve by 27.2% compared with the uncorrected before. Therefore, before the step of laser point introduction, it is recommended to use other control data to improve image plane positioning accuracy, which can effectively improve the improvement effect of laser point in elevation.
The adjustment results after using both LAPs and PCPs are 3.38 m, 6.30 m, 7.15 m, and 2.29 m in the east, north, plane, and elevation directions, respectively, which are comparable to the accuracy of 7.13 m obtained by using PCPs alone in plain, while the elevation accuracy is slightly improved compared to 3.09 m obtained by using LAPs alone. This indicates that the simultaneous use of elevation and plane control points can improve the accuracy in both plane and elevation directions, and can also maximize the use of control information provided by multi-source control data.  Firstly, it can be seen that the elevation accuracy can be improved with the introduction of LAPs, regardless of whether the plane coordinates were corrected by adjustment with PCPs, from the original 15.30 m to 4.24 m without corrected plane accuracy and 3.09 m after correction, respectively, with 72.2% and 79.8% improvement. When LAPs collected before adjustment using PCPs are introduced as elevation control, the accuracies of the east, north, plane, and elevation directions are 7.31 m, 11.82 m, 13.90 m, and 4.24 m, respectively, while the adjustment results using LAPs collected after adjustment using PCPs are 6.85 m, 4.19 m, 9.02 m, and 3.09 m. It can also be found that, compared with the elevation accuracy after uncorrected plane coordinates, the introduction of LAPs as elevation control after corrected plane accuracy can reduce the actual

Combined Block Adjustment Test for Mapping Satellite-1 Images
To further validate the method proposed in this paper, the Mapping Satellite-1 image test was conducted according to the same scheme in the Section 3.2, and a total of 13,213 qualified LAPs were obtained after laser point selection in the Hubei region, and their distribution is shown in Figure 9b. Using the extracted LAPs and PCPs, the block adjustment was performed, and the results of the ICPs after adjustment are shown in Table 3, while the distribution of the residual maps is shown in Figure 11.  The elevation accuracy was improved from 8.37 m to 3.56 m in case of using LAPs before adjustment with PCPs and to 3.51 m after adjustment with PCPs, which improved 57.3% and 58.1%, respectively. Whether the plane accuracy is improved or not before introduction of LAPs does not have much effect on the elevation accuracy improvement of Mapping Satellite-1 image after adjustment, probably because the resolution of Mapping From the results above, it can be seen that the accuracies of block adjustment without LAPs and PCPs are 5.94 m, 9.68 m, 11.35 m, and 8.37 m in east, north, plane, and elevation directions, respectively, and after introducing PCPs into adjustment first, the accuracies are 5.46 m, 7.14 m, 9.00 m, and 5.91 m in the east, north, plane, and elevation directions, respectively. Compared with the completely free network adjustment, the plane direction is improved from 11.35 m to 9.00 m, increased by 20.7%, and the result is the same as the result of ZY-3, which indicates that the plane control point can effectively improve the positioning accuracy in the plane direction. Meanwhile, the elevation accuracy is basically unchanged because the elevation direction is not constrained during the adjustment, and the elevation results of the intersection after each iteration are retained.
The elevation accuracy was improved from 8.37 m to 3.56 m in case of using LAPs before adjustment with PCPs and to 3.51 m after adjustment with PCPs, which improved 57.3% and 58.1%, respectively. Whether the plane accuracy is improved or not before introduction of LAPs does not have much effect on the elevation accuracy improvement of Mapping Satellite-1 image after adjustment, probably because the resolution of Mapping Satellite-1 image is lower compared with that of ZY-3, so the improvement of plane positioning accuracy does not have a significant effect on the elevation accuracy improvement. At the same time, the planimetric accuracy is partially improved compared with the free network adjustment, which reveals that the elevation accuracy can help improve the planimetric accuracy and make the ground point coordinates obtained by the stereo intersection model closer to the real ground point coordinates when the forward intersection is performed.
After using both LAPs and PCPs, the geometric positioning accuracies were improved to 5.45 m, 7.12 m, 9.00 m, and 3.51 m in east, north, planar, and elevation directions, which is basically the sum of the effects of using LAPs and PCPs separately, proving that using both control points can effectively improve the positioning accuracy of images after adjustment. The results are consistent with the conclusions drawn from the test results of ZY-3.

Discussion
From the experimental results, the following summary can be obtained.
(1) Using LAPs as elevation control can effectively improve the elevation positioning accuracy of the image after adjustment, and at the same time, if the accuracy of the plane position is firstly improved before introducing LAPs, then the plane error of the laser control points can be partially eliminated, thus further improving the results of its elevation positioning accuracy. The results of ZY-3 show that after adding LAPs, the elevation accuracy is improved from the original 15.30 m to 4.24 m without corrected plane accuracy and 3.09 m after corrected plane accuracy, which is 72.2% and 79.8% respectively, while the accuracy after corrected plane position is 27.2% higher than that without correction, which indicates that corrected plane positioning accuracy is effective for improving elevation accuracy. (2) The PCPs obtained by using SAR registration can effectively improve the planar positioning accuracy of the image after adjustment, but the accuracy depends on the sampling interval of the SAR DOM and the positioning accuracy, and if the resolution of the optical image is high and the resolution of the SAR DOM is enough high, taking the HOPC registration accuracy of 1.5 pixels as an example, a SAR DOM with a sampling interval of 5 m can provide a maximum positioning accuracy of 7.25 m for the optical image, which can effectively ensure the planar positioning accuracy for areas lacking GPCs in the global mapping process. However, for higher accuracy positioning requirements, it is necessary to provide a higher resolution and higher accuracy SAR DOM to meet the accuracy needs. This result shows that the combination of plan and elevation control points provides the best plain and elevation positioning accuracy.

Conclusions
This paper proposes a combined optical satellite stereo image adjustment method using spaceborne SAR image and laser altimetry data assistance, which first extracts the PCP by optical stereo image and SAR-DOM registration, and selection the LAPs from original laser points. Then, the observation equations of LAPs, PCPs, and tie points are established, respectively, and finally the orientation correction parameters of optical stereo images are obtained by iterative solution to improve the plane and elevation accuracy. The results of ZY3 using both LAPs and PCPs are 51.1% and 85.1% higher than those of free network adjustment in planar and elevation directions, respectively, reaching 7.15 m and 2.29 m, which can provide technical guarantee for stereo mapping of global large-scale satellite images without GCPs.