A Novel Three-Dimensional Block Adjustment Method for Spaceborne InSAR-DEM Based on General Models

Digital elevation model (DEM) generated by bistatic synthetic aperture radar interferometry (InSAR) usually has systematic plane and elevation errors. However, the DEM adjustment method based on the function model can only correct systematic elevation error, affecting the plane accuracy of DEM mosaic products. Given this issue, this article proposes a 3-D InSAR-DEM block adjustment method based on general models with a separate adjustment strategy, that is, first correct vertical error using the adjustment method based on the function model and then correct horizontal error using the method based on rational function model after compensating for elevation errors. In order to extract the tie points (TPs) and ground control points (GCPs) used in adjustment, we propose a normalized correlation coefficient matching method based on complex slope maps (NCC-CSM). We conduct experiments using 29 TanDEM-X coregistered single look slant range complex data. The accuracies of adjustment results are verified using ICESat-2 ATL08 data and selected GCPs, respectively. The adjustment experiments show that the elevation accuracy of InSAR-DEM is improved from 1.70 to 1.35 m, the absolute and relative plane accuracies are improved by 5.58 and 11.53 m, respectively, and the horizontal and vertical geometric discontinuity in the overlapping area disappear. The matching experiments of TPs and GCPs show that the proposed NCC-CSM method is superior to the traditional NCC method based on DEM (NCC-DEM), especially in flat areas.

Due to the uncertainty of satellite orbit, space baseline, and slant range, the DEM generated by InSAR inevitably has systematic errors, including vertical and horizontal errors. The vertical error is mainly caused by the parallel baseline error [5], which is expressed as the tilt of the elevation surface along with the ground distance [6], [7]. The horizontal error is caused by the errors of geolocation parameters [8], including satellite orbit, slant range, imaging time, and elevations of objects.
Accurate positioning accuracy is an essential prerequisite for DEM application. The block adjustment method can correct the systematic errors of InSAR-DEM and improve the relative and absolute accuracies. According to the error model, the existing InSAR-DEM adjustment methods can be divided into two categories: 1) the methods based on the rigorous [9], [10], [11] and 2) general model [6], [7], [12]. The former refers to the adjustment method based on the "Range-Doppler-phase" model [10]. Its core idea is to correct the geometric parameters of the InSAR system, including baseline length and inclination, so it can also be called a parameter-level adjustment method. Currently, this method is often used in the joint calibration of airborne InSAR systems [13]. The latter refers to the DEM adjustment method based on the function model [6]. By analyzing the cause mechanism, this method models the systematic height error of interferometric DEMs as a polynomial relationship about image coordinates and solves the adjustment parameters through tie point (TP) chips [14] and height control points (HCPs) [15]. Because of the simplicity of polynomials, after adjustment, the function model can directly be used to compensate for DEM elevation error. It has been applied in the mosaicking and calibration processor [16] of TanDEM-X, so it can also be called the product-level adjustment method. However, the existing adjustment methods still have some problems. First, after parameterlevel adjustment, it is necessary to use the corrected system parameters to regenerate DEM. This step is time-consuming and unsuitable for the batch DEM production process of the spaceborne bistatic InSAR system. Second, the product-level adjustment method does not integrate the vertical and horizontal error, which is not conducive to the seamless DEM [17] production.
The errors of geolocation-related parameters, such as slant range and imaging time, can affect the horizontal accuracy of geocoded DEM [18]. Generally speaking, there is still a lack of a general model-based 3-D adjustment method for Spaceborne InSAR-DEM.
From the perspective of processing flow, block adjustment technology consists of two steps: 1) control data extraction and 2) adjustment model construction. The control data can constrain remote sensing products' absolute and relative accuracies. Ground control points (GCPs) and TPs are their most widely used forms. For large-scale area adjustment technology, image matching is crucial for acquiring TPs and GCPs. The existing DEM matching methods can be roughly divided into two categories. One is the matching method based on the minimum height difference, i.e., Least Z-difference (LZD) [19], [20]. LZD method models the offset between DEMs as seven parameters, including 3-D rotation, translation, and elevation scaling, and uses mean square elevation difference as a loss function to solve these parameters iteratively. The other is the method based on digital image processing (DIP), such as the mutual information (MI) [21] and the normalized cross-correlation (NCC) method using DEMs (NCC-DEM) [22]. The two methods take the DEM mutual information and cross-correlation coefficient as the measurement function, respectively, and find the peak value of the measurement function to determine the offset of DEMs. However, the existing methods still have three problems. First, the calculation is complex. For example, the MI method needs to calculate the joint histogram between two images. Second, the matching effectiveness becomes worse in areas where the topographic relief is not apparent, such as the NCC-DEM method. Third, these matching methods does not set up a reliability measurement criterion to eliminate mismatched points, which is a necessary step in InSAR-DEM block adjustment. Therefore, it is necessary to develop an effective DEM matching method for InSAR-DEM block adjustment.
In summary, it is necessary to develop a 3-D adjustment method and a DEM matching method to extract control data for InSAR-DEM adjustment. In response to the above problems, this article makes the following innovations. First, we propose a 3-D InSAR-DEM block adjustment method based on general models with a separate adjustment strategy. The core idea of this method is to carry out elevation adjustment first using the method based on the function model and then execute plane adjustment using the method based on rational function model (RFM) after correcting the systematic height errors. Second, we develop a matching method based on complex slope maps (CSMs), which contain the slope and aspect of DEMs, to extract the TPs and GCPs used in adjustment. The method realizes DEM matching by calculating CSM's normalized correlation coefficient maps. Third, we develop the corresponding screening criteria according to the characteristics of correlation coefficient maps to select TPs and GCPs. We use 29 TanDEM-X Coregistered single-look slant-range complex (CoSSC) data covering Henan province, China, to verify the proposed method. The adjustment results show that the absolute and relative elevation accuracies of DEM have been improved by 0.35 and 1.19 m, respectively. The absolute and relative plane accuracies have been improved by 5.58 Fig. 1. Plane error caused by elevation error. R, h, and θ are the slant range, elevation, and incidence angle of a target, respectively. Δh is the height error, and Δx is the geolocation error caused by Δh. and 11.53 m, respectively. The geometric discontinuity between different scenes has disappeared. Matching experiments show that the performance of the proposed normalized correlation coefficient matching method based on the complex slope maps (NCC-CSM) method is better than the traditional NCC-DEM method, especially in flat areas.

