A Three-Dimensional Block Adjustment Method for Spaceborne InSAR Based on the Range-Doppler-Phase Model

Abstract

The block adjustment method can correct systematic errors in the bistatic Synthetic Aperture Radar Interferometry (InSAR) satellite system and effectively improve the accuracy of the InSAR-generated Digital Elevation Model (DEM). Presently, non-parametric methods, which use the polynomial to model the systematic errors of InSAR-generated DEMs, are most frequently used in spaceborne InSAR-DEM adjustment. However, non-parametric methods are not directly related to the physical parameters in the InSAR imaging process. Given the issue, this paper conducts adjustments in the parameter domain and proposes a three-dimensional block adjustment method for spaceborne bistatic InSAR systems based on the Range-Doppler-Phase (RDP) model. First, we theoretically analyze the sensitivities of spatial baseline, azimuth time, and slant range to the RDP geolocation model and confirm the analysis method with a simulated geolocation result. Second, we use total differential and differential geometry theories to derive adjustment equations of available control data based on sensitivity analysis. Third, we put forward an iterative solution strategy to solve the corrections of parallel baseline, azimuth time, and slant range to improve the plane and elevation accuracies of InSAR-generated DEMs. We used 29 scenes of TanDEM-X Co-registered Single look Slant range Complex (CoSSC) data to conduct simulated and real data experiments. The simulated results show that the proposed method can improve the accuracies of baseline, range, and timing to 0.05 mm, 0.1 m, and 0.006 ms, respectively. In the real data experiment, the proposed method improves the plane and elevation accuracies to 4.14 m and 1.34 m, respectively, and effectively suppresses the fracture phenomenon in the DEM mosaic area.

1. Introduction

The bistatic Synthetic Aperture Radar Interferometry (InSAR) technology [1,2] uses two SAR satellites in formation to observe the ground and obtains the phase difference through interference processing to inverse the Digital Elevation Model (DEM). The TerraSAR-X/TanDEM-X [3,4] developed by German Aerospace is a typical representative of spaceborne bistatic InSAR systems.

The interference factors, such as measurement noise, can affect the accuracy of InSAR geometric parameters, thus leading to systematic errors in inversed DEMs [2]. The InSAR block adjustment technology [5,6] can jointly correct these errors. The core idea is to list the loss function using control data and solve adjustment parameters using optimization algorithms. According to the adopted models, the current adjustment methods can be divided into two categories: the method based on general models [5,7], called the non-parametric method, and the method based on strict models [6,8,9], called the parametric method. The former models the systematic elevation errors of DEMs as polynomials about image coordinates and the systematic errors of each DEM can be directly compensated by polynomial coefficients. The latter method is developed based on the physical geolocation model, the Range-Doppler-Phase (RDP) model [10], of the bistatic InSAR system to solve the measurement errors of slant range, baseline, and other parameters. After adjustment, subscribers need to eliminate the measurement errors from the original parameters to re-inverse DEMs. The non-parametric method is widely employed in spaceborne InSAR systems [11]. In 2010, González, Jaime Hueso et al. verified the general adjustment method using simulated experiments [12]. Huber, Martin et al. put forward the concept of tie point chips to constrain the height accuracy of the DEM overlapping area [13]. In 2012, Gruber, Astrid et al. obtained the first result of TanDEM-X DEM using the general method, proving that the elevation accuracy after adjustment can reach 2.0 m. At present, this method has been applied in TanDEM-X Mosaic and Calibration Processor (MCP) [14]. The parametric method is usually used to adjust airborne InSAR systems [15,16]. In 2011, Jin Guowang et al. jointly calibrated the baseline error of airborne systems using the block adjustment method [16]. In 2012, Ma Jing et al. proposed an airborne three-dimensional block adjustment method by extending the Range-Doppler model [8]. After adjustment, the plane and elevation accuracies are 1.60 m and 2.22 m, respectively.

However, the current adjustment methods still have room for improvement. On the one hand, the non-parametric method can not be applied to the plane adjustment because it only models the elevation error [17]. On the other hand, the adjustment parameters (i.e., polynomial coefficients) in the non-parametric method do not have strict physical significance and cannot correspond to the physical parameters in the InSAR inversion process. Most parametric methods are derived from airborne geometry without considering the characteristics of spaceborne InSAR systems. For example, the airborne adjustment usually uses manually-deployed corner reflectors as external control points, which are strong reflection points in SAR images with high recognition and incredible three-dimensional accuracy [18,19,20]. However, due to the limitation of cost, accessibility, and other issues, the feasibility of this strategy is extremely low for spaceborne large-area mapping. It undoubtedly brings challenges to spaceborne InSAR adjustment.

In response to the above problems, this paper proposes a three-dimensional spaceborne InSAR-DEM block adjustment method based on the RDP geolocation model. In Section 2, we analyze the sensitivities of physical parameters, derive the mathematical expressions between geolocation errors and the uncertainty of physical parameters by the total differential and differential geometry theories, and determine the slant range, spatial baseline, and timing parameters as adjustment parameters. In Section 3, we derive the adjustment equations of available control data, including control points and tie points, establish the spaceborne InSAR-DEM adjustment model based on the RDP model, and propose an iterative solution strategy. In Section 4, we use 29 scenes of CoSSC data to design the simulated and real data experiments to verify the proposed adjustment method. The simulated results show that the proposed method can accurately correct the slant range, parallel baseline, and timing parameters, and the corresponding uncertainties are 0.1 m, 0.05 mm, and 0.006 ms, respectively. In the data experiment, the proposed method improves the plane and elevation accuracies of InSAR-generated DEM to 4.14 m and 1.34 m, respectively, where the height adjustment accuracy is consistent with the popular non-parametric adjustment method.

2. Analysis of the Range-Doppler-Phase Model

From the solution perspective, sensitivity analysis and parameter adjustment are opposite processes. The former is to calculate the geolocation errors with known parameters’ errors, while the latter is to solve the corrections of InSAR parameters with known geolocation errors. In this section, we use total differential and differential geometry theories to analyze the sensitivities of InSAR geometric parameters, then reveal the mathematical relationship between the parametric and geolocation errors, and verify the correctness of the analysis method.

2.1. Range-Doppler-Phase Model

The rigorous geolocation model of the spaceborne bistatic InSAR system includes equations of slant range, Doppler frequency, and interferometric phase, called the RDP model [10], as Equation (1).

