A Novel DEM Block Adjustment Method for Spaceborne InSAR Using Constraint Slices

The lack and uneven distribution of Ground Control Points (GCPs) will lead to the deterioration of Digital Elevation Model (DEM) block adjustment results in the bistatic Interferometric Synthetic Aperture Radar (InSAR) system. Given this issue, we first explain the relationship between the stability of adjustment parameters and the GCP distribution pattern theoretically using matrix perturbation theory. Second, we put forward the Constraint Slices (CSs) concept and first introduce CSs into the adjustment optimization model as constraint conditions rather than actual values as GCPs. Finally, we propose a novel DEM block adjustment method for spaceborne InSAR using CSs based on an optimization model with nonlinear constraints. The simulated experiment shows the instability of the conventional method and validates the proposed method under different parallel baseline errors. Four groups of real experiments were carried out according to the size of the uncontrolled area using twelve Co-registered Single-look Slant–range Complex (CoSSC) datasets for Henan Province, China. The adjustment results verified by the ICESat-2 ATL08 data demonstrate that the performance of the proposed method is better than the conventional method in the uncontrolled area; the corresponding improvements in adjustment accuracies compared with the conventional method are 0.13 m, 1.02 m, 2.12 m, and 8.18 m, respectively. At the same time, the proposed method can enhance the height consistency in overlapping areas, which is vital for seamless DEM production.