II. ADJUSTMENT STRATEGY AND MODEL
In this section, we introduce the strategy, model, and process of the proposed InSAR-DEM adjustment method, where the separate adjustment strategy is named as "first elevation and then plane," and the adjustment model is based on general models, including the function and RFM model.

A. Adjustment Strategy
The core idea of the "first elevation and then plane" adjustment strategy proposed in this article is to correct the vertical errors of InSAR-DEM and then correct their horizontal errors. The adjustment strategy is proposed primarily for the following reasons. First, the generation mechanism of elevation and plane errors of InSAR-DEM are different. The elevation error is mainly caused by parallel baseline error, while the plane error is affected by other geometric parameters, such as slant range. Second, due to the side-view imaging characteristics of SAR systems, the elevation error will affect the plane accuracy of DEMs, the mechanism of which is shown in Fig. 1. The plane error Δx caused by elevation error Δh can be expressed as Therefore, carrying out the elevation adjustment first can avoid the impact of elevation error on the plane geolocation accuracy. Third, globally distributed control data have different characteristics. Laser altimetry data of Ice, Cloud, and land Elevation Satellite (ICESat) is often used as elevation control data, and can provide height control points better than 1.5 m. However, it does not have visual image information, so it cannot be used as horizontal control data. In contrast, the global public DEM is also commonly used as the control data in adjustment. Although it has visual image information, its elevation accuracy is usually inferior to InSAR-DEMs. Therefore, the separate correction of elevation and plane can give full play to the respective advantages of global control data. The proposed adjustment strategy includes two steps: 1) elevation adjustment and 2) plane adjustment, which are connected through systematic elevation error correction. The specific process are are as follows.
1) Select TP chips and HCPs. Here, the TP chips refer to the grids in the InSAR-DEM overlapping area, and the HCPs are from ICESat global land surface altimetry data (GLAH14). In Section III, we will introduce the extraction and selection method of TP chips and HCPs. 2) TP chips and HCPs are used for elevation adjustment to solve the compensation parameters of DEM data. 3) Compensation for systematic elevation error of InSAR-DEMs. 4) Use the compensated DEM to match the TPs and GCPs.
After steps 1)-3), the systematic elevation error of InSAR-DEM disappears. Therefore, compared with the original InSAR-DEM, the TPs and GCPs matched by the compensated DEM have higher elevation accuracy. In addition, the TPs and GCPs here are obtained using the NCC-CSM matching method proposed in this article, which is introduced in Section III. 5) Use TPs and GCPs for plane adjustment. Fig. 2 is the flow chart of this strategy, where the green and blue boxes represent the elevation and plane adjustment, respectively. These two parts are connected through the elevation error compensation step.