$$\left\{\begin{array}{cc}\hfill & |\overrightarrow{P_m}-\overrightarrow{T}{|}^2=R^2\hfill \\ \hfill & (\overrightarrow{P_m}-\overrightarrow{T})·\overrightarrow{V_m}=-\frac{1}{2}\lambda f_{dc}R\hfill \\ \hfill & |\overrightarrow{P_m}+\overrightarrow{B}-\overrightarrow{T}{|}^2={(R+\Delta R)}^2\hfill \end{array}\right.$$

(1)

where $\lambda$ is the electromagnetic wave length. R and $f_{dc}$ are the slant distance and Doppler frequency of the master satellite. $\overrightarrow{P_m}=(P_{m,x},P_{m,y},P_{m,z})$ and $\overrightarrow{V_m}=(V_x,V_y,V_z)$ are the position and velocity of the master satellite, respectively, and $\overrightarrow{B}=(B_x,B_y,B_z)$ is the spatial baseline, representing the spacial vector formed by Antenna Phase Centers (APCs) of the master and slave satellites, i.e., $B=\overrightarrow{P_s}-\overrightarrow{P_m}$, where $\overrightarrow{P_s}$ is the position of the slave satellite. $\overrightarrow{P_m}$, $\overrightarrow{V_m}$, and $\overrightarrow{B}$ can be expressed as polynomial functions of azimuth time t. $\Delta R$ represents the wave path difference between the master and slave satellites, which can be calculated using the absolute interference phase, i.e., $\Delta R=\frac{\lambda {\varphi }_{abs}}{2\pi }$. $\overrightarrow{T}=(X,Y,Z)$ is the geographical coordinates of the ground objects and satisfies the Earth’s ellipsoidal Equation [21], as Equation (2).

$$\frac{X^2+Y^2}{{(R_a+h)}^2}+\frac{Z^2}{{(R_b+h)}^2}=1$$

(2)

where $R_a$ and $R_b$ are the semi-major and semi-minor axises of the ellipsoidal model, respectively, and h is the elevation of the ground target. Therefore, after calculating $\overrightarrow{T}$ through Equation (1), the target’s elevation h can be solved by Equation (2). It should be noted that Equations (1) and (2) are based on the Earth-Centered Earth-Fixed (ECEF) coordinate system.

2.2. Analysis for RDP Model

In positioning model Equation (1), the uncertainty of parameters such as spatial baseline will affect the three-dimensional positioning results. In this section, we use the total differential method and differential geometry theory to analyze the influence of these parameter errors on geolocation results, including horizontal and vertical accuracies. In order to visually illustrate their impacts, we choose a specific set of geometric parameters to draw the sensitivity curves, as

Table 1. The parameters used for sensitivity analysis.

Parameters	Values
Slant range (km)	682.70
Satellite position (km)	(−2706.04, 5088.03, 3765.30)
Satellite velocity (m/s)	(−524.50, 4389.90, −6288.98)
Baseline vector (m)	(0.62, −266.42, 430.36)
Absolute interferometric phase (rad)	−174.18
Target coordinate (km)	(−2070.03, 4915.39, 3486.46)
Eastern unit vector (m)	(−0.91, −0.42, 0.0)
Northern unit vector (m)	(0.23, −0.51, 0.83)
Upper unit vector (m)	(−0.35, 0.75, 0.56)

. The ENU coordinate system is necessary when determining the vertical and horizontal components of three-dimensional geolocation error. Appendix A deduces the conversion relationship between the ENU and ECEF coordinate systems. Table 1 also shows the northern, eastern, and upper unit direction vectors in the ECEF coordinate system.

2.2.1. Spatial Baseline

The spatial baseline is an essential parameter in the inversion process, representing the relative position of the APCs between the master and slave satellites. Figure 1 shows the vector decomposition of spatial baseline in a bistatic InSAR system.

It is generally believed that the parallel baseline error will affect the vertical accuracy of DEMs [22,23,24,25]. In fact, due to the side-view imaging pattern, the elevation error will affect the plane accuracy, which means that the baseline error will affect the three-dimensional accuracy of DEMs. The differential of Equation (1) for $\overrightarrow{B}$ and $\overrightarrow{T}$ is

$$\left\{\begin{array}{cc}\hfill & (\overrightarrow{T}-\overrightarrow{P})·d\overrightarrow{T}=0\hfill \\ \hfill & \overrightarrow{V_m}·d\overrightarrow{T}=0\hfill \\ \hfill & (\overrightarrow{T}-\overrightarrow{P_s})·d\overrightarrow{T}=(\overrightarrow{T}-\overrightarrow{P_s})·d\overrightarrow{B}\hfill \end{array}\right.$$

(3)

where $d\overrightarrow{B}={[d{B}_x,d{B}_y,d{B}_z]}^{\mathrm{T}}$ is the error of baseline $\overrightarrow{B}$. The symbol · stands for inner product operation of vectors. Equation (3) establishes the relationship between the three-dimensional baseline error $d\overrightarrow{B}$ and the resulting geolocation error $d\overrightarrow{T}$, and we can further solve this relationship through matrix transformation. Let

$$\left\{\begin{array}{cc}\hfill \overline{\overline{D}}& ={[\overrightarrow{T}-\overrightarrow{P_m},\overrightarrow{V_m},\overrightarrow{T}-\overrightarrow{P_s}]}^{\mathrm{T}}=\left[\begin{array}{ccc}X-P_{m,x}& Y-P_{m,y}& Z-P_{m,z}\\ V_x& V_y& V_z\\ X-P_{s,x}& Y-P_{s,y}& Z-P_{s,z}\end{array}\right]\hfill \\ \hfill \overrightarrow{u_B}& ={[0,0,(\overrightarrow{T}-\overrightarrow{P_s})·d\overrightarrow{B}]}^{\mathrm{T}}\hfill \end{array}\right.$$

(4)

It can be found from Equation (4) that $\overline{\overline{D}}$ is composed of two line-of-sight (LOS) vectors and a velocity vector. According to the geometric relationship of the bistatic InSAR system [1,2], the three vectors are linearly uncorrelated, which means that $\overline{\overline{D}}$ is a full-rank matrix. Therefore, the geolocation error can be expressed as Equation (5),

$$d\overrightarrow{T}={\overline{\overline{D}}}^{-1}\overrightarrow{u_B}$$

(5)

Then, we can calculate the normalized direction vectors in point $(X,Y,Z)$ (we have shown this process in Appendix A) and solve the components of geolocation error in elevation $d{T}_u$, east $d{T}_e$, and north $d{T}_n$, as follows.

$$\left\{\begin{array}{cc}\hfill d{T}_u& =d\overrightarrow{T}·\overrightarrow{n_u}\hfill \\ \hfill d{T}_e& =d\overrightarrow{T}·\overrightarrow{n_e}\hfill \\ \hfill d{T}_n& =d\overrightarrow{T}·\overrightarrow{n_n}\hfill \end{array}\right.$$

(6)

According to the currently published literature, the baseline uncertainty of the bistatic InSAR system is usually less than 2 mm [7,23,24,26], and that of some round-pass systems is on the order of centimeters [27]. In order to vividly show the impact of spatial baseline error, we have changed the errors of along-track, parallel and perpendicular baselines, respectively, and calculated the resulting three-dimensional geolocation errors according to Equations (5) and (6). At the same time, we calculated the actual geolocation errors according to the RDP model Equation (1) to verify the proposed analysis method. We use dotted and solid lines to express the real and analyzed positioning errors, respectively, in Figure 2. It can be seen that our analysis results coincide with the actual geolocation errors, which demonstrates the effectiveness of the proposed analysis method. The parallel baseline error has a significant impact on the geolocation results. In contrast, the impact of the perpendicular and the along-track baseline errors can be ignored.

2.2.2. Slant Range

The measurement accuracy of the slant range is affected by various factors, including systematic delay, atmospheric delay, and solid earth tide. Currently, the slant range accuracy can generally reach meters. After accurate corrections, it can reach the level of decimeters or centimeters [28,29]. We can follow Equations (3)–(6) to analyze the sensitivity of slant range. Here we only list the differential form of Equation (1) for range error $dR$ and the resulting geolocation error $d\overrightarrow{T}$.

$$\left\{\begin{array}{cc}\hfill & (\overrightarrow{T}-\overrightarrow{P_m})·d\overrightarrow{T}=RdR\hfill \\ \hfill & \overrightarrow{V_m}·d\overrightarrow{T}=-\frac{1}{2}\lambda f_{dc}dR\hfill \\ \hfill & (\overrightarrow{T}-\overrightarrow{P_s})·d\overrightarrow{T}=(R+\Delta R)dR\hfill \end{array}\right.$$

(7)

We set the range error as −15 to 15 m and calculate the resulting three-dimensional geolocation error according to Equation (7), as shown in Figure 3. It can be seen that the geolocation errors in the elevation and east directions are the same magnitude as $dR$.

2.2.3. Orbit and Timing Parameter

The orbit information of the master satellite is another essential parameter in Equation (1). At present, the position and velocity vectors of SAR satellites are highly accurate. For example, the position accuracy of TerrSAR-X can reach 3.0 cm [30], and the velocity and position accuracies of GF-3 orbit data can reach 5 cm and 0.05 mm/s [31], respectively. These errors have little impact on geolocation results, so they can be ignored in adjustment.

The satellite state vectors are usually fitted as polynomials of azimuth time. Therefore, the timing delay will cause inaccurate position and velocity, leading to the plane offset of geolocation results. For example, the timing delay of TerrSAR-X is about 0.18 ms, which will cause a plane error of about 1.27 m [30]. Here we list the differential form of Equation (1) for timing error $dt$ and the resulting geolocation error $d\overrightarrow{T}$.

$$\left\{\begin{array}{cc}\hfill & (\overrightarrow{T}-\overrightarrow{P_m})·d\overrightarrow{T}=(\overrightarrow{T}-\overrightarrow{P_m})·\frac{d\overrightarrow{P_m}}{dt}dt\hfill \\ \hfill & \overrightarrow{V_m}·d\overrightarrow{T}=(\overrightarrow{P_m}-\overrightarrow{T})·\frac{d\overrightarrow{V_m}}{dt}dt+\overrightarrow{V_m}·\frac{d\overrightarrow{P_m}}{dt}dt\hfill \\ \hfill & (\overrightarrow{T}-\overrightarrow{P_s})·d\overrightarrow{T}=(\overrightarrow{T}-\overrightarrow{P_s})·\frac{d\overrightarrow{P_s}}{dt}dt=(\overrightarrow{T}-\overrightarrow{P_s})·\frac{d(\overrightarrow{P_m}+\overrightarrow{B})}{dt}dt\hfill \end{array}\right.$$

(8)

where the symbol $\frac{d·}{dt}$ represents the derivatives for the azimuth time t. All the derivatives can be easily calculated because of the polynomial relationships. Take $\overrightarrow{P_m}$ as an example,

$$\left\{\begin{array}{cc}\hfill & \overrightarrow{P_m}=[P_{m,x},P_{m,y},P_{m,z}]=\left[\sum _{n=0}^{N_p}{p}_{m,x}^n{t}^n,\sum _{n=0}^{N_p}{p}_{m,y}^n{t}^n,\sum _{n=0}^{N_p}{p}_{m,z}^n{t}^n\right]\hfill \\ \hfill & \frac{\overrightarrow{P_m}}{dt}=[\frac{d{P}_{m,x}}{dt},\frac{d{P}_{m,y}}{dt},\frac{d{P}_{m,z}}{dt}]=\left[\sum _{n=1}^{N_p}n{p}_{m,x}^n{t}^{n-1},\sum _{n=1}^{N_p}n{p}_{m,y}^n{t}^{n-1},\sum _{n=1}^{N_p}n{p}_{m,z}^n{t}^{n-1}\right]\hfill \end{array}\right.$$

(9)

where N is the order of polynomials, and $p_{m,x}^n$, $p_{m,y}^n$, and $p_{m,z}^n$ are the fitting coefficients.

Figure 4 shows the sensitivity curve of the geolocation error versus $dt$. It can be seen that $dt$ mainly affects the positioning error in the north (approximately azimuth) direction but has little impact in the east and elevation.

3. Method

We can draw two conclusions from the discussion in Section 2. First, the proposed analysis method combining differential geometry and total differential theory is correct. The sensitivity curves demonstrate that our theoretical analysis method can accurately describe the impact of error sources on InSAR three-dimensional geolocation results. Second, the errors of parallel baseline, slant range, and azimuth time are the main physical parameters that affect the positioning results of spaceborne InSAR. The influences of parallel baseline and slant range are concentrated in height and east–west directions, while that of azimuth time is concentrated in the north–south direction.

According to the above two conclusions, we can apply the differential geometry and total differential theory to the field of InSAR adjustment and correct the spacial baseline, slant range, and timing parameters to improve the three-dimensional accuracy of InSAR-generated DEMs. Figure 5 shows the derivation process workflow of the proposed method. Our derivation process can be divided into three steps. The first step is to use the total differential theory to calculate the explicit expression of InSAR three-dimensional positioning errors concerning the errors of physical parameters in the ECEF coordinate system, similar to Equation (5). The second step is to use the differential geometry theory to calculate the northern, eastern, and upper three-dimensional unit vectors corresponding to the current ECEF coordinates of ground object points and project the geolocation errors obtained in the first step in these three directions, the control directions of various control data. This step is designed to enable the proposed method to utilize the advantages of various control data comprehensively. For example, the Height Control Points (HCPs) have accurate elevation information and can be used to correct the systematic elevation error of DEMs. Consequently, the InSAR geolocation errors in the ECEF coordinate need to be converted to the elevation directions to be consistent with the control directions of HCPs. The third step is to list adjustment equations according to the characteristics of various control data. For example, we can discard the adjustment equations in their planes (northern and eastern directions) for height control points and only use their height (in upper directions) equations.

In this section, according to Figure 5, we first derive the explicit expression of the three-dimensional geolocation error, i.e., $d\overrightarrow{T}=f(dR,dt,d\overrightarrow{B})$. Then, we further derive adjustment equations for HCPs, Plane Control Points (PCPs), Height Tie Points (HTPs), and Plane Tie Points (PTPs) employing the vector projection method. Finally, we propose an iterative scheme to solve adjustment parameters.

3.1. Three-Dimensional Geolocation Error Equation

Before analyzing various control data, we need to establish the fundamental adjustment equation, solving the integrated three-dimensional geolocation error of Equation (1) caused by all the parameter errors. First, the total differential of Equation (1) was calculated, as follows.

$$\left\{\begin{array}{cc}\hfill & (\overrightarrow{T}-\overrightarrow{P_m})·d\overrightarrow{T}=RdR+(\overrightarrow{T}-\overrightarrow{P_m})·d\overrightarrow{P_m}\hfill \\ \hfill & \overrightarrow{V_m}·d\overrightarrow{T}=(\overrightarrow{P_m}-\overrightarrow{T})·d\overrightarrow{V_m}+\overrightarrow{V_m}·d\overrightarrow{P_m}+\frac{1}{2}\lambda f_{dc}dR\hfill \\ \hfill & (\overrightarrow{T}-\overrightarrow{P_s})·d\overrightarrow{T}=(R+\Delta R)dR+(\overrightarrow{T}-\overrightarrow{P_s})·(d\overrightarrow{P_m}+d\overrightarrow{B})\hfill \end{array}\right.$$

(10)

Based on the analysis in Section 2, we model the parallel baseline error $d{B}_{\Vert }$ as a polynomial about azimuth time t. Then the three-dimensional baseline error $d\overrightarrow{B}$ can be expressed as the product of $d{B}_{\Vert }$ and the unit LOS vector of the master satellite $\overrightarrow{l_m}$.

$$\left\{\begin{array}{cc}\hfill d{B}_{\Vert }& =\sum _{n=0}^{N_b}{b}_n{t}^n\hfill \\ \hfill d\overrightarrow{B}& =d{B}_{\Vert }\overrightarrow{l_m}=\overrightarrow{l_m}\sum _{n=0}^{N_b}{b}_n{t}^n\hfill \end{array}\right.$$

(11)

where $\overrightarrow{l_m}$ is the unit direction vector of $\overrightarrow{T}-\overrightarrow{P_m}$. $N_b$ and $b_n$ are polynomial order and coefficients of $d{B}_{\Vert }$, respectively. In order to avoid the Runge phenomenon, we can select a low-order polynomial to fit the baseline error, for example, $N_b=1$. Then, let $\overrightarrow{L_m}=\overrightarrow{T}-\overrightarrow{P_m}$, $\overrightarrow{L_s}=\overrightarrow{T}-\overrightarrow{P_s}$ and $R_s=R+\Delta R$, Equation (10) can be further deduced as

$$\left\{\begin{array}{cc}\hfill & \overrightarrow{L_m}·d\overrightarrow{T}=RdR+\overrightarrow{L_m}·\frac{d\overrightarrow{P_m}}{dt}dt\hfill \\ \hfill & \overrightarrow{V_m}·d\overrightarrow{T}=\frac{1}{2}\lambda f_{dc}dR+\left(\overrightarrow{V_m}·\frac{d\overrightarrow{P_m}}{dt}-\overrightarrow{L_m}·\frac{d\overrightarrow{V_m}}{dt}\right)dt\hfill \\ \hfill & \overrightarrow{L_s}·d\overrightarrow{T}=R_{s}dR+\overrightarrow{L_s}·\left(\frac{d\overrightarrow{P_m}}{dt}+\overrightarrow{l_m}\sum _{n=1}^{N_b}n{b}_n{t}^{n-1}\right)dt+\overrightarrow{L_s}·\overrightarrow{l_m}\sum _{n=0}^{N_b}{t}^{n}d{b}_n\hfill \end{array}\right.$$

(12)

Appendix B.1 (in Appendix B) shows part of the derivation of Equation (12). Similar to Equation (3), Equation (12) can still be expressed as $\overline{\overline{D}}·d\overrightarrow{T}=\overrightarrow{u}$. Let

$${\overline{\overline{D}}}^{-1}={\left[\begin{array}{ccc}X-P_{m,x}& Y-P_{m,y}& Z-P_{m,z}\\ V_x& V_y& V_z\\ X-P_{s,x}& Y-P_{s,y}& Z-P_{s,z}\end{array}\right]}^{-1}=\left[\begin{array}{ccc}{d}_{11}& d_{12}& d_{13}\\ d_{21}& d_{22}& d_{23}\\ d_{31}& d_{32}& d_{33}\end{array}\right]$$

(13)

then we can get the specific form of $d\overrightarrow{T}=(dX,dY,dZ)$, as Equation (14). On the whole, the adjustment parameter is $\mathbf{x}={(d{R}_1,d{T}_1,d{b}_{1,0},\cdots ,d{b}_{1,N_b},d{R}_2,\cdots )}^{\mathrm{T}}$.

$$\left\{\begin{array}{cc}\hfill & dX=\left(d_{11}R-d_{12}\frac{1}{2}\lambda f_{dc}-d_{13}{R}_s\right)dR+d_{13}\overrightarrow{L_s}·\overrightarrow{l_m}\sum _{n=0}^{N_b}{t}^{n}d{b}_n\hfill \\ \hfill & +\left[d_{11}\overrightarrow{L_m}·\frac{d\overrightarrow{P_m}}{dt}+d_{12}\left(\overrightarrow{V_m}·\frac{d\overrightarrow{P_m}}{dt}-\overrightarrow{L_m}·\frac{d\overrightarrow{V_m}}{dt}\right)+d_{13}\overrightarrow{L_s}·\left(\frac{d\overrightarrow{P_m}}{dt}+\overrightarrow{l_m}\sum _{n=1}^{N_b}n{b}_n{t}^{n-1}\right)\right]dt\hfill \\ \hfill & dY=\left(d_{21}R-d_{22}\frac{1}{2}\lambda f_{dc}-d_{23}{R}_s\right)dR+d_{23}\overrightarrow{L_s}·\overrightarrow{l_m}\sum _{n=0}^{N_b}{t}^{n}d{b}_n\hfill \\ \hfill & +\left[d_{21}\overrightarrow{L_m}·\frac{d\overrightarrow{P_m}}{dt}+d_{22}\left(\overrightarrow{V_m}·\frac{d\overrightarrow{P_m}}{dt}-\overrightarrow{L_m}·\frac{d\overrightarrow{V_m}}{dt}\right)+d_{23}\overrightarrow{L_s}·\left(\frac{d\overrightarrow{P_m}}{dt}+\overrightarrow{l_m}\sum _{n=1}^{N_b}n{b}_n{t}^{n-1}\right)\right]dt\hfill \\ \hfill & dZ=\left(d_{31}R-d_{32}\frac{1}{2}\lambda f_{dc}-d_{33}{R}_s\right)dR+d_{33}\overrightarrow{L_s}·\overrightarrow{l_m}\sum _{n=0}^{N_b}{t}^{n}d{b}_n\hfill \\ \hfill & +\left[d_{31}\overrightarrow{L_m}·\frac{d\overrightarrow{P_m}}{dt}+d_{32}\left(\overrightarrow{V_m}·\frac{d\overrightarrow{P_m}}{dt}-\overrightarrow{L_m}·\frac{d\overrightarrow{V_m}}{dt}\right)+d_{33}\overrightarrow{L_s}·\left(\frac{d\overrightarrow{P_m}}{dt}+\overrightarrow{l_m}\sum _{n=1}^{N_b}n{b}_n{t}^{n-1}\right)\right]dt\hfill \end{array}\right.$$

(14)

3.2. Adjustment Equations of Various Control Data

After obtaining Equation (14), we have completed the first step in Figure 5, i.e., solving the relationship between three-dimensional geolocation errors and physical parameters’ errors in the ECEF coordinate system. Equation (14) is linear about each error source because of the total differential theory used in the solution process and can be written as follows.

$$d\overrightarrow{T}=\left(\begin{array}{c}dX\\ dY\\ dZ\end{array}\right)=\left(\begin{array}{c}{c}_{11}dR+c_{21}dt+c_{31}\sum _{n=0}^{N_b}{t}^{n}d{b}_n\\ c_{12}dR+c_{22}dt+c_{32}\sum _{n=0}^{N_b}{t}^{n}d{b}_n\\ c_{13}dR+c_{23}dt+c_{33}\sum _{n=0}^{N_b}{t}^{n}d{b}_n\end{array}\right)=\overrightarrow{c_1}dR+\overrightarrow{c_2}dt+\overrightarrow{c_3}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_n\right)$$

(15)

where $\overrightarrow{c_i}=(c_{i,1},c_{i,2},c_{i,3}),(i=1,2,3)$ is the vector representation of coefficients in Equation (14), which is known during the adjustment. For the above equation, we need to make a few additional explanations. First, the derivations of Equation (14) are all in the ECEF coordinate system and do not involve the conversion of different coordinate systems. Second, Equation (14) is actually the adjustment equation for high-precision three-dimensional control points, such as corner reflectors. Third, as we discussed earlier, the unknown parameters in Equation (14) are the physical parameter errors, including $dR$, $dt$, and $d{b}_n$, which need to be determined by control data.

Generally speaking, the available control data have different characteristics. For example, the data with high elevation but low plane accuracy can be employed as HCPs. Therefore, Equation (14) needs to be further modified by geometric projection, the second and third steps in Figure 5.

3.2.1. Equation of HCPs

The HCP is a kind of data with accurate elevations, which is used to improve the absolute elevation accuracy of interferometric DEMs. The global public laser altimetry data of the Ice, Cloud, and land Elevation Satellite (ICESat) is an essential source of HCPs, the elevation accuracy of which is better than 1.0 m [32], has been widely used in the field of non-parametric adjustment [5,7,9]. However, this data does not have recognizable image information and cannot be matched with InSAR image or DEM data, so it cannot provide plane control information. Therefore, we hope that the height errors (projection of three-dimensional positioning errors in height direction) caused by the InSAR physical parameters’ errors equal deviations between InSAR-generated DEMs and HCPs, as Equation (16).

$$\begin{array}{cc}\hfill \Delta T_{u,hc}& =h_{HCP}-h_{DEM}=d\overrightarrow{T}·\overrightarrow{n_u}\hfill \\ \hfill & =\overrightarrow{n_u}·\overrightarrow{c_1}dR+\overrightarrow{n_u}·\overrightarrow{c_2}dt+\overrightarrow{n_u}·\overrightarrow{c_3}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_n\right)\hfill \end{array}$$

(16)

where $h_{HCP}$ and $h_{DEM}$ represent the elevation of HCP and the corresponding pixel in DEM, respectively, and $\overrightarrow{n_u}$ is the upper unit vector at a ground target T (please refer to Appendix A). Equation (16) can be further summarized in the form of an error equation, as,

$$v_{hc}=\overrightarrow{n_u}·\overrightarrow{c_1}dR+\overrightarrow{n_u}·\overrightarrow{c_2}dt+\overrightarrow{n_u}·\overrightarrow{c_3}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_n\right)-\Delta T_{u,hc}$$

(17)

where $v_{hc}$ is the residual error, which is 0 in the ideal case. For multi-scene data and multiple groups of HCPs, Equation (17) can be further written in matrix form.

$${\mathbf{V}}_{hc}={\mathbf{A}}_{hc}\mathbf{x}-\Delta {\mathbf{T}}_{hc}$$

(18)

where ${\mathbf{A}}_{hc}$ is the design matrix composed of coefficient $\overrightarrow{n_u}·\overrightarrow{c_i}(i=1,2,3)$, and its specific form is shown in Appendix B.2. Assuming that the number of adjustment InSAR tracks is N and the polynomial order of baseline error is $N_b$, then the parameters to be adjusted, namely, the corrections of slant range, azimuth time, and parallel baseline in each scene data, is $\mathbf{x}=[d{R}_1,d{t}_1,d{b}_{1,0},\cdots ,d{b}_{1,N_b},d{R}_2,\cdots ,d{b}_{N,N_b}]$. In the following text, $\mathbf{x}$ has the same meaning as Equation (18). $\Delta {\mathbf{T}}_{hc}$ represents the elevation deviation vector composed of $\Delta T_{h,hc}$ in Equation (16). ${\mathbf{V}}_{hc}$ is the residual vector. In the subsequent text, we will omit the error equation as Equation (17).

3.2.2. Equation of PCPs

The PCP is a kind of control data with high plane accuracy, which is used to improve the absolute plane accuracy of DEM. We can obtain PCPs from globally distributed DEM/DSM data (such as AW3D30) [33,34] or optical image data (such as Google Maps) [35] by homologous or heterogeneous matching methods. The elevation quality of such data is usually inferior to interferometric DEMs, so it is difficult to use as elevation control data. Therefore, we hope that the plane errors (projection of three-dimensional positioning errors in northern and eastern directions) caused by the InSAR physical parameters’ errors equal deviations between InSAR-generated DEMs and PCPs, as

$$\left\{\begin{array}{cc}\hfill \Delta T_{n,pc}& =y_{PCP}-y_{DEM}=d\overrightarrow{T}·\overrightarrow{n_n}\hfill \\ \hfill & =\overrightarrow{n_n}·\overrightarrow{c_1}dR+\overrightarrow{n_n}·\overrightarrow{c_2}dt+\overrightarrow{n_n}·\overrightarrow{c_3}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_n\right)\hfill \\ \hfill \Delta T_{e,pc}& =x_{PCP}-x_{DEM}=d\overrightarrow{T}·\overrightarrow{n_e}\hfill \\ \hfill & =\overrightarrow{n_e}·\overrightarrow{c_1}dR+\overrightarrow{n_e}·\overrightarrow{c_2}dt+\overrightarrow{n_e}·\overrightarrow{c_3}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_n\right)\hfill \end{array}\right.$$

(19)

where $(x,y)$ is the plane coordinate of DEMs or PCPs. Similar to Equation (18), Equation (19) can also be written as

$${\mathbf{V}}_{pc}={\mathbf{A}}_{pc}\mathbf{x}-\Delta {\mathbf{T}}_{pc}$$

(20)

where ${\mathbf{A}}_{pc}$ is the design matrix composed of coefficient $\overrightarrow{n_n}·\overrightarrow{c_i}$ and $\overrightarrow{n_e}·\overrightarrow{c_i}(i=1,2,3)$. $\Delta {\mathbf{T}}_{pc}$ represents the horizontal deviation vector composed of $\Delta T_{n,pc}$ and $\Delta T_{e,pc}$ in Equation (19). ${\mathbf{V}}_{pc}$ is the residual vector.

3.2.3. Equation of HTPs

In order to ensure the relative elevation accuracy of elevation DEM, M. Huber et al. put forward the concept of Tie Point chips [13], that is, the elevation chips located in different images corresponding to the same area. The core idea of TP chips is to calculate the median elevation in an overlapping area as the observation value to avoid height outliers [5,7]. In this article, we call such data as HTPs. They can be obtained only through geolocation calculation without image-matching methods. Therefore, such data cannot guarantee the relative plane accuracy of DEMs. We hope the corrected DEMs have identical elevation values in overlapping areas during adjustment, that is,

$$h_{DEM,J}+d\overrightarrow{T_J}·\overrightarrow{n_n}=h_{DEM,I}+d\overrightarrow{T_I}·\overrightarrow{n_n}$$

(21)

where I and J are the indexes of DEMs. Then Equation (21) can be further transformed into

$$\begin{array}{cc}\hfill \Delta T_{u,ht}& =h_{DEM,J}-h_{DEM,I}=d\overrightarrow{T_I}·\overrightarrow{n_u}-d\overrightarrow{T_J}·\overrightarrow{n_u}\hfill \\ \hfill & =\overrightarrow{n_u}·\overrightarrow{c_{1,I}}d{R}_I+\overrightarrow{n_u}·\overrightarrow{c_{2,I}}d{t}_I+\overrightarrow{n_u}·\overrightarrow{c_{3,I}}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_{n,I}\right)\hfill \\ \hfill & -\overrightarrow{n_u}·\overrightarrow{c_{1,J}}d{R}_J-\overrightarrow{n_u}·\overrightarrow{c_{2,J}}d{t}_J-\overrightarrow{n_u}·\overrightarrow{c_{3,J}}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_{n,J}\right)\hfill \end{array}$$

(22)

Here, $\Delta T_{h,ht}$ represents the elevation deviation of the same ground point in different InSAR-generated DEMs. For multi-scene data and multiple groups of HTPs, Equation (22) can be further written in matrix form.

$${\mathbf{V}}_{ht}={\mathbf{A}}_{ht}\mathbf{x}-\Delta {\mathbf{T}}_{ht}$$

(23)

where ${\mathbf{A}}_{ht}$ is the design matrix composed of the coefficient $\overrightarrow{n_u}·\overrightarrow{c_i}(i=1,2,3)$ of different points corresponding to a same ground object, and its specific form is shown in Appendix B.2. $\Delta {\mathbf{T}}_{ht}$ represents the elevation deviation vector composed of $\Delta T_{h,ht}$ in Equation (22). ${\mathbf{V}}_{ht}$ is the residual vector.

3.2.4. Equation of PTPs

The PTPs are the tie points in our general sense: the pixel points located in different images and corresponding to the same ground object [36,37,38]. It is used to constrain the relative plane accuracy of DEMs and can be obtained by homologous-image matching methods. However, such points are usually located at the boundary area of imaging data. Compared with the internal area, their elevations may be polluted by heavier noise [7,13], so they cannot be used to ensure the relative elevation accuracy of DEMs. After adjustment, we hope different DEMs can get consistent plane positioning results at the PTP position, that is,

$$\left\{\begin{array}{cc}\hfill & x_{DEM,I}+d\overrightarrow{T_I}·\overrightarrow{n_e}=x_{DEM,J}+d\overrightarrow{T_J}·\overrightarrow{n_e}\hfill \\ \hfill & y_{DEM,I}+d\overrightarrow{T_I}·\overrightarrow{n_n}=y_{DEM,J}+d\overrightarrow{T_J}·\overrightarrow{n_n}\hfill \end{array}\right.$$

(24)

The above equation can be further transformed into

$$\left\{\begin{array}{cc}\hfill \Delta T_{e,pt}& =x_{DEM,J}-x_{DEM,I}=d\overrightarrow{T_I}·\overrightarrow{n_e}-d\overrightarrow{T_J}·\overrightarrow{n_e}\hfill \\ \hfill & =\overrightarrow{n_e}·\overrightarrow{c_{1,I}}d{R}_I+\overrightarrow{n_e}·\overrightarrow{c_{2,I}}d{t}_I+\overrightarrow{n_e}·\overrightarrow{c_{3,I}}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_{n,I}\right)\hfill \\ \hfill & -\overrightarrow{n_e}·\overrightarrow{c_{1,J}}d{R}_J-\overrightarrow{n_e}·\overrightarrow{c_{2,J}}d{t}_J-\overrightarrow{n_e}·\overrightarrow{c_{3,J}}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_{n,J}\right)\hfill \\ \hfill \Delta T_{n,pt}& =y_{DEM,J}-y_{DEM,I}=d\overrightarrow{T_I}·\overrightarrow{n_n}-d\overrightarrow{T_J}·\overrightarrow{n_n}\hfill \\ \hfill & =\overrightarrow{n_n}·\overrightarrow{c_{1,I}}d{R}_I+\overrightarrow{n_n}·\overrightarrow{c_{2,I}}d{t}_I+\overrightarrow{n_e}·\overrightarrow{c_{3,I}}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_{n,I}\right)\hfill \\ \hfill & -\overrightarrow{n_n}·\overrightarrow{c_{1,J}}d{R}_J-\overrightarrow{n_n}·\overrightarrow{c_{2,J}}d{t}_J-\overrightarrow{n_n}·\overrightarrow{c_{3,J}}\left(\sum _{n=0}^{N_b}{t}^{n}d{b}_{n,J}\right)\hfill \end{array}\right.$$

(25)

Here, $\Delta T_{e,pt}$ and $\Delta T_{n,pt}$ are the horizontal deviations of the same ground point in different InSAR-generated DEMs. Similar to Equation (23), Equation (25) can also be written as

$${\mathbf{V}}_{pt}={\mathbf{A}}_{pt}\mathbf{x}-\Delta {\mathbf{T}}_{pt}$$

(26)

where ${\mathbf{A}}_{pt}$ is the design matrix composed of coefficient $\overrightarrow{n_n}·\overrightarrow{c_i}$ and $\overrightarrow{n_e}·\overrightarrow{c_i}(i=1,2,3)$ of different points corresponding to a same ground object. $\Delta {\mathbf{T}}_{pt}$ represents the horizontal deviation vector composed of $\Delta T_{n,pt}$ and $\Delta T_{e,pt}$ in Equation (25). ${\mathbf{V}}_{pt}$ is the residual vector.

3.3. Solution Strategy

In the previous subsection, we have derived the adjustment equations of various control data and given the corresponding matrix expressions. It can be found that the residual vectors in Equations (18), (20), (23) and (26) represent the vertical or horizontal positioning residuals, the units of which are meters. Therefore, we consider unifying all equations into a loss function to solve the adjustment parameter $\mathbf{x}$. We give the loss function as follows,

$$\begin{array}{cc}\hfill L\left(\mathbf{x}\right)& =q_{hc}\parallel {\mathbf{V}}_{hc}{\parallel }_2^2+q_{pc}\parallel {\mathbf{V}}_{pc}{\parallel }_2^2+q_{ht}\parallel {\mathbf{V}}_{ht}{\parallel }_2^2+q_{pt}\parallel {\mathbf{V}}_{pt}{\parallel }_2^2+\mu {\parallel \mathbf{x}\parallel }_2^2\hfill \\ \hfill & =q_{hc}\parallel {\mathbf{A}}_{hc}\mathbf{x}-\Delta {\mathbf{T}}_{hc}{\parallel }_2^2+q_{pc}{\parallel {\mathbf{A}}_{pc}\mathbf{x}-\Delta {\mathbf{T}}_{pc}\parallel }_2^2\hfill \\ \hfill & +q_{ht}\parallel {\mathbf{A}}_{ht}\mathbf{x}-\Delta {\mathbf{T}}_{ht}{\parallel }_2^2+q_{pt}\parallel {\mathbf{A}}_{pt}\mathbf{x}-\Delta {\mathbf{T}}_{pt}{\parallel }_2^2+\mu {\parallel \mathbf{x}\parallel }_2^2\hfill \end{array}$$

(27)

where $\mathbf{x}=[d{R}_1,d{t}_1,d{b}_{1,0},\cdots ,d{b}_{1,N_b},d{R}_2,\cdots ,d{b}_{N,N_b}]$ is the parameters to be adjusted and N is the number of adjustment InSAR tracks. ${\mu \parallel \mathbf{x}\parallel }_2^2$ is a regular term in the loss function to avoid the rank defect problem in the optimization process. When the control data are sufficient and there is no rank deficiency, $\mu$ can be set as zero. Otherwise, $\mu$ can be determined by the ridge mark method [39]. $q_{hc}$, $q_{pc}$, $q_{ht}$, and $q_{pt}$ are the weights of control data, which are determined by the number $N_{hc}$, $N_{pc}$, $N_{ht}$, $N_{pt}$, and the prior accuracy ${\sigma }_{hc}$, ${\sigma }_{pc}$, ${\sigma }_{ht}$, ${\sigma }_{pt}$ of control data, as follows

$$\left\{\begin{array}{cc}\hfill q_{hc}& ={\left({\sigma }_{hc}^2{N}_{hc}\right)}^{-1}\hfill \\ \hfill q_{pc}& ={\left({\sigma }_{pc}^2{N}_{pc}\right)}^{-1}\hfill \\ \hfill q_{ht}& ={\left({\sigma }_{ht}^2{N}_{ht}\right)}^{-1}\hfill \\ \hfill q_{pt}& ={\left({\sigma }_{pt}^2{N}_{pt}\right)}^{-1}\hfill \end{array}\right.$$

(28)

The control data should be as much as possible to ensure the absolute and relative accuracy of DEMs. The adjustment parameter $\mathbf{x}$ can be solved by minimizing the loss function, equivalent to calculating the gradient zero point of $L\left(\mathbf{x}\right)$, as follows.

$$\underset{\mathbf{x}}{\min }L\left(\mathbf{x}\right)\iff {\nabla }_{\mathbf{x}}L\left(\mathbf{x}\right)=\mathbf{0}$$

(29)

The gradient function of $L\left(\mathbf{x}\right)$ is

$$\begin{array}{cc}\hfill {\nabla }_{\mathbf{x}}L\left(\mathbf{x}\right)& =q_{hc}{\mathbf{A}}_{hc}^{\mathrm{T}}{\mathbf{A}}_{hc}\mathbf{x}+q_{pc}{\mathbf{A}}_{pc}^{\mathrm{T}}{\mathbf{A}}_{pc}\mathbf{x}+q_{ht}{\mathbf{A}}_{ht}^{\mathrm{T}}{\mathbf{A}}_{ht}\mathbf{x}+q_{pt}{\mathbf{A}}_{pt}^{\mathrm{T}}{\mathbf{A}}_{pt}\mathbf{x}+\mu \mathbf{x}\hfill \\ \hfill & -q_{hc}{\mathbf{A}}_{hc}^{\mathrm{T}}\Delta {\mathbf{T}}_{hc}-q_{pc}{\mathbf{A}}_{pc}^{\mathrm{T}}\Delta {\mathbf{T}}_{pc}-q_{ht}{\mathbf{A}}_{ht}^{\mathrm{T}}\Delta {\mathbf{T}}_{ht}-q_{pt}{\mathbf{A}}_{pt}^{\mathrm{T}}\Delta {\mathbf{T}}_{pt}\hfill \end{array}$$

(30)

Let

$$\left\{\begin{array}{cc}\hfill \mathbf{A}& =q_{hc}{\mathbf{A}}_{hc}^{\mathrm{T}}{\mathbf{A}}_{hc}+q_{pc}{\mathbf{A}}_{pc}^{\mathrm{T}}{\mathbf{A}}_{pc}+q_{ht}{\mathbf{A}}_{ht}^{\mathrm{T}}{\mathbf{A}}_{ht}+q_{pt}{\mathbf{A}}_{pt}^{\mathrm{T}}{\mathbf{A}}_{pt}+\mu \mathbf{I}\hfill \\ \hfill \mathbf{L}& =q_{hc}{\mathbf{A}}_{hc}^{\mathrm{T}}\Delta {\mathbf{T}}_{hc}+q_{pc}{\mathbf{A}}_{pc}^{\mathrm{T}}\Delta {\mathbf{T}}_{pc}+q_{ht}{\mathbf{A}}_{ht}^{\mathrm{T}}\Delta {\mathbf{T}}_{ht}+q_{pt}{\mathbf{A}}_{pt}^{\mathrm{T}}\Delta {\mathbf{T}}_{pt}\hfill \end{array}\right.$$

(31)

where $\mathbf{I}$ is unit matrix. Then the Least-Square (LS) solution of Equation (29) is

$$\mathbf{x}={\mathbf{A}}^{-1}\mathbf{L}$$

(32)

Since rigorous model Equation (1) is nonlinear, we adopt the iterative least squares method to solve adjustment parameters. The specific steps are as follows:

i.

Initialize adjustment parameter $\mathbf{x}$, i.e., ${\mathbf{x}}_0=\mathbf{0}$.

ii.

Set tolerance limits $to{l}_B$, $to{l}_R$, and $to{l}_T$ for the increment of baseline, range, and timing parameters. According to the analysis in Section 2, $to{l}_B$, $to{l}_R$, and $to{l}_T$ should be below millimeters, meters, and 0.1 milliseconds, respectively. For example, $to{l}_B=0.1$ mm, $to{l}_R=0.1$ m, and $to{l}_T=0.01$ ms.

iii.

List the loss function $L_k\left(\mathbf{x}\right)$ as Equation (27), and solve the adjustment parameter increment $\delta {\mathbf{x}}_k$ using Equation (32), i.e., $\delta {\mathbf{x}}_k={\mathbf{A}}_k^{-1}{\mathbf{L}}_k$, where k represents the kth iteration.

iv.

Sperate the corrections of slant range, azimuth time, and spatial baseline according to their track, add the corrections to the result of the previous iteration, i.e., ${\mathbf{x}}_k={\mathbf{x}}_{k-1}+\delta {\mathbf{x}}_k$, to update the range, timing, and baseline parameters.

v.

Judge whether the increment of each parameter is less than the tolerance limit. If less than, end iterative calculation and take ${\mathbf{x}}_k$ as the solution of adjustment parameters. Otherwise, re-execute step iii.

4. Experiment and Discussion

4.1. Experiment Data

We conduct experiments with twenty-nine scenes of TerraSAR-X/TanDEM-X Co-registered Single look Slant range Complex (CoSSC) data, including twenty ascending data and nine descending data. The specific information on InSAR data is shown in

Table 2. Specific information of CoSSC data. The table is listed by track.

Time	Orbit	Incidence (${}^{\circ }$)	Resolution (m)	Baseline Length (m)
26 January 2012	ASC	33.63	1.34 × 1.98	506.15
6 February 2012	ASC	34.16	1.36 × 2.03	577.90
17 February 2012	ASC	33.58	1.36 × 2.19	607.98
2 November 2012	ASC	33.91	1.36 × 2.04	485.38
5 December 2012	ASC	33.86	1.36 × 1.99	657.02
16 December 2012	ASC	33.83	1.36 × 1.71	612.14
24 July 2013	ASC	33.95	1.36 × 1.85	559.86
27 January 2018	DEC	34.27	1.36 × 2.19	479.02
27 February 2019	DEC	38.35	1.36 × 1.83	459.50
10 March 2019	DEC	33.72	1.36 × 1.98	480.65

. The experimental area covers Xuchang and Zhengzhou, Henan Province, China, covering an area of about 56,000 square kilometers with rich geomorphic features. The Yellow River flows through the north, the Songshan Mountains are located in the west, and most experimental areas are plains. Figure 6 shows the location and elevation of experimental data, respectively. The HCPs used in the experiment are ICESat laser altimetry data GLAH14 selected according to the extreme criteria proposed in [9], with an elevation accuracy of about 0.2 m. The height checkpoints are from ICESat-2 ATL08 data (Land and Vegetation Height), with a vertical accuracy of 2.0 m and 0.2 m, respectively, in mountainous areas and flat lands [40]. The PCPs and PTPs are matched using the DEM matching method [33,34], where PCPs are from the public SRTM-DEM.

4.2. Simulated Experiment

The simulated experiment is designed to verify the performance of the proposed method in the noise environment. The primary benchmarks are the adjustment parameters and three-dimensional accuracy after adjustment. The specific steps of the simulated experiment are as follows:

i.

Generate absolute interferometric phase. We select the external DEM in the experimental area as the simulated terrain and the geometric parameters of InSAR data, such as orbit and baseline, as the accurate geometric parameters. The absolute interferometric phase is calculated through the reverse geolocation method.

ii.

Generate checkpoints and control points. We select a series of geographic points in the simulated terrain, part of which are checkpoints, and the other points are added with Gaussian white noise as the HCPs and PCPs, respectively.

iii.

Generate DEMs. We add systematic error into baseline, slant range, and timing parameters and then generate DEMs according to the strict model Equation (1). Then, we determined the observation values of heights corresponding to HCPs and PCPs extracted in step ii.

iv.

Extract tie points. As described in Section 3.2.3 and Section 3.2.4, HTPs can be selected by the geolocation method, while PTPs need to be extracted by matching methods.

v.

Adjust and verify accuracy. Unlike the real data experiment, the accuracy verification of the simulated experiment includes not only the plane and elevation accuracies of adjusted DEMs but also the difference between the simulated error and the adjustment calculation value. The latter is used to verify whether the proposed method can accurately calculate the systematic errors of bistatic InSAR systems.

Figure 7 shows the flow chart of the simulated experiment. Generally speaking, the simulated data experiment consists of three parts: InSAR processing, control data extraction, and adjustment.

In practice, high-precision height control points (ICESat data) are more accessible, while plane control points are more challenging to obtain. Therefore, the simulated experiment takes the noise of PCPs as a variable, controls its standard deviation within 0 to 10.0 m, and conducts ten groups of experiments. Each group is composed of thirty Monte Carlo simulated experiments. That is, under the same noise level, we have conducted thirty simulated experiments. We use the Root-Mean-Square Error (RMSE) to measure the adjustment accuracy. Figure 8 shows the results of the simulated experiment. Figure 8a,b show the height and plane accuracies after adjustment. Figure 8c–e, respectively, show the errors of parallel baseline, slant range, and timing parameters obtained by our proposed method. In Figure 8, the blue lines are the curves of median RMSEs versus noise standard deviation after adjustment and the point represent the RMSEs of all Monte Carlo experiments.

When the PCP noise changes from 0.0 m to 10.0 m, the RMSE median of each parameter changes steadily, demonstrating that our proposed methods can obtain stable results in the simulated environment. From the perspective of median RMSE, the adjustment errors of spacial baseline, slant range, and timing parameters after adjustment are, respectively, stable at about 0.05 mm, 0.1 m, and 0.006 ms, and the plane and elevation accuracies are, respectively, stable at about 0.2 m. Under the same noise level, some simulated experiments may use the outlier control data for adjustment so that the results will have a fluctuation range. Even so, the adjustment accuracy is still high, and RMSEs of elevation and plane are stable below 1.0 m, which indicates that the method proposed by us has an excellent adjustment effect.

In order to verify the influence of the number of control data on the adjustment accuracy, we designed a simulated experiment. The core idea of this experiment is to gradually increase the number of PCPs and observe the corresponding adjustment accuracy. The experimental framework is consistent with Figure 7. Figure 9 shows the simulated experiment results for the number of PCPs. It can be seen that when the number of PCPs contained in each track changes from 1 to 10, the estimated errors of the adjustment parameters, that is, baseline, range, and azimuth time, gradually decrease. For example, when the number of PCPs is 1, the baseline estimated error is about 0.14 mm. When PCPs exceed 5 per track, the estimation error is stable at about 0.06 mm. With the increase in the number of PCPs, the plane accuracy of DEM has gradually improved from 0.6 m to 0.2 m. Since the number of HCPs has not changed, the elevation accuracy of ten groups of experiments is stable at about 0.15 m, and there is no clear law with the number of PCPs.

In addition, Figure 10 shows the convergence curve of an iterative calculation process. In general, after three to five iterations, the increments of spatial baseline, slant range, and timing parameters all tend to zero, which proves that the proposed method has strong convergence.

4.3. Real-Data Experiment

In this part, we first introduce the elevation and plane fracture phenomenon in DEM mosaic maps caused by systematic geolocation errors before adjustment. Then, we show the three-dimensional accuracies of InSAR-generated DEMs after adjustment and compare them with the non-parametric adjustment method.

4.3.1. Fracture in DEM Mosaic Maps

We use Global Mapper to mosaic twenty-nine InSAR-generated DEMs. Figure 11a is the rendered mosaic map before adjustment. In order to visually illustrate the fracture phenomenon of DEM in the overlapping area, we have selected two typical areas in Figure 11a for magnification. Figure 11b is the gradient map of marked area (b), and Figure 11c is an enlarged view of marked area (c). Here, the gradient map highlights the elevation difference of DEMs in the mosaic area and clearly shows the elevation fracture phenomenon. The slope of DEMs are calculated as follows:

$$s={tan}^{-1}\sqrt{h_x^2+h_y^2};$$

(33)

$h_x$ and $h_y$ are gradients in two orthogonal directions, which can be calculated by filtering operators [41]. It can be seen that elevation fracture will occur in the DEM overlapping area, which is shown as a straight line of the mosaic boundary in the slope map. A similar phenomenon occurs in the plane. In Figure 11c, a railway runs through the southeast and northwest directions, so the elevation in this area should be continuous on the plane. However, there are apparent fractures in the red-marked area in Figure 11c, which shows that the eastern DEM has shifted a certain distance to the north. In Figure 11c, the red arrows indicate the position and direction of the plane offset. Therefore, it is necessary to adopt the three-dimensional adjustment method to eliminate systematic errors of interferometric DEMs.

4.3.2. Adjustment Effect

We obtain 1225 HCPs, 63 PCPs, 938 pairs of HTPs, and 947 pairs of PTPs to conduct the real data experiment. The maximum shooting duration of the CoSSC data is about 22 s. During this period, the baseline error can be considered a constant [22,23,24]. Therefore, we set the order of the baseline error polynomial to zero. We use RMSEs of height and plane checkpoint to measure three-dimensional geolocation accuracy, and the expressions are as follows.

$$\left\{\begin{array}{cc}\hfill {\mathrm{RMSE}}_h& =\sqrt{\sum _{i=1}^{N_{hckp}}\frac{{(h_{DEM,i}-h_{CKPs,i})}^2}{N_{hckp}}}\hfill \\ \hfill {\mathrm{RMSE}}_p& =\sqrt{\sum _{i=1}^{N_{pckp}}\frac{{(x_{DEM,i}-x_{CKPs,i})}^2+{(y_{DEM,i}-y_{CKPs,i})}^2}{N_{pckp}}}\hfill \end{array}\right.$$

(34)

where ${\mathrm{RMSE}}_h$ and ${\mathrm{RMSE}}_p$ are the vertical and horizontal RMSEs, respectively. h is the elevation of height checkpoints or DEM points. $(x,y)$ are the plane coordinates of plane checkpoints or DEM points. In addition, we use the nonparametric method to adjust the elevation of InSAR-generated DEMs as a comparative experiment of the proposed method.

Table 3. The adjustment results of InSAR parameters. The table is listed by track.

Data	Range (m)	Timing (ms)	Basline Offset (mm)
26 January 2012	5.47	0.109	2.03
6 February 2012	4.40	0.054	1.51
17 February 2012	5.42	−0.28	2.34
2 November 2012	2.01	−0.641	1.99
5 December 2012	−14.84	−0.459	1.26
16 December 2012	0.75	−0.040	0.25
24 July 2013	2.61	1.40	0.10
27 January 2018	2.47	0.937	2.12
27 February 2019	1.87	0.540	0.20
10 March 2019	0.95	0.734	1.14

shows the parameter correction after adjustment. The corrections of parallel baseline are less than 3 mm, which is consistent with the current TanDEM-X baseline calibration results. Figure 12a,b are the elevation and plane accuracy counted by tracks before and after adjustment, respectively. The cyan bars in Figure 12a represent the adjustment accuracy obtained by the general method. Before the adjustment, the elevation accuracy is 1.70 m. The non-parametric and proposed methods can decrease the overall elevation error to 1.35 m and 1.34 m, respectively. Moreover, the RMSEs of the non-parametric and proposed methods in Figure 12a are roughly consistent, demonstrating that our proposed method can get the same good results as the non-parametric method regarding elevation accuracy. Figure 12b is the plane accuracy statistics diagram. Figure 13 is the error vector diagram of plane checkpoints, wherein circles are the checkpoint locations, and arrows represent the size and direction of plane errors. Before adjustment, the plane error is about 14.31 m, and our method can improve the plane accuracy to 4.14 m. Besides, it can be found that the systematic plane error of each orbit disappears after the adjustment, which shows that the proposed method can effectively improve the horizontal geolocation accuracy of the spaceborne InSAR-generated DEM.

In order to visually display the adjustment results, we compensated the adjustment parameters to the original parameters and re-produced DEMs, as shown in Figure 14b,c are the adjusted maps corresponding to Figure 11b,c, respectively. In Figure 14b, the height fracture phenomenon disappears compared with Figure 11b; that is, the abnormal slope values at the mosaic boundary have disappeared. Furthermore, the plane fracture phenomenon shown in Figure 11c also disappears, as Figure 14c. To clearly illustrate the improvement of the relative plane accuracy after adjustment, we further enlarge the critical areas of Figure 11c and Figure 14c, as shown in Figure 15. The dotted line indicates the rail’s shape. Before adjustment, the two tracks of DEM located at the mosaic area had different degrees of plane errors, resulting in the fracture of the rail, as Figure 15a. After adjustment, the rail area presents a continuously changing straight line, as Figure 15b. The disappearance of the plane fracture phenomenon demonstrates that the plane geolocation accuracies of the two DEMs tend to be consistent. Therefore, the proposed method effectively enhances the relative geolocation accuracy of InSAR-generated DEMs.

In addition, Figure 16 shows the change trends of adjustment parameter increments with iterations. After three to four iterations, the adjustment increment tends to zero, consistent with the simulated experiments, indicating that the proposed method has a strong convergence characteristic in natural environments.

5. Conclusions

This paper proposed a three-dimensional adjustment method for spaceborne bistatic InSAR based on the RDP model. We use the total differential theory to linearize the RDP model to deduce the sensitive equations and obtain the adjustment equations of available control data, including HCPs, PCPs, HTPs, and PTPs, with the help of differential geometry theory. The proposed adjustment method can correct the errors of spatial baseline, slant range, and timing parameters and jointly improve the three-dimensional accuracy of InSAR-generated DEMs. We conduct experiments using 29 scenes of TanDEM-X CoSSC data covering the Henan Province, China. The simulated data experiment shows that the proposed method can accurately calculate the systematic errors of physical parameters under different noise levels. The adjustment accuracies of the parallel baseline, slant range, and timing parameters are 0.05 mm, 0.1 m, and 0.006 ms, respectively. In the real data experiment, the proposed method improves the plane and elevation accuracies to 4.14 m and 1.34 m, respectively. It effectively overcomes the horizontal and vertical fracture phenomenon in DEM mosaic maps and improves the relative accuracy of InSAR-generated DEMs.