Introduction
The Bistatic Interferometric Synthetic Aperture Radar (InSAR) system [1] can effectively suppress the phase components caused by atmospheric delay and other factors to maintain the high coherence of interferograms, which is the leading technology to generate high-precision global DEM. TerraSAR-X/TanDEM-X is a typical example of the bistatic system, and has a baseline accuracy of about 1-2 mm [2,3].
The generation of interferometric DEM usually includes SAR image coregistration, interferometric processing, interferogram filtering, removing reference phases, phase unwrapping, elevation inversion, and geocoding [1,4]. The evaluation indicators of generated DEM include vertical accuracy and plane accuracy. Vertical accuracy is principally affected by the uncertainty of the spatial baseline, while plane accuracy is limited to the uncertainty of the parameters related to SAR imagery geolocation. Referring to the method of SAR bundle adjustment, J.Ma et al. [5] has proposed the Range-Doppler-Phase (RDP) model to correct key geometric parameters used in the R-D model and DEM inversion, including near range, imaging start time, orbit vector, and baseline. However, in recent years the geolocation accuracy of SAR satellites has improved from the hundreds of meters [6] to the decimeter order [7] with the development of satellite orbit determination and SAR calibration technology. For example, the orbital accuracy of TerraSAR-X is 20 cm, and the positioning accuracy can reach the order of decimeters [8,9]. Therefore, it is unnecessary to adjust geolocation parameters for spaceborne SAR systems with high positioning accuracy.
Because the SAR instrument calibration and the baseline determination cannot reduce all errors, there is a height error in the interferometric DEM caused by baseline uncertainty, instrument drifts, and other factors called the remaining systematic error [10]. Compared with the other factors, the parallel baseline error has a more significant impact on the accuracy of DEM generation, as follows: where R is the slant range, θ is the incidence angle, B ⊥ is the perpendicular baseline, and B err is the parallel baseline error. In order to correct the remaining systematic error, B.Wessel et al. put forward a method of DEM calibration to correct the systematic error of DEM products [11], which was the core principle of the Mosaicking and Calibration Processor (MCP) of TanDEM-X [12]. In this method, the systematic error of the generated DEM is expressed as a polynomial function of the image coordinates and the function model is solved employing the external high-precision Ground Control Points (GCPs), that is, the Global Land Surface Altimetry Data (GLAH14) of the Ice, Cloud, and land Elevation Satellite (ICESat) [13]. J. Hueso González et al. [14] focused on the quality assessment of ICESat GLAH14 elevation data and put forward basic and extreme selection criteria according to the number of echo peaks and bandwidth of measurement points. The characteristics of the systematic error were analyzed in [15], and the simulation data were used to illustrate the necessity and reliability of DEM calibration. Astrid Gruber et al. [10] described the concept of the Tie-Point chip and the selection of ICESat data, as well as validating the function model using external lidar DEM. GCPs are usually required to be sufficiently and evenly distributed to ensure the accuracy of DEM adjustment parameters [5,16,17]. Nevertheless, the GLAH14 data points are distributed along the ICESat orbit, which cannot meet these two requirements, especially the latter (the global distribution of ICESat data is shown on https://search.earthdata.nasa. gov/search?fp=ICESat&as [platform][0]=ICESat, accessed on 21 February 2022). Therefore, the adjustment results will deteriorate in certain cases and it is necessary to introduce external control data to ensure the accuracy of DEM adjustment parameters. However, this topic has seldom been discussed in previous studies on InSAR-DEM adjustment. The global public DEM can be regarded as reliable external control data. As an essential data source in geoscience, the global public DEM has high precision and wide coverage characteristics. The current global DEM/DSM can be divided into two categories according to the production technology; one is generated based on electromagnetic interference, such as Shuttle Radar Topography Mission (SRTM) [18] and the other is based on optical stereo pairs [19] such as ALOS World 3D (AW3D30) and the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM). In addition, the National Elevation Dataset project has provided LiDAR DEM covering almost the whole of the United States [20]. Many scholars have evaluated the accuracy of global public DEM/DSM. R.Talchabhadel et al. [21] used thousands of ground-based points near the basin of the West Rapti River, Nepal, to access the vertical accuracy of AW3D30 and SRTM. The AW3D30 showed the highest accuracy of 3.41 m in the lower region of the study area. H.Li et al. [22] used ICESat data, SRTM, and ASTER GDEM to evaluate AW3D30 comprehensively. Their results showed that the horizontal offset of AW3D30 is approximately 5.0 m and 10.0 m, respectively, relative to SRTM and ASTER GDEM2, while it has a higher vertical accuracy of 4.81 m than the others. However, it should be noted that the vertical accuracy of public DEM is generally less than that of GLAH14 laser altimetry data, which means that the public DEM cannot be directly regarded as GCPs in adjustment as otherwise it will affect the correction effect of GLAH14 data on the generated DEMs. Therefore, it is necessary to develop an appropriate form to introduce the global public DEM into DEM block adjustment in order to solve the deterioration caused by the lack or uneven distribution of GCPs.
Based on the above analysis, this paper first analyzes the influence of ICESat distribution on adjustment results using matrix perturbation theory. Then, we propose the concept of constraint slices (CS) and the corresponding extraction method from external public DEM. Finally, taking the relative accuracy of CSs as constraint conditions, we propose a DEM block adjustment method assisted by CSs. A simulated experiment was designed according to different spatial baseline errors, and the overall adjustment accuracies of the simulated experiment improved by 0.41 m to 7.01 m. Four groups of real data experiments were designed according to the number of GCPs, and the overall adjustment accuracies improved by 0.02 m to 7.86 m.

Conventional Method and Shortcomings
The conventional method [10,11,14] models the remaining systematic height error of InSAR DEMs as a polynomial function and uses Tie Point (TP) chips and GCPs for joint correction.

TP Chips
In the overlapping area of two DEMs, pixels corresponding to the same ground object are called tie points [11]. This concept is similar to the connections/tie points in SAR imagery block adjustment [23], though their extraction methods are different. Gaussian white noise usually contaminates the generated DEM, and the pixels near the edge may even be affected by colored noise [11]. To ensure the relative accuracy and reliability of connection points, Huber M. et al. put forward the concept of TP chips and developed a set of extraction processes [24]. First, the image coherence threshold is set to eliminate the pixels with low coherence. At the same time, the areas with higher height differences compared with the public external DEM are removed to eliminate the invalid areas, such as shadow, water, and overlap areas. Then, a series of grids called TP chips are delimited along the longitude and latitude direction in the overlapping area. A median height value is assigned to a TP chip after calculating the histogram of effective pixels in each grid. The geographic coordinate of a TP chip is that of the center point in the grid. Compared with feature matching, this method is more time-saving [10,11,24]. It is able to realize time savings of more than 95.0% thanks to the extraction of TP chips mainly using logical operations such as establishing grids and masking invalid areas. The extraction diagram of TPs is shown in Figure 1.

GCPs
ICESat data provide globally distributed accurate elevation measuring points, and has been used as an absolute height reference for DEM adjustment [10,11]. ICESat data are measured with lasers, which means the measurement points are affected by clouds, vegetation, and other factors and cannot reflect the actual ground elevation. Therefore, J.H. González et al. [14] focused on ICESat data screening in DEM adjustment and put forward extreme and basic selection criteria. The specific process is as follows. First, the ICESat points with significant height differences (greater than 200 m) are eliminated through comparison with external DEM/DSM. Then, data points with one peak for the reflected signal and a time bandwidth of fewer than 3.2 ns (extreme selection criteria) or 8.0 ns (basic selection criteria) [14] are screened to ensure a relatively concentrated energy distribution to the greatest possible extent. Finally, all DEM pixels in the footprint of the ICESat point are averaged using two-dimensional Gaussian weights.

Function Model
According to [10,15], the systematic error of generated DEM can be expressed as a polynomial of image coordinates, as in Equation (2) g I (x, y) = a I + b I x + c I y (2) where I is the index of DEM image, a I is the height offset, and b I and c I are the tilts in ground distance and azimuth, respectively. x, y are the image coordinates, i.e., ground distance and azimuth with respect to a reference point. The function model can adopt higher-order polynomials, although Runge's phenomenon, where the model fluctuates violently at the boundary of the adjustment region, may need to be taken into account. For TPs and GCPs, the observation equations [10,11] are where t and p are the indexes of TPs and GCPs, h is the height extracted from the generated DEM, H p is the height of GCP, and l is the residual. Equation (3) is the error equation of a TP chip and represents the elevation difference of different DEMs in the overlapping area. Equation (4) is the error equation of a GCP, and represents the elevation difference between DEM and absolute elevation reference data. Equations (3) and (4) can be further sorted as follows: For multiple GCPs and TP chips, Equations (5) and (6) can be further sorted into the matrix form as Equation (7): where A P = diag(A 1 P , · · · , A I P , · · · , A N P ) are the coefficient matrices of TP chips and GCPs, respectively, x = (a 1 , b 1 , c 1 , ..., ..., a N , b N , c N ) T is the vector of unknown parameters, N is the number of generated DEMs, I ,J, and K are indexes of generated DEMs, respectively, ∆H is the vector of height differences, and L is the residual vector. The least square solution of Equation (7) is

Shortcomings of the Function Model
The distribution and number of GCPs are the key factors affecting the adjustment results. On the one hand, the distribution of GCPs should be even for the adjustment results to fully reflect the variation in the elevation error with two-dimensional coordinates for ground distance and azimuth. On the other hand, all the generated DEMs involved in the adjustment should contain enough GCPs to prevent the errors of adjacent DEMs from spreading to the whole area through the TPs, leading to deterioration of the adjustment results in the uncontrolled area. In 2011, Ma Jing [5] discussed the influence of the distribution of GCPs on airborne InSAR adjustments. She explained that when the distribution of GCPs is relatively concentrated, the error can be transmitted to the uncontrolled area, resulting in deterioration of the adjustment results.
However, GLAH14 data have difficulty meeting the above two requirements, which is not conducive to obtaining reliable and stable adjustment results. Section 3 shows the distribution of GLAH14 data. It can be seen that the distribution pattern of ICESat points is a line along the satellite's flight trajectory. Conversely, ICESat data does not meet the distribution requirement of GCPs, such as rectangular shape or cross shape in adjustment, as shown in Figure 2. On the other hand, the significant distance between ICESat orbits causes certain images to be adjusted to contain few or no GCPs. Appendix A theoretically analyzes the influence of linear GCP distribution patterns on large-scale DEM block adjustment. First, we demonstrate that all eigenvalues of the coefficient matrix A should be large enough to obtain a stable solution x. Then, we prove that when all GCPs concentrate near a straight line, A has at least a relatively small eigenvalue, which will lead to the instability of x. Here, we take the calculation of the ground tilt b as an example, as shown in Figure 3, to provide a visual explanation for the analysis in Appendix A. The estimated value is b = ∆h err /∆x with standard deviation σ b = σ ∆h err /∆x, where σ ∆h is the uncertainty of the height error and ∆x is the distance between two different GCPs. It can be seen that σ b increases with the decrease of ∆x, which means that when the GCPs concentrate at a certain ground distance, that is, ∆x is small, the estimated value of b will fluctuate greatly. For example, assuming that the swath width of generated DEM is 30 km and σ ∆h err is 0.5 m, and assuming that the uncertainty of ground distance tilt b is less than 0.6 m within the width of 30 km, then ∆x = σ ∆h err /σ b = 0.5/(0.6/30,000) = 25 km, which means that the distance between two GCPs in the direction of ground distance should be greater than 25 km. However, ICESat data cannot always meet this strict condition, where the yellow dots show the location of GLAH14 data). It can be seen that the ICESat data for the adjustment area are concentrated in the middle, around the red solid frame area. The northwest area (i.e., the red dotted frame area) is a large uncontrolled area. Estimation of range tilt. The dotted line represents the actual terrain, the solid line represents the generated DEM, and the longitudinal axis represents the height difference between the generated DEM and the actual terrain; x 1 and x 2 are the ground distance coordinates of two different GCPs, h err,1 and h err,2 are their systematic elevation errors, and ∆h err = h err,2 − h err,1 and ∆x = x 2 − x 1 .

Constraint Slice (CS) and Constraint Conditions
As an essential data source in geoscience the global public DEM have high precision and wide coverage characteristics, as represented by the Shuttle Radar Topography Mission (SRTM) and ALOS World 3D (AW3D30). To solve the problem mentioned in Section 2.1, we tried to introduce the external public global DEM into DEM block adjustment. It should be noted that the vertical accuracy of global public DEM is often far less than that of GLAH14 data and generated DEM. For example, the nominal vertical accuracy of ALOS DEM is 5 m [25] and therefore it cannot be directly added to the adjustment process in the form of GCP.
As discussed in Section 2.1.4, it is difficult to ensure the reliability and stability of adjustment results in the conventional model when GCPs are linearly distributed. Therefore, we tried to overcome this weakness by introducing external terrain data.
Based on the above analysis, we propose a DEM block adjustment method assisted by Constraint Slices (CSs). This method improves the conventional model via two approaches, namely, the concept and extraction approach of Constraint Slices (CSs) and the addition of the elevation variances of CSs to the adjustment optimization model as constraint conditions. The specific process is described as follows.

The Concept and Extraction of Constraint Slices
CSs are a form of external DEM participating in DEM adjustment. CSs are added to the adjustment optimization model as nonlinear constraint conditions and are extracted as slices, inspired by the TP chips described in Section 2.1.1, which represent the two crucial differences in CSs compared with ordinary GCPs and TPs. A flowchart of the process of extracting CSs is provided in Figure 4. The specific steps are as follows: i.
Geocode the external DEM into the SAR coordinate system (i.e., forward geocoding). This step is usually performed after phase unwrapping because this results in an integer multiple of 2π between the absolute phase and the unwrapped phase [4], which needs to be corrected by external DEM. ii. Generate the slope map of the reference DEM in the SAR coordinate system. Wang [26] provides a method for generating a slope map using DEM through calculating the gradients in the range and azimuth directions ∂h/∂rg, ∂h/∂az using the Sobel operator and then converting them into a slope map according to Equation (9).
iii. Calculate the elevation differences between the generated and external DEM, and set the pixels with large elevation differences (e.g., greater than 50 m) as invalid pixels, including phase unwrapping error areas and low coherence areas in generated DEMs or elevation anomaly area in external DEMs. iv. Divide the geocoded DEM into several blocks at an interval of about 1 km. There are two sizes of TP chips in the existing literature, 1 km [10,24] and 500 m [11]. However, because the resolution of the external public DEM is often lower than that of CoSSC data, we assume that a grid can contain more elevation pixels, ensuring the reliability of the CSs. Therefore, inspired by the concept of TP chips and combined with the characteristics of public DEMs, we set the size of the CSs to 1 km. Then, calculate the histogram of each block in the generated and reference DEMs taking the median values of the histograms as their elevations and assign the center coordinates to each block. This part draws on lessons from the extraction method of TPs in Section 2.1.1 in order to avoid the influence of noise elevation outliers. v.
Calculate the average slope of each block and divide the CSs into flat and mountainous areas, as the public DEM often shows different elevation accuracies in mountainous areas and plain areas [21,22]. When the experimental area contains complex terrain, considering only the flat area with the highest accuracy while ignoring the mountainous area will lead to insufficient CSs.

Transformation of Adjustment Model
Taking the relative height accuracy of CSs as nonlinear constraints, we transform the adjustment optimization model into Equation (10): where f and m denote CSs of the flat and mountainous region, respectively, σ f and σ m are the upper limits of the residual standard deviation (which can be empirically obtained from previous studies [18,21,22] or calculated using the GLAH14 data of the experimental area according to the DEM accuracy evaluation method [21,22,27,28]), A f and A m are the coefficient matrices of the CSs, and ∆H f and ∆H m are the elevation differences between the CSs and generated DEMs. In Appendix B, we use the Lagrange multiplier method [29] to deduce the solutions of the improved model.
It should be noted that the RMSEs of CSs, i.e., A f x − ∆H f 2 2 and A m x − ∆H m 2 2 , are not used as constraints here. This is because RMSE represents the absolute height difference between the generated DEM and the external DEM, and using RMSE means taking the external DEM as absolute elevation reference data; in this way, the external DEM is equivalent to the GCPs, as per Equation (4). However, as described in Section 1, the vertical accuracy of the public DEM is usually worse than GLAH14 data, and using RMSEs relative to the external DEM as constraints thus results in lost adjustment accuracy.  [14]. The location relationship between GLAH14 and CoSSC data is shown in Figure 5, in which the serial numbers are arranged according to their acquisition time. The figure shows that the eighth and ninth images do not contain any GCPs, and are thus called the uncontrolled area, while the distribution of GCPs in the other images is uneven. The experimental region contains various geomorphic features with a maximum height difference of about 1400 m. In the southeast is a large plains region, to the northwest are the Songshan Mountains, and to the north lies the Yellow River.

Experimental Data
The reference data used were ICESat-2 Land and Vegetation Height data (ATL08), whit a vertical uncertainty of 0.2 m for flat terrain and 2 m for mountainous terrain [30].

Simulated Data
The purpose of this experiment was to illustrate that (1) the traditional methods can lead to inaccurate estimation of adjustment parameters in the case of uneven distribution and insufficient quantity of ICESat data, and (2) the improved method can effectively stabilize the estimated parameters. The central idea of the experiment was to simulate multiple groups of InSAR DEM and compare the performance of the conventional model and proposed model under different baseline errors. The overall flow chart of the simulated experiment is shown as Figure 6, and the specific process is described as follows: i Simulate the terrain of the experimental area using fractal theory [1], as in Figure 7. ii Extract ICESat elevation data and generate external DEM. We used the actual geographical coordinates of the ICESat points to extract elevation data from the simulated terrain in step i and added 0.5 m Gaussian noise. Similarly, we added 5.0 m Gaussian noise to the actual terrain as the external DEM. iii Geocode the simulated terrain. The simulated terrain is geocoded into the slant-range coordinate system using the real satellite orbit information. The geocoded DEM data are used to generate absolute interferometric phase and as the test set of adjustment results. iv Generate the absolute terrain phase using actual orbit information and the geocoded data from Step iii. v Carry out InSAR processing, including removing the reference plane phase and generating DEMs. To the master and slave satellite orbit data used in this step were added baseline errors ranging from 0.5 mm to 5.0 mm. In addition, 1.0 m Gaussian noise was added to the generated DEMs. The simulated baseline and random height errors were set based on the actual situation. At present, the baseline accuracy of the bistatic InSAR system is usually on the order of millimeters [2,3]. The published global TanDEM-X DEM has a vertical accuracy of about 1.0 m [31]. Therefore, we chose 0.5 to 5 mm as the baseline error range and 1.0 m as the random noise of the simulated generated DEM. vi Carry out DEM block adjustment using the conventional and improved method. In the improved method, the CSs are extracted from the external DEM in Step ii. In addition, there is no distinction between mountains and flat land when adding noise in Step ii; thus, the CSs are not distinguished. vii Accuracy check; checkpoints are evenly selected from the geocoded data from Step iii.
To verify the effectiveness of the proposed method under different degrees of baseline error, we gradually increased the parallel baseline error at intervals of 0.5 mm to generate ten groups of simulation data. We used the root-mean-square error (RMSE) to measure the adjustment accuracy, as follows: where h DEM and h check are the elevations of the generated DEMs and check points, respectively.   The RMSEs of ten sets of data are shown in Figure 8, where Figure 8a-c are the RMSEs of all twelve DEMs, controlled DEMs, and uncontrolled DEMs, respectively. The red, blue, and green lines represent the height RMSEs of the generated DEMs, the adjustment results using the conventional method (7), and the adjustment results using the proposed method (10), respectively. Before adjustment, the elevation errors of the generated DEM are approximately proportional to the baseline errors, which is consistent with Equation (1). Overall, as shown in Figure 8a, the vertical accuracy of traditional methods is unstable and sometimes worse than the original DEM. The RMSEs of the improved method are stable at about 1.0 m, which is very close to the energy of the manually added Gaussian noise in step iii, demonstrating that the improved method nearly eliminates all the systematic errors. In the controlled area, as shown in Figure 8b, both methods show good performance; however, the proposed model is slightly better than the function model. In the uncontrolled area, as in Figure 8c, the adjustment results of traditional methods cannot reflect the actual systematic height error. In addition, the changing trend of RMSEs in Figure 8c is consistent with that in Figure 8a, which confirms that the deterioration of vertical accuracy in uncontrolled areas seriously affects the overall adjustment accuracy. After being constrained by the external DEM, the RMSEs are stable at about 1.0 m. All in all, the experimental results show that after adding constraint slices the overall adjustment accuracy is significantly improved, which is due to the stability of the adjustment results in the uncontrolled area on the one hand and the improvement of the adjustment accuracy in the controlled area on the other. To further explain these two reasons, we show the range tilt b and azimuth tilt c of two different orbits in Figure 9. Figure 9a,b show the mean coefficients of the track containing the eighth to ninth CoSSC data (the uncontrolled area), and Figure 9c,d show the mean coefficients of the track containing the first to third CoSSC data (the controlled area). The red, blue, and green lines are the actual values of b or c, the estimated values of the conventional method, and the estimated values of the proposed method, respectively. As there are no GCPs, the coefficients of the eighth and ninth images fluctuate considerably without constraint slices, as shown in the green lines in Figure 9a,b. The RMSEs of the estimated values are 3.08 × 10 −4 and 5.07 × 10 −5 , respectively. However, after introducing constraint slices the adjustment results are stable (shown as blue lines in Figure 9a,b) and the RMSEs of the estimated values are 2.33 × 10 −5 and 6.73 × 10 −6 , respectively, which is an order of magnitude lower than the conventional method. The first, second, and third images contain enough GCPs; however, they are distributed unevenly, especially in the ground distance, as shown in Figure 5, which directly results in the offsets of the range and azimuth tilts in Figure 9c,d. The RMSEs of the estimated values using the conventional method (7) are 3.44 × 10 −5 and 8.80 × 10 −6 , respectively. After introducing constraint slices, the RMSEs decrease to 1.03 × 10 −5 and 6.89 × 10 −6 .

Real Data
To verify the effectiveness of the proposed method in the actual situation, we changed the size of the uncontrolled area and designed four groups of data experiments. Table 1 shows the data contained in the uncontrolled area in each group of experiments. We generated InSAR DEMs in order to mosaic the geocoded DEMs. Figure 10 displays the geocoded height map of the original DEMs before adjustment. The Yellow River passes through the north of the experimental area, and the phase unwrapping in this area is therefore wrong, resulting in abnormal elevation values. The CSs were extracted from AW3D30 DEM according to the extraction process described in Section 2.2. Figure 11 shows the distribution and slope of the constraint slices in terms of slant-range coordinates. The color represents the slope of CSs. We set the slope threshold to 10 • in order to distinguish mountain CSs and flat CSs. The ICESat-2 ATL08 data were used as check points to verify the accuracy of adjustment results. Unfortunately, the tenth image does not contain any proper check points. Figure 12 shows RMSEs of all images. Table ?? shows RMSEs of the conventional and proposed methods in the controlled, uncontrolled, and general experimental areas.

Performance of Conventional Method in Uncontrolled Areas
The conventional method can maintain a stable adjustment effect when the uncontrolled area is small. Taking Ex.1 as an example, the elevation accuracy of the uncontrolled area before adjustment is 3.08 m and is improved to 2.12 m after adjustment by the conventional method. With the gradual increase of uncontrolled areas, the adjustment results of the conventional method become worse. Taking Ex.4 as an example, prior to adjustment the elevation accuracy of the uncontrolled area is 2.37 m, while after adjustment RMSE deteriorates to 10.41 m.

Performance of Conventional Method in Controlled Areas
The conventional method shows excellent adjustment performance in the controlled area. Taking the first scene data as an example, as shown in Figure 12, the RMSE is 1.96 m before adjustment. After adjustment with the conventional method, the vertical accuracy of the four groups of experiments is stable at about 1.53 m, about 0.4 m better than the original RMSE.

Performance of the Improved Method in Uncontrolled Areas
From Figure 12, the overall adjustment results of the improved method are stable. In Ex.1, the RMSE of the generated DEMs in the uncontrolled area is 3.08 m. The proposed method using CSs brings the RMSE to 1.99 m, and the overall accuracy is improved by about 1.0 m. In Ex.4, before adjustment the vertical accuracy of the uncontrolled area is 2.37 m. After adjustment by the proposed method, the RMSE decreases to 1.96 m, about 8.4 m better than the conventional method. This shows that the improved method can effectively stabilize the adjustment results in the uncontrolled area and improve the vertical accuracy of DEM products.

Performance of the Improved Method in Controlled Areas
In the four experiments, the conventional method improves the vertical accuracies of controlled areas by 0.35 m, 0.20 m, 0.17 m, and 0.43 m, respectively, while the proposed method improves the vertical accuracy by 0.37 m, 0.11 m, 0.20 m, and 0.43 m respectively. Therefore, the performance of the proposed method is similar to the conventional method in the controlled area.
In Table ??, the vertical error of the generated DEM is 2.32 m. The adjustment result of the conventional method is greatly affected by the uncontrolled area. With the increase of the uncontrolled area, the RMSE of the traditional method deteriorates from 1.91 m to 9.78 m. The adjustment results of the proposed method are stable at around 2.0 m, about 0.3 m higher than the original height error. All in all, the data experiment shows that the improved method significantly improves the adjustment results in uncontrolled areas and maintains the adjustment results in controlled areas at the level of the conventional method; thus, the adjustment results tend to be stable and reliable. Figure 13 shows the geocoded DEM after adjustment where Figure 13a-d are the results of the traditional method and Figure 13e-h are adjusted by the proposed method. It can be seen that the deterioration in the adjustment results in uncontrolled areas leads to significant elevation differences in overlapping areas, resulting in the blocky effect in the geocoded DEM shown in Figure 13b-d. The proposed method can sufficiently suppress the elevation difference, mainly because the public DEM used in the proposed method is seamless. Figure 14 shows the height difference map of Figure 13a,e compared with AW3D30. The suppression effect of the proposed method on elevation difference is reflected significantly in the marked area.
However, our method does has disadvantages. The core idea of our method is to carry out DEM adjustment by introducing global public DEM as auxiliary data. Nevertheless, the measurement time of the external DEM is not consistent with that of the InSAR data. Therefore, elevation changes will inevitably occur when there is a significant time interval between the acquisition of the public DEMs and the CoSSC data, and the resulting CSs then include the systematic elevation error and changes, affecting the final adjustment accuracy. In addition, certain DEMs may have systematic errors, which can affect the accuracy of the CSs. For example, research [25] has pointed out that AW3D30 suffers from the mismatch between adjacent scenes and horizontal and oblique stripping for certain areas. Consequently, to resolve this problem we can introduce multi-source public elevation data, such as SRTM, GDEM, or LiDAR DEM, to develop a robust CS extraction scheme considering data acquisition time and external DEM accuracy. For example, we can assign weights to various external public DEMs according to their acquisition time and obtain the weighted elevation to ensure the reliability of CSs.

Conclusions
This paper proposed a DEM block adjustment method for spaceborne InSAR using Constraint Slices to solve the deterioration of DEM adjustment results caused by the lack or uneven distribution of GCPs. We selected Constraint Slices from the external public DEM and introduced them into the conventional DEM calibration method in the form of nonlinear constraints. A simulation experiment explained the influence of the lack and uneven distribution of GCPs on the adjustment results and proved that the proposed method could obtain more stable and accurate adjustment results under different degrees of baseline errors. The real data experiment with twelve CoSSC datasets demonstrated that the performance of the proposed method is better than the conventional method in uncontrolled areas and can enhance the elevation consistency of InSAR DEMs in overlapping areas. In the simulated experiment, the overall adjustment accuracy was increased by 0.41 m to 7.01 m. Compared with the conventional method, the adjustment accuracies of the four real data experiments in the uncontrolled area were improved by 0.13 m, 1.02 m, 2.12 m, and 8.45 m, respectively, and the overall adjustment accuracies were improved by 0.02 m, 0.25 m, 1.01 m and 7.86 m, respectively. However, our method still has disadvantages. In the future, it might be possible to introduce multi-source public elevation data in DEM block adjustment to develop a robust CS extraction scheme that considers data acquisition time and external DEM accuracy.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Influence of GCPs Distribution on Eigenvalues of Coefficient Matrix in the Conventional Model
First, we use the Singular Value Decomposition (SVD) method to reveal the relationship between the solution (8) and the eigenvalues of coefficient matrix A. Assuming that the coefficient matrix A is accurate, the uncertainty of observation vector ∆H is δH and the SVD of A ∈ R m×n is A = UΣV T , whereΣ = [Σ 1 , O] T and Σ 1 = diag{σ 1 , σ 2 , · · · , σ n } = diag{λ 2 1 , λ 2 2 , · · · , λ 2 n }, σ are the singular values of A and λ represents the eigenvalues of A T A. U and V are the Unitary matrices. Then, the error of x caused by δH is where Σ = [Σ −1 1 , O] T and Σ −1 1 = diag{1/λ 2 1 , 1/λ 2 2 , 1/λ 2 3 }; (A1) shows that the stability of solution (8) completely depends on the singular values of A. If A has a relatively small eigenvalue, or even an eigenvalue approximately equal to zero, there will be a large error in solution (8).
Assuming that the Jth generated DEM contains n GCPs located on the same satellite trajectory and that the geographical coordinates of each GCP are (x i , y i ), then the coefficient matrix is and According to the distribution characteristics of ICESat points, it can be assumed that the geographical coordinates x i and y i of GCPs are approximately distributed on the same straight line in each generated DEM, i.e., where k and b are the slope and intercept of the fitted line, respectively, y i is the fitting value corresponding to x i , and ε i is the offset between real coordinates and the fitted line.
In fact, as shown by the yellow points in Figure 5, ε i is tiny relative to x i and y i because of the stable satellite orbit. Replace y i with Equation (A4), and Equation (A3) can be expanded as Equation (A5), where δB J is called perturbation of B J in perturbation theory [32]. It can be proved that B J is rank-deficient, that is, B J has at least one zero eigenvalue.
Suppose the eigenvalues of A T J A J are {σ * J,1 , σ * J,2 , σ * J,3 }, σ * J,1 ≥ σ * J,2 ≥ σ * J,3 ≥ 0 and the eigenvalues of B J are {σ J,1 , σ J,2 , σ J,3 }, σ J,1 ≥ σ J,2 ≥ σ J,3 = 0; then, where tr(A T J A J ) is the trace of A T J A J . Then, according to the Weyl monotonicity theorem [33,34], where δB J spec is the spectrum norm of B J , that is, the maximum eigenvalue of B J , as Equation (A8): In order to make the analysis more clear, we assume that ε is a random variable with zero expectation and independent of the geographical coordinates of ICESat points (which, in fact, it is). Then, Equations (A6) and (A8) are arranged into statistical form as where E(·) is an expectation operator. Then, combining Equations (A9) and (A10), we can obtain σ *

Appendix B. Solution of the Constraint Model
The Lagrange function corresponding to the constraint model Equation (10) is where A f ∈ R n 1 ×3 and A m ∈ R n 2 ×3 n 1 and n 2 represent the numbers of flat and mountain CSs, respectively. Var(·) is the variance operator and can be expanded as where I n×1 is an all-1 vector, 1 n×n is an all-1 matrix, and m f = (A f x − ∆H f ) T I n 1 ×1 is the mean value of vector A f x − ∆H f . Var(A m x − ∆H m ) can be expanded as Equation (A13). Then, the gradient function of Equation (A12) with respect to x is For simplicity, we replace (E n 1 − 1 n 1 I n 1 ×n 1 )A f and (E n 2 − 1 n 2 I n 2 ×n 2 )A m in Equation (A14) with A f and A m , respectively; E is an unit matrix.
Finally, we list the Karush-Kuhn-Tucker conditions (KKT conditions) [29], as follows: which can be solved iteratively by classical constrained optimization algorithms such as the Interior Point algorithm [35]. The iteration result of the kth step is