B. Adjustment Model
As shown in Fig. 2, the overall 3-D adjustment process includes plane adjustment and elevation adjustment. Here, the DEM adjustment method based on the function model is used for elevation adjustment, and the method based on RFM is used for plane adjustment. Both use the polynomial model to express systematic errors of InSAR-DEMs.
1) Elevation Adjustment: Since the instrument calibration cannot reduce all errors, systematic height error is still caused by baseline uncertainty, instrument drifts, and other factors in the generated DEM. Compared with the other factors, the parallel baseline error has a more significant impact on the generated DEM, the mechanism of which is shown as where B err is the paralell baseline error, B ⊥ is perpendicular baseline, θ is the incidence angle, and R is the slant range. Jaime Hueso González et al. [6] expressed the systematic elevation error as a polynomial relationship about image coordinates (X, Y ), where c 0 is the height offset, c 1 and c 2 are the tilts in ground distance and azimuth, respectively. When the data length is long, for example, greater than 2000 km, we can consider using higher order polynomials to model the systematic elevation error [6], [12]. The observation equations of TP chips and HCPs can be listed according to (3) where t and p are the indexes of TPs and GCPs, I and J are the indexes of DEM data, h is the height extracted from the generated InSAR-DEM, and h HCP is the height of HCPs. Equation (4) can be further sorted as follows: Then they can be sorted into the form of matrixes where L and H are the vectors of l and h, respectively. A is the design matrix. vector of height adjustment parameters. According to (7), we define the loss function as follows: Then the least-square (LS) solution of (8) is The RFM [23], [24] is a high-precision fitting of the rigorous geolocation model of SAR sensors, which expresses the normalized pixel coordinates of the SAR image as the ratios of cubic polynomials about normalized geographical coordinates, as follows: where (x, y) is the normalized coordinate of image coordinate (X, Y ), and (B, L, H) is normalized coordinate of geographical coordinate (latitude, longitude, height). P 3,i (i = 1, 2, 3, 4) is a cubic polynomial about (B, L, H). The adjustment method based on the RFM models [25] the systematic horizontal error as a polynomial relationship about image coordinates (X, Y ), as follows: where (ΔX, ΔY ) is the parametric error model. a 0 and b 0 are translation parameters. a 1 , b 1 , a 2 , and b 2 are linear parameters. When modeling systematic horizontal errors, we can also use higher order polynomials, such as affine transformation models [25]. In comparison, the linear model (11) is more widely used in adjustment. The adjustment model based on RFM has been widely used in block adjustment of remote sensing imagery [26], [27]. The adjustment model (11) is solved by inputting GCPs and TPs, which are used to correct the absolute and relative systematic plane error, respectively. For TPs and GCPs, we can list two sets of equations about range and azimuth, as follows: where t and p are the indexes of TPs and GCPs, respectively. I and J are the indexes of InSAR data. (X, Y ) is the pixel coordinate of TPs or GCPs obtained by data matching. Then, we can obtain the corresponding error equations. For TPs, For GCPs, Equation (14) and (15) can be sorted into the form of matrixes as follows: where y = (a 0,1 , a 1,1 , a 2,1 , . . . , a 0,N , a 1,N , a 2,N ) is the vector of plane adjustment parameters about pixel coordinate X. According to (16), we define the loss function as Then the LS solution is . (18) Following the process from (12) to (18), we can get the solution of

III. INPUT DATA OF THREE-DIMENSIONAL ADJUSTMENT
In addition to strategy and model, the input data are another critical point for InSAR-DEMs. As shown in Fig. 2, the input data in the adjustment process are TPs chips, HCPs, TPs, and GCPs. In adjustment, TPs or TP chips correspond to the same ground object in different DEMs, and they are used to constrain the relative plane and elevation accuracy of InSAR-DEMs, respectively. GCPs and HCPs are used as external reference data to improve the absolute plane and elevation accuracy. This article proposes the NCC-CSM to extract TPs and GCPs. In this section, we first review the extraction and selection methods of TP chips and HCPs and then introduce the proposed extraction and selection methods of TPs and GCPs.

A. Extraction and Selection of TP Chips
Interferometric DEM is usually contaminated by Gaussian white noise, and elevation points close to the boundary may even be affected by colored noise [12]. To ensure the relative accuracy and reliability of TPs, Huber et al. [14] put forward the concept of TP chips, the core idea of which is to calculate the median elevation in the overlapping area as the observation values to avoid the influence of outliers, and have developed a process for extracting TP chips. First, eliminate invalid pixels, such as shadow, water, and low-coherence areas. Second, set a series of grids called TP chips along the longitude and latitude direction in the overlapping area. There are two typical sizes of the grids, 500 and 1000 m, where the latter are widely used. Finally, assign the median height value to a TP chip after calculating the histogram of valid pixels in each grid. The geographic coordinate of the TP chip is that of the center point in the grid. The elevations of TP chips depend on the elevation histograms in grids. Therefore, when TP chips are located in an area with sizeable topographic relief, the inaccurate horizontal location may cause a significant change in the elevation histogram, thus affecting the elevation accuracy of TP chips.

B. Selection of HCPs
ICESat data provides globally distributed accurate elevation measured points and has been used as an absolute height reference, HCPs, for InSAR-generated DEMs [7], [12]. After adjustment, the vertical positions of DEMs will approach HCPs. However, some survey points will be affected by clouds and vegetation so that the data cannot reflect the actual ground elevation. Therefore, González et al. [15] have put forward extreme and basic selection criteria to select reliable HCPs used in DEM adjustment. First, eliminate the outliers compared with external DEM. Then, screen the data points with one peak of reflected signal and time-bandwidth of fewer than 3.2 ns (extreme selection criteria) or 8.0 ns (basic selection criteria) [15] to ensure a relatively concentrated energy distribution as much as possible. Finally, average all DEM pixels in the footprint of the ICESat point using 2-D Gaussian weight. It should be noted that the elevation datums of ICESat and InSAR-DEM are the geoid and the ellipsoidal surface, respectively. Therefore, it is necessary to use the parameter d _ gdHt the geoid height above the reference ellipsoid, provided by ICESat officials to convert the elevation datum of measurement points.

C. Extraction and Selection of TPs
The TPs are the pixel point corresponding to the same object in the two scene data. Matching TPs usually include three processing steps. First, unify the coordinate system to reduce the influence of distortion on data matching. Second, match data. Third, eliminate the unreliable matching results. This section puts forward innovative methods for the latter two aspects. 1) Matching Method: As mentioned above, the existing DEM matching methods have problems with computation and mismatch. And we hope to apply the NCC method, commonly used in InSAR data processing, to TPs matching. In this section, we first propose the concept of the CSM and then design an NCC matching method based on the CSM (NCC-CSM).
In the spatial domain, the terrain usually does not change suddenly, so the texture features of DEM are far inferior to remote sensing images, such as SAR and optical imagery. Therefore, we consider using a gradient operator, such as Sobel operator, to filter the DEMs before matching, which is usually called the image enhancement method in DIP. After filtering, there will be two gradient maps in the range and azimuth directions, and it is necessary to consider an appropriate scheme to combine them. One strategy is to solve the slope map directly, as follows: where s is the slope of DEMs, h x and h y are gradients in the range and azimuth direction, respectively. However, the slope map obtained by (19) no longer contains gradient information in two directions, which reduces the number of features to be matched. Given this, we expect to use the idea of the complex function to integrate the two gradient maps. In interferometric DEMs, the elevation h can be regarded as a function of the complex variable z = x + iy, where (x, y) is the plane coordinate of DEMs and i = √ −1 is imaginary unit. According to the Cauchy-Riemann equation, the differential of the complex function h to z can be expressed as The above formula has unified the gradient maps h x and h y by the complex function theory. We define h z as a CSM and have proved in Appendix A that the slope and aspect of terrain are unified in (20). After generating the CSM of DEM, we use the NCC method [22] for matching. The core idea of NCC method is to solve the sliding correlation coefficient of two complex images, and the peak position of the normalized correlation coefficient map is their 2-D offset. This method has been widely used in InSAR and DInSAR technologies [28], [29]. Equation (21) shows the principle of solving the NCC map of two complex images. S 1 and S 2 are the CSMs. C and R are cross-correlation and autocorrelation operation, respectively, 2) Selection Criteria: Compared with the NCC-DEM method [22], [30], the normalized correlation coefficient obtained by the CSM is lower, mainly due to the following two reasons. First, calculating gradient is the process of image enhancement, which enlarges the difference between pixels. Second, the CSM contains more terrain information, as the proof in Appendix A, which increases the difference between matched image pairs. Therefore, only using the correlation coefficient to screen GCPs is unsuitable for the proposed NCC-CSM method. Given this issue, we propose a screening criterion based on the peak-side-lobe ratio (PSLR). The meaning of PSLR is the ratio of the maximum values of the primary and secondary peaks of NCC maps, as follows: where ρ p and ρ s are the correlation coefficients of the primary and secondary peaks, respectively. Therefore, a greater PSLR means a more significant correlation between the image pairs corresponding to the peak than other matching results, and thus, the higher matching reliability. The critical step to calculating PSRL is determining the bandwidth of the prominent peak. In this article, we expect to use a 2-D Gaussian function, as follows, to fit the prominent peak and then determine the bandwidth through the parameters of the fitted Gaussian surface Besides the similarity between the Gaussian surface and NCC peak shape, there are two reasons for choosing the Gaussian surface to fit the prominent peak. First, the Gaussian surface gives the bandwidth of peaks, i.e., σ x and σ y . Second, leastsquare method (LSM) can be directly used to solve parameters (k, σ x , σ y , σ xy , x 0 , y 0 ) in (23) after a simple transformation. Appendix B describes the methods of Gaussian surface fitting, after which we can eliminate the pixels within the main peak bandwidth through the logic calculation to calculate the maximum value of other pixels, then get the PSLR using (22).

D. Extraction and Selection of GCPs
GCPs are the reference of the absolute plane accuracy of InSAR-generated DEMs. After adjustment, the horizontal positions of DEMs will approach GCPs. Like the matching method of TPs, this article attempts to match the GCPs from the external DEM using the proposed NCC-CSM method. It should be noted that the resolutions of external public DEMs are usually lower than that of interferometric DEMs. For example, the resolutions of SRTM DEM and AW3D30 DEM currently are about 30 m [17], [31]. Therefore, interferometric and external DEMs cannot be directly coregistered at the same scale. Given this issue, we modify the matching strategy and the calculation method of the CSM and put forward novel selection criteria for GCPs.

1) Calculation of CSM:
The slope map of external DEM is still calculated directly using the gradient operator, while the slope operator of InSAR-DEM needs to be expanded. Fig. 3 shows the reasons and details of the extension process. The red circle represents the external DEM grid points, and the blue circle represents the InSAR-DEM grid points. It can be seen that the resolution of interferometric DEM is twice that of external DEM, which means that the operation between adjacent pixels of external DEM (R 1 and R 2 ) is not the same as the operation between adjacent pixels of InSAR-DEM (B 1 and B 2 ) but the interval pixels (B 1 and B 3 ). Therefore, it is necessary to expand the operator template. Assuming that the size of the original operator is M × N , and the resolution of InSAR-DEM is K times that of the external DEM, then the size of the extended operator is [(M − 1)K + 1] × [(N − 1)K + 1], the number with coordinate ((i − 1)K + 1, (j − 1)K + 1) in the extended operator is equal to the number with coordinate (i, j) in the original operator, and the values of rest of the positions in the extended operator are 0. Fig. 3 shows the extended result of the Sobel operator when K = 2.
2) Matching Strategy: The CSMs of external DEM and InSAR-DEM obtained by the above method are still at their original resolutions, so we can consider a sliding registration scheme to realize their coregistration. Assuming that the resolution of InSAR-DEM is K times that of the external DEM,  N) to form the matching template of InSAR-DEM and to take out the pixels with coordinates (u, v) to form the matching template of external DEM. Second, we can calculate the correlation coefficient ρ(x, y) using (21) and get a set of correlation coefficient graphs by changing (x, y) evenly. Finally, according to the peak position of the correlation coefficient map, we can obtain a matching GCP. Fig. 4 shows this matching strategy. The red dots indicate the elevation pixels in the external DEM. The blue dots with green borders indicate the current matching template of InSAR-DEMs and with yellow borders indicate the matching template of InSAR-DEMs the next time. The arrow indicates the sliding direction of the matching template.
3) Selection Criteria: Due to the low resolution of external DEM, the adjacent matching templates in InSAR-DEM will have high similarity to the corresponding template in external DEM, which means that the prominent peak of the NCC map has a large bandwidth. Therefore, all points in the NCC map may be within the main peak bandwidth. In this case, the PSLR-based filtering strategy proposed in Section III-C will fail. Given this issue, we propose a joint screening criterion based on the main peak bandwidth and the Gaussian surface difference (GSD). The core idea is to screen the matching results with a narrow bandwidth and a slight GSD. A small bandwidth (σ x , σ y ) means that the image pairs corresponding to the maximum NCC have a more significant relevance than other pairs. The GSD here refers to the root-mean-square (rms) difference between the NCC map and the fitting surface. The expression is as follows: where (x 0 , y 0 ) is the peak position of the fitted Gaussian surface in (23). ρ g u,v and ρ u,v represent the fitted and original correlation coefficients, respectively. r x and r y limit the solution range of GSD. Therefore, a slight GSD means that all pixels near the prominent peak stably approach (i.e., with slight fluctuation) the NCC maximum value. It should be noted that when selecting TPs and GCPs, areas with severe geometric distortions, especially those with many shadows and extreme foreshortening distortions, need to be avoided to ensure matching accuracy. Because these geometric distortions will lead to errors between the actual and the inverted elevations of ground objects, thus increasing the difference between the matched DEM pairs.

A. Experiment Data
This experiment selects 29 CoSSC data of TanDEM-X covering Henan Province, China. Twenty scenes are ascending-orbit, and nine are descending-orbit. The data are acquired from 2013 to 2019, and most of them are measured in winter. Table I shows the information about the experimental CoSSC data, including the acquisition time, the number of images contained in each track, the resolution in azimuth and range direction, and the ascending (A) or descending (D) orbit. The GCPs used in the experiment are selected from GLAH14 data of ICESat according to the criterion in [15]. The height checkpoints are from land and vegetation height data (ATL08) of ICESat-2, the vertical uncertainty of which is 0.2 m for flat terrain and 2.0 m for mountainous terrain [32]. The external DEM selected in the experiment is ALOS World 3D-30 m (AW3D30), with a resolution of about one arcsec and a nominal vertical accuracy of about 5.0 m [31]. Fig. 5 shows the location information of experimental data and the DEM in the experimental area, respectively, where all the CoSSC data are numbered according to their acquisition time for analysis. The experimental region contains various geomorphic features, including plains, mountains, and rivers, with a maximum height difference of 2200 m.
The experiment is carried out according to our proposed method in Section II. We focus on the steps of TPs extraction, GCPs extraction, and 3-D adjustment.

B. Extraction of TPs
This experiment extracts all TPs using the method proposed in Section III-C. This section focuses on the performance of resampling polynomials, matching methods, and screening criteria in practice.  Fig. 6 shows four groups of matching results, i.e., NCC maps, of TPs. The first column shows the terrain to be matched, and the second to fourth columns are the NCC maps of CSM, DEM, and slope map, respectively. The first to fourth rows are the matching results of the same orbit data, mountainous and flat areas of different tracks, and ascending/descending data. Table III shows the evaluation indicators of the matching results shown in Fig. 6, including PSLR and peak bandwidths.

1) Matching Results:
From Fig. 6, there are still convex peaks in Fig. 6(c), (g), and (o), but it disappears in Fig. 6(k), which means that the NCC-DEM method can still get a specific matching effect in mountainous areas, but it will fail in flat areas. It is because the elevation of mountainous areas fluctuates wildly, e.g., Fig. 6(e) contains more valley areas, which provides more texture features for matching image pairs. However, flat areas have small fluctuations, e.g., the maximum elevation difference in Fig. 6(i) is only about 30 m, and the texture segments are not apparent. On the other hand, most cities are located in plain areas. Human activities will change the ground elevation and increase the difference between matching templates, failing the NCC-DEM method. In contrast, the NCC-CSM method can get good matching results in these areas, and Fig. 6(f), (j), and (n) have an apparent prominent peak. Because calculating CSMs is a process of image enhancement (i.e., high-pass filtering of DEMs), the NCC-CSM method actually matches the changing trend of DEM in the spatial domain. Similarly, the NCC-SM method can also obtain a specific matching effect, as shown in Fig. 6(j). Nevertheless, the prominent peaks of NCC-SM are not as apparent as that of NCC-CSM because slope maps only include slope information but ignore aspects. The same slope may correspond to different aspects. For example, both vectors (h x , h y ) and (−h x , h y ) correspond to the slope s = tan −1 h 2 x + h 2 y . Therefore, using slope alone is not enough to describe elevation points. As proved in Appendix A, CSM contains slope and aspect information, which increases the difference between matching image pairs. Therefore, Fig. 6(j) and (n) have apparent prominent peaks. In terms of digital indicators, as shown in Table III, the coefficient map bandwidths of NCC-CSM are the smallest, usually less than 10 pixels, while that of NCC-DEM are the largest, in the order of 10 2 to 10 3 . At the same time, the NCC-CSM method has the highest PSLR. Therefore, on the whole, the proposed matching method is better than the NCC-DEM and NCC-SM methods, especially in flat areas.
The size of image pairs used to match TPs is 600 × 600 pixels, and the platform is MATLAB R2020a on Intel(R) Core(TM) i5-10300H CPU at 2.50 GHz. We list the time consumption of TPs in Table II, where the CSM, SM, and NCC represent the calculating time of the CSM, SM, and NCC map, respectively. It can be seen from Table II that the calculation time of CSMs is slightly longer than that of SMs, and the calculation time of the correlation coefficient map of the NCC-CSM method is longer than that of NCC-DEM and NCC-SM methods. It is because the CSM is a complex image containing gradient graphs in two orthogonal directions, and the algorithm of the complex number is more complex than the real number.
Moreover, both Fig. 6(i) and (m) are located in flat areas, but the NCC map Fig. 6(k) contains no peak. These two images are located in flat areas with little topographic relief and contain some urban or rural buildings. The significant difference between the two images is that the elevation structure of Fig. 6(i) is relatively complex, including multiple roads and buildings. Because of the similarity of building layouts in cities or villages, Fig. 6(i) has low identification in the preset matching area, thus leading to no peak in its correlation coefficient map when using the NCC-DEM method. In contrast, the structure in Fig. 6(m) is relatively simple, with only one road and a few buildings. Its recognition degree is relatively higher in the preset matching area, and when matching with the NCC-DEM method, its correlation coefficient map Fig. 6(o) shows a peak.
2) Selection Criteria: In the matching experiment of TPs, the proposed screening criteria ensure the reliability of the matching results. We have set the screening threshold of PSLR as 1.5 and rejected the results less than the threshold. Fig. 7 shows some matching results, in which Fig. 7(g) to (l) are screening results, Fig. 7(a) to (f) are culling results, and fourth column pictures are corresponding fitted Gaussian surfaces. Table IV shows the PSLR and GSD of the matching results. It should be noted that the GSD listed here is used to illustrate the fitting accuracy of Gaussian surface, not as a screening criterion. It can be seen that the Gaussian surface can accurately fit the prominent peak of the correlation map, and the GSDs of fitting results are in the order of 10 −2 , which guarantees the calculation of PSLR. Due to the unique [e.g., Fig. 7(d)] or inconspicuous [ Fig. 7(a)] geomorphic features of the matching image pairs, Fig. 7(e) and (b) do not contain an apparent prominent peak, thus leading to low PSLRs, as shown in Table IV. The threshold of PSLR can effectively eliminate these results.

C. Matching of GCPs From External DEMs
This experiment extracts all GCPs using the method proposed in Section III-D. This section focuses on the performance of the proposed matching method and screening criteria in practice.
1) Matching Results: Fig. 8 shows some results of GCPs matching. The first column shows the elevation of matching areas, and the second to third columns are the normalized  III  INDICATORS OF TPS MATCHING RESULTS   TABLE IV  NUMERICAL INDICATORS Table VI shows the numerical indicators of these matching results. Considering that PSLR calculation will be invalid due to the large bandwidth, we have listed the maximum relative difference (MRD) of NCC maps in Table VI to measure the prominence of prominent peaks. MRD is calculated as follows: where ρ max and ρ min are the maximum and minimum in the NCC map, respectively. From Table VI, the MRDs of Fig. 8(b) and (e) are 0.612 and 0.981, while that of Fig. 8(c) and (f) are 0.002 and 0.001, and the peak bandwidths of Fig. 8(b) and (e) are far narrow than Fig. 8(c) and (f). Therefore, the proposed method can  obtain prominent NCC peaks with the largest MRD and least bandwidth, which is consistent with the matching results of TPs. Comparing Fig. 8(b) with Fig. 8(e), we can see that the correlation coefficient in flat areas is lower than in mountainous areas. There are many reasons behind this phenomenon. On the one hand, the gradient characteristics of flat areas are less significant. On the other hand, the elevation change of plain areas is usually more significant than that of mountainous areas due to artificial activities. Compared with the matching results of TPs, the correlation coefficient maps of GCPs have larger prominent peak bandwidths. Table VI shows that the bandwidths of NCC-CSM are usually more than 10 pixels and sometimes more than 20 pixels, which has exceeded the preset maximum offset of 20 pixels. It is because the resolution of external DEM AW3D30 is about 25.0 m (six to eight times that of InSAR-DEM), which leads to the high similarity between adjacent sliding matching templates of InSAR-DEM. Table V shows the time consumption of GCPs, and the size of image pairs used to match GCPs is 200 × 200 pixels. The time consumption of GCPs matching is similar to that of TPs; that is, the time consumption of the NCC-CSM method is longer than that of the NCC-DEM and NCC-SM methods. The reason can still be attributed to the calculation of complex and real numbers.
In addition, whether matching TPs or GCPs, the bandwidth of the NCC-CSM methods is far less than that of the NCC-DEM method, which means that the matching image pairs of NCC-CSM are more similar than that of NCC-DEM. As mentioned above, calculating the CSM is an enhancement process of DEMs. Furthermore, the CSM contains the slope and aspect information of DEM, which increases the difference of adjacent matching templates, resulting in the significant similarity of the matching pairs corresponding to the primary peak in an NCC map.
2) Selection Criteria: The selection strategy for calculating PSLR fails because all the pixels in a correlation coefficient map are likely within the main peak. In this case, the proposed selection strategy based on GSD can screen out reliable GCPs. In the experiment, we set the threshold value of GSD to 0.02, and the matching results that are greater than the threshold value will be rejected. Fig. 9 shows the selected and culled matching results of GCPs, respectively. Table VII shows the corresponding numerical indicators. It can be seen that Fig. 9(b) and (e) have large fluctuations in the peak area, and their corresponding GSDs are 0.0257 and 0.0650, respectively. However, the variation trends in the peak area of Fig. 9(h) and (k) are stable, with 0.0056 and 0.0132 GSDs, which are significantly smaller than that of Fig. 9(b) and (e), thus can be selected as GCPs used in the following adjustment process.

D. Adjustment of Plane and Elevation
We use the method proposed in Section II for 3-D adjustment. First, we take the TP chips and elevation GCPs as input data to adjust elevation using the method based on the function model. After compensating for elevation error, we take the matching TPs and GCPs as inputs for plane adjustment using the method based on RFM.
1) Elevation Adjustment: We use the method based on the function model for elevation adjustment. HCPs are selected from ICESat GLAH14 data according to extreme selection criteria [15]. The 16th and 17th data do not contain any HCPs, called uncontrolled data, and the rest are called controlled data. TP chips are extracted using the method described in Section III, and elevation checkpoints are from ICESat-2 ATL08 data.
First, we show the mosaic rendered map of InSAR-DEMs before adjustment, as Fig. 10. Due to inconsistent systematic elevation error in two adjacent scenes, the mosaic map has an apparent elevation step phenomenon. In order to vividly show this phenomenon, we show the 3-D elevation views of the marked areas in Fig. 10(b) and (c). The maximum and minimum heights of Fig. 10(b) are about 73 and 66 m, respectively, and that of Fig. 10(c) are about 60 and 48 m, respectively. Subsequently, we carried out height adjustment. We use the root-mean-square   error (RMSE) to measure the elevation accuracy before and after adjustment, as where h DEM and h ck are the elevations of InSAR-DEM and checkpoints, respectively, N ck is the number of checkpoints. Table VIII shows the elevation accuracies of the controlled, uncontrolled, and whole areas and the relative accuracy of InSAR-DEM, i.e., the rms elevation differences of TP chips, before and after adjustment. Fig. 12 shows the elevation accuracy of each track before and after adjustment. Overall, the elevation accuracy of the controlled area and the uncontrolled areas are 1.70 and 1.80 m before adjustment. After adjustment, they are improved to 1.36 and 1.25 m, and the improvements are 0.34 and 0.55 m, respectively. The overall adjustment accuracy is improved by about 0.35 m, and the elevation accuracies of each track data have been improved to varying degrees, as shown in Fig. 12. Diverse factors affect the adjustment effect, such as the original accuracies of DEMs, the terrain factors, and the number of control points. First, its accuracy improvement space is narrow for DEM with high accuracy, and the adjustment effect may be insignificant. Second, some areas inevitably lack control points due to the uneven distribution of ICESat data, and their accuracy improvement may be lower than the area with sufficient control points. In addition, the vertical accuracies of ICESat and DEM data in mountainous areas are generally lower than in flat areas. So even after adjustment, the vertical errors of the first to third tracks are still relatively high. After adjustment, the difference of TP chips decreases from 1.67 to 0.48 m. Fig. 11 shows 3-D view of the marked position in Fig. 10(a) after elevation adjustment. It can be seen that the elevation step effect disappears. The above experimental results prove that our proposed adjustment method has an excellent effect on improving the relative and absolute elevation accuracy of InSAR-DEMs.
2) Plane Adjustment: As mentioned in the previous part, InSAR-DEM inevitably has systematic plane errors due to the errors of geolocation parameters, which cannot be eliminated only through elevation adjustment. We selected a specific area in interferometric DEMs to illustrate the impact of plain error, as shown in Fig. 13, where Fig. 13(a) is the DEM mosaic map after elevation adjustment, and Fig. 13(b) is the corresponding optical image. The west part in Fig. 13(a) belongs to 12th data, and the east part belongs 9th data. It can be seen that there is a railway passing through this area, so the elevation map of the corresponding position shows a straight line from southeast to northwest. However, because the two scene data have different degrees of plane systematic error, the straight line has an obvious fracture near the mosaic line, and the three red notes mark the fracture area. Therefore, it is necessary to correct the systematic plane error of InSAR-DEM data.
In the plane adjustment, we select the points matched from external DEMs with small GSDs in each scene data as the check points (CKPs) to measure the absolute plane accuracy, as follows: where P CKP,i and P DEM,i are the plane coordinates of CKPs and corresponding points in InSAR-DEMs, respectively. And we use the points matched between InSAR-DEMs with high PSLR to measure the relative plane accuracy. Before adjustment, the relative and absolute plane accuracies are 14.70 and 12.22 m, respectively. After adjustment, they are improved to 3.17 and 6.44 m, respectively. Fig. 14 shows the systematic plane error of each track. Fig. 15 shows the error vector map of CKPs, in which the red circles are the real coordinates  of CKPs, and the arrows represent the difference between the geolocation and the real coordinates of CKPs. It can be seen that before the adjustment, each scene image had different degrees of plane systematic error. For example, the plane error vector in the 13th DEM mainly deviates to the northeast along with the range direction, likely caused by the range error. After adjustment, the systematicness of CKP error disappears, as shown in Fig. 15(b), and the plane accuracy of each track data has been significantly improved. Fig. 16 shows the DEM mosaic rendered map after plane adjustment. We separately list the mosaic map of the  railway area in Fig. 13, as shown in Fig. 17(b). We also list the mosaic map before the plane adjustment in Fig. 17(a) for the convenience of comparison. It can be seen that the terrain of this area has roughly shown as a continuous straight line and the fracture phenomenon in this area has disappeared. In general, the plane adjustment effectively improves the absolute and relative plane accuracy of DEM.

V. CONCLUSION
This article proposes a novel 3-D block adjustment method for spaceborne bistatic InSAR system based on general models to correct the horizontal and vertical systematic errors of InSAR-DEMs. The core idea of the proposed method is first to correct systematic elevation errors using the DEM adjustment method based on the function model and then to correct systematic plane errors using the method based on the RFM model after elevation correction. In addition, we put forward the NCC-CSM matching method to extract TPs and GCPs used in 3-D adjustment. We used 29 TanDEM-X CoSSC data covering Henan Province, China, in experiments and focused on three parts: 1) TPs matching and screening, 2) GCPs matching and screening, and 3) 3-D adjustment. The matching results showed that the proposed NCC-CSM method could effectively match TPs and GCPs, especially in flat areas, and the screening criteria can accurately eliminate invalid matching results. The adjustment experiment showed the horizontal and vertical geometric discontinuities of InSAR-DEM and proved that the proposed 3-D adjustment strategy could effectively improve the 3-D accuracy of InSAR-DEM and make up for the deficiency of the DEM adjustment method based on the function model. After adjustment, the elevation accuracy of InSAR-DEM was improved from 1.70 to 1.35 m, and the elevation step effect of the overlapping area was significantly ameliorated. The absolute and relative plane accuracies were improved by 5.58 and 11.53 m, respectively, and the plane geometric discontinuity in the overlapping area disappeared. In the future, we can utilize more CoSSC data from other regions, including different landforms and countries, especially regions with high-precision 3-D control points, to further confirm our method.

APPENDIX A PROOF OF THE PHYSICAL MEANING OF COMPLEX GRADIENT
Assuming that the elevation surface h = h(x, y) is a real function of 2-D image coordinates (x, y), thus the outer normal vector of the point (x, y, h) can be expressed as n = (h x , h y , −1), where h x = ∂h ∂x and h y = ∂h ∂y . According to [33], the terrain slope is the change rate of the vertical elevation with the distance from the water plane, and the aspect refers to the projection direction of the external normal vector on the horizontal plane. In Fig. 18, n is the external normal vector of point p, and s and a represent the slope and aspect, respectively. Therefore, the slope and aspect of the point can be expressed as Comparing (29) with (20), it can be seen that the norm of h z is actually the tangent value of the slope, and the phase of h z is aspect, as follows: Furthermore, (20) can be expressed as where s t is the tangent value of the slope.

APPENDIX B FITTING OF GAUSSIAN SURFACE
Since the Gaussian surface function is nonlinear, using (23) will increase the computational complexity. Therefore, we perform a logarithmic operation on (23) to obtain In the above equation, some parameters are coupled together (e.g., σ x and x 0 ), so it is challenging to use the LSM to solve directly. Therefore, we simplify (31) as follows: Then, the six unknown parameters of the Gaussian surface [k, x 0 , y 0 , σ x , σ y , σ xy ] T are converted to q = [q 0 , q 1 , q 2 , q 3 , q 4 , q 5 ] T by (32) and (33). Let a = [1, x, y, x 2 , y 2 , xy] T , then (31) can be written as log(g) = a T q.
Then, we select a series of points (x i , y i ) near the NCC peak as the observation data, and the error equation can be listed as We can use the LSM to solve q. Finally, we can determine the range and azimuth bandwidth of the prominent peak through σ x and σ y . For simplicity, we can directly make them equal to σ x and σ y , i.e.,