SAR Interferometric Baseline Reﬁnement Based on Flat-Earth Phase without a Ground Control Point

: Interferometric baseline estimation is a key procedure of interferometric synthetic aperture radar (SAR) data processing. The error of the interferometric baseline a ﬀ ects not only the removal of the ﬂat-earth phase, but also the transformation coe ﬃ cient between the topographic phase and elevation, which will a ﬀ ect the topographic phase removal for di ﬀ erential interferometric SAR (D-InSAR) and the accuracy of the ﬁnal generated digital elevation model (DEM) product for interferometric synthetic aperture (InSAR). To obtain a highly accurate interferometric baseline, this paper ﬁrstly investigates the geometry of InSAR imaging and establishes a rigorous relationship between the interferometric baseline and the ﬂat-earth phase. Then, a baseline reﬁnement method without a ground control point (GCP) is proposed, where a relevant theoretical model and resolving method are developed. Synthetic and real SAR datasets are used in the experiments, and a comparison with the conventional least-square (LS) baseline reﬁnement method is made. The results demonstrate that the proposed method exhibits an obvious improvement over the conventional LS method, with percentages of up to 51.5% in the cross-track direction. Therefore, the proposed method is e ﬀ ective and advantageous.


Introduction
During the past two decades, the interferometric synthetic aperture radar (InSAR) has been gradually matured in terms of theory and widely used in topographic mapping [1] and surface deformation monitoring [2]. For example, the shuttle radar topography mission (SRTM) data forms a popular digital elevation model (DEM, a surface model representing ground elevation as an ordered set of numerical arrays), which is acquired by airborne SAR [3]. Recently, the application fields of InSAR have been expanded from earthquakes, volcanic eruptions, ground subsidence, and other large-scale macroscopic monitoring to infrastructure monitoring, such as bridges, dams, and buildings, as well as high-accuracy 3-D urban reconstruction [4][5][6], the prospect of which is very broad.
However, the accuracy of InSAR is affected by a variety of errors, such as atmospheric delay error, DEM error, random phase error, and baseline estimation error. Among such errors, interferometric baseline estimation error systematically and seriously affects the accuracy of InSAR measurement [7]. The error of the interferometric baseline affects not only the removal of the flat-earth phase, but also the transformation coefficient between the topographic phase and elevation, which will affect the topographic phase removal for D-InSAR and the accuracy of the finally generated DEM product for InSAR. The present methods of initial baseline estimation make it difficult to obtain high-accuracy baseline estimations due to factors such as the precision of the satellite orbit state vector.
The baseline estimation method based on the satellite orbit state vector utilizes the difference in instantaneous satellite coordinates to obtain the interferometric baseline estimation, while it is directly Remote Sens. 2020, 12,233 3 of 17 iterations and then computed the baseline estimation based on LS. However, the accuracy of the iterative calculation of GCP is limited.
In recent years, some researchers have combined SAR interferometric baseline refinement with multi-temporal InSAR and proposed corresponding methods [27,28]. However, when calculating radar vectors in line of sight (LOS), the baseline refinement models they have used are highly dependent on the precision of ground elevation information, which affects the accuracy of the baseline refinement results to some extent.
In order to obtain an interferometric baseline with a high accuracy and avoid the dependence of the refinement model on ground elevation, this paper firstly investigates the geometry of InSAR imaging and analyzes the rigorous geometric relationship between the interferometric baseline and flat-earth phase in detail. Then, a baseline refinement method without GCPs is proposed and the geometry is converted from precise baseline estimation to an InSAR flat-earth phase. We develop the theoretical model and resolving method. Synthetic and real SAR datasets are used in the experiments. Compared with the conventional LS method, the proposed method is more effective and advantageous.

SAR Interferometric Baseline Model
Since the orbits of an SAR satellite do not coincide with each other at the time of imaging, we refer to the spatial vector consisting of the satellite platform orbital positions as the SAR interferometric baseline, i.e., where O(·) = [X(·), Y(·), Z(·)] is the 3-D spatial coordinate vector of the SAR sensor platform at ground point P; t M and t S are the acquisition times of point P in the master and slave image, respectively; and B x , B y , and B z are the baseline components in X, Y, and Z directions, respectively. As shown in Figure 1, the baseline can be represented as the parallel component along the LOS of the radar and the perpendicular component perpendicular to the LOS, i.e., It should be noted that the precise look angle θ is required for this baseline decomposition, while its calculation depends on the ground elevation. Here, we introduce an alternative baseline representation independent of ground elevation for subsequent theoretical derivation.
We adopt an SAR platform-fixed coordinate system, TCN (see bibliography [22,26] for details, T, Along-Track, the flight direction of satellite; C, Cross-Track, perpendicular to the plane consisting of T and N to meet the left-handed coordinate system; N, Normal-Direction, the direction from the center of the earth towards the satellite). Then, the transformation matrix from the earth-fixed coordinate to the TCN coordinate is formed by where → t , → c , and → n , represent the unit vectors of the T axis, C axis, and N axis in the earth-fixed coordinate, respectively [26]. The baseline can be converted to the TCN coordinate by Remote Sens. 2020, 12, 233 4 of 17 As we can see from Formula (5), the decomposition of the baseline in the TCN coordinate does not depend on the ground elevation, which will be adopted for the subsequent theoretical derivation.
As we can see from Formula (5), the decomposition of the baseline in the TCN coordinate does not depend on the ground elevation, which will be adopted for the subsequent theoretical derivation. Since the SAR interferometric baseline in the whole interferogram is not a fixed value, we can model the baseline components using linear formulas, i.e., where ,0 c B and ,0 n B respectively represent the components in the directions C and N of the SAR interferometric baseline model at the reference time, namely the baseline constants. We usually set the image centers as the time reference points (then  and is no longer within consideration.

Coefficient of Variation
According to SAR imaging geometry (see Figure 1), the unwrapped flat-earth phase can be expressed by O(t S ) are the coordinates of the master and slave satellites in the geocentric coordinate system, respectively; t M and t S are the respective imaging time of the master and slave images; r 1 and r 2 denote the respective slant ranges of target P in the master and slave image; and are the nominal slant ranges with respect to the ellipsoid reference surface in the master and slave acquisitions, respectively; is the tilt angle of the baseline; θ and θ 0 are the look angle and nominal look angle, respectively;B ⊥ and B are the perpendicular and parallel baseline, respectively.
Since the SAR interferometric baseline in the whole interferogram is not a fixed value, we can model the baseline components using linear formulas, i.e., where B c,0 and B n,0 respectively represent the components in the directions C and N of the SAR interferometric baseline model at the reference time, namely the baseline constants. We usually set the image centers as the time reference points (then B c,0 and B n,0 are the components of the image center baseline in the directions C and N, respectively). α c and α n are the baseline change rates in the directions C and N, respectively. t is the time at which the point was acquired relative to the reference time. B c,0 , B n,0 , α c , and α n are unknown model parameters. Here, taking the image co-registration into account, the baseline component in T has been eliminated, which means B t = 0 and is no longer within consideration.

Coefficient of Variation
According to SAR imaging geometry (see Figure 1), the unwrapped flat-earth phase can be expressed by where λ is the radar wavelength, where ∆ is the error vector of observation L. Equation (11) is a non-linear equation, and in this paper, it will be transformed into a linear one. By expanding Equation (11) into Taylor Series at , retaining zero-and first-order items, and writing it as an error equation, we can get 2,i ) are the nominal slant ranges of master and slave acquisition, respectively (i.e., with respect to the ellipsoid reference surface, see Figure 1).
Then, Equation (12) can be rewritten as a matrix form: The corresponding normal equation is where P = σ 2 0 D −1 is the weighted matrix of the observation vector L. Assuming that the observation values are mutually independent, let the prior variance matrix be where σ 2 L k represents the variance of observation value L N . Then, we can utilize the residual unwrapped differential interferometric phase φ res,unw without system signals to estimate σ 2 L k . From Equation (14), we can see that five or more observations, i.e., N ≥ 5, are needed to estimate the five unknowns. The best linear unbiased estimation of x iŝ where P is the weighted matrix. Therefore, the best linear unbiased estimation of X iŝ By analyzing design matrix A, if the initial parameter vector X 0 is close to the true vector X, the numerator of each element in the first four columns of A will become very small. After being divided by , the first four columns will become highly linearly dependent on each other, i.e., the design matrix A will tend to be multicollinear. Therefore, the normal equation, Equation (14), will become an ill-condition. The least-square estimates obtained with Equation (16) are not reliable in this situation. Therefore, we can introduce ridge estimation [29] to solve this problem and Equation (14) can be rewritten as where k is a positive constant value named the ridge parameter and I is the identity matrix. It should be noted that ridge estimation is used to improve the ill pose condition of Equation (14), but the small value k still leaves A T PA + kI singular. In this context, truncated singular value decomposition [30] (TSVD) is used to compute the pseudo-inverse of A T PA + kI. Let SVD of A T PA + kI be where U and V are 5 × 5 unitary matrixes, and the diagonal matrix S = Diag(s 1 , s 2 , s 3 , s 4 , s 5 ) contains non-negative singular values s 1 ≥ s 2 ≥ s 3 ≥ s 4 ≥ s 5 ≥ 0. The threshold of TSVD is defined with where max{s i } is used for calculating the largest singular value. The pseudo-inverse of S with the condition of S r > thres is Then, the pseudo-inverse of A T PA + kI is Therefore, the ultimate estimation of X iŝ Remote Sens. 2020, 12, 233 7 of 17 By substituting Equation (23) into Equation (17), the estimation of X can be resolved.

Model Resolving
According to the previous analysis, we can easily solve the model problem above and get the estimation X with the determination of the weight matrix P and ridge parameter k, which are not given precisely at first. Simultaneously, some other problems still exist in Equation (11). Therefore, iterative parameter estimation is adopted to solve these problems. The detailed steps are as follows: Step 1: Model initialization. Let the initial weighted matrix P = I be an identity matrix, ridge parameter k = 10 −3 , which will not cause large deviations, and first four values are initial baseline parameters calculated from the orbital state vector. Furthermore, let the initial reference point phase be φ 0 = 0; Step 2: Expand Equation (11) into Taylor Series at X 0 to compute the design matrix A and vector l to establish the normal equation Equation (14). To evaluate the quality of the model, obtain the sum of weighted squared residuals (SWST) by where Step 3: Add k to the diagonal of coefficient matrix in normal equation Equation (14), establish the ridge normal equation Equation (18), and calculate the corresponding pseudo-inverse A T PA + kI † ; Step 4: Calculatex based on Equation (23). Expand Equation (11) into Taylor Series at X 0 +x, calculate the new design matrix A and vector l , and establish the new SWST new by Equation (24); Step 5: Compare SWST new with SWST old , determine whether to accept the parameter vectorx, and output the sum of weighted squared residuals x 2 new = min{SWST old , SWST new }. If SWST new < SWST old , then calculateX with Equation (17), and update X 0 +x, as well as the design matrix A = A and vector l = l , reducing ridge parameter k = 0.1·k. Update the weighted matrix P ii with P ii = P ii / ∆L 2 i + z , where z is a positive number, such as 10 −3 , preventing the new weights from being too large, and then let SWST old = SWST new . Conversely, keepX, X 0 , A, l, and P unchanged, and update k = 0.1·k. Actually, the SWST changes slowly with the magnitude of k. Therefore, scaling the magnitude of k up or down at the ratio of ten will accelerate convergence [31][32][33]; Step 6: Determine the terminative condition of the iteration: if it has just iterated once and Step 3 to for a second iteration; if it has iterated more than once, then the terminative condition is which meets twice or the maximum iteration number is greater than 20. If the condition has been met, then output the result parameter vectorX. Alternatively, return to Step 3 and go on with the iterations. As the number of iterations increases, the ridge parameter k gradually becomes unbiased.

Algorithm Flowchart of Baseline Refinement
As we can see from Section 3.1, SAR interferometric baseline refinement can be strictly performed as long as the flat-earth phase φ f lat,unw is accurately obtained. There are two steps for obtaining an accurate flat-earth phase. First, calculate the initial baseline using the orbital state vector, as well as the corresponding flat-earth phase φ f lat,init . Due to the inaccurate baseline, φ f lat,init is not equal to φ f lat,unw . Furthermore, process differential interferometry using the initial baseline combined with external DEM and unwrap the interferogram. Then, calculate the flat-earth phase residue with planar polynomial fitting, which will be verified in the experiment part. Add φ f lat,res and φ f lat,init to make φ f lat,unw . Therefore, a detailed flow of the proposed algorithm for baseline refinement is as follows:

1.
Register master and slave images, and estimate the initial baseline using the orbital state vector to get B 0 c,0 , B 0 n,0 , α 0 c , and α 0 n ; 2.
Process differential interferometry of SAR data using initial baseline results, and calculate flat-earth phase φ f lat,res from the unwrapped interferogram using planar polynomial fitting; 3. Calculate flat-earth phase φ f lat,init using initial baseline estimation and add φ f lat,res and φ f lat,init obtained from step 2 to make a "true" unwrapped flat-earth phase φ f lat,unw ;

4.
Refine the interferometric baseline strictly and obtain an accurate estimation of the SAR interferometric baseline. Firstly, establish adjustment models of baseline parameters, including a functional model with φ f lat,unw obtained from step 3 based on the method introduced in Section 3.1 and the random model. Then, calculate the model parameters according to the algorithm in Section 3.2. By expanding Equation (11) into Taylor Series to build the normal equation, simultaneously calculate the sum of weighted squared residual SWST. The iteration mainly exists in ridge estimation and truncated singular value decomposition (TSVD). The magnitude relationship between x 2 old and x 2 new resulting from the iteration above is regarded as the terminative condition of iteration. If the condition is met, end the iteration and output the accurate SAR interferometric baseline estimation.
The specific flowchart is shown in Figure 2.

Polynomial Fitting of Interferometric Baseline Error
According to InSAR imaging geometry and the corresponding theoretical formulas, the error of the interferometric baseline affects not only the removal of the flat-earth phase, but also the transformation coefficient between the topographic phase and elevation, which will affect the topographic phase removal for D-InSAR and the accuracy of the finally generated DEM product for InSAR. In order to quantitatively analyze the effect of baseline error, especially the impact on the removal of the flat-earth phase, an experiment is carried out by a random synthetic dataset. At first, synthesize 50 sets of baselines with different lengths in the range of 50-2500 m, and then simulate the DEM is. Calculate the "true" flat-earth phase and topographic phase combined with the parameters of the ALOS PALSAR imaging system. Following this, add 100 groups of random errors with 0 means to each set of baseline parameters (the standard deviations of the added baseline errors are

Polynomial Fitting of Interferometric Baseline Error
According to InSAR imaging geometry and the corresponding theoretical formulas, the error of the interferometric baseline affects not only the removal of the flat-earth phase, but also the transformation coefficient between the topographic phase and elevation, which will affect the topographic phase removal for D-InSAR and the accuracy of the finally generated DEM product for InSAR. In order to quantitatively analyze the effect of baseline error, especially the impact on the removal of the flat-earth phase, an experiment is carried out by a random synthetic dataset. At first, synthesize 50 sets of Remote Sens. 2020, 12, 233 9 of 17 baselines with different lengths in the range of 50-2500 m, and then simulate the DEM is. Calculate the "true" flat-earth phase and topographic phase combined with the parameters of the ALOS PALSAR imaging system. Following this, add 100 groups of random errors with 0 means to each set of baseline parameters (the standard deviations of the added baseline errors are σ B c,0 = 1.3 m, σ B n,0 = 0.9 m, σ α c = 3 mm·s −1 , and σ α n = 2 mm·s −1 , respectively, and eliminate the baseline error more than two times the standard deviation), and use the baseline with errors to simulate the flat-earth phase and topographic phase. Next, calculate the difference value between the flat-earth phase simulated with baseline error and that without baseline error, and conduct the same operation for the topographic phase. Then, we get the residuals of the flat-earth phase φ f lat,res and residuals of the topographic phase φ top,res caused by the baseline error. Since we only add the baseline error into the simulation experiment, the resulting residuals of the flat-earth phase and topographic phase will only be related to the baseline estimation, independent of other errors (e.g., DEM errors), so the sum of residuals is also called the baseline error phase. Next, we can analyze the baseline error by polynomial fitting, for which we can use the quadratic polynomial model, that is, where φ f it is the baseline error phase obtained by fitting; x and y are normalized coordinates in the range and azimuth, respectively; and α i (i = 0, 1, · · · , 5) is the unknown model parameters; (see bibliography of [22] for more details). Then, we differ φ f lat,res and φ f it as φ f it , and calculate the corresponding root-mean-square errors (RMSEs). In this way, we can obtain the result of quadratic polynomial fitting in 5000 groups of experiments, and the RMSEs are shown as Figure 3.

Baseline Refinement of Simulation Data
On the basis of the experiments in Section 4.1.1, we first added the random errors with the range of [0,16] meters to the "true" DEM. Then, we simulated the DEM and the topographic phase by 5000 groups of the above-mentioned baselines with and without errors, to obtain the interferogram with baseline errors and DEM errors, as well as the "true" interferogram without any errors. Next, we differed the "true" interferogram and the one with errors to get the initial differential interferograms. To make the experimental results more reliable, we added the decorrelation error, atmospheric delay error, and random phase noise to the initial interferogram [34], so as to generate 5000 scenes of differential interferograms.
In the experiment, we refined the baselines of 5000 scenes of interferograms using the LS method. In order to reduce the effect of random phase noise, we filtered the differential interferograms by the Goldstein filter (filtering window is 32 × 32 pixels, filtering factor = 0.5, and the size of overlap is 14 pixels) [35], and then unwrapped the filtered interferograms by minimum cost flow (MCF) [13] to obtain the differential interferograms. In the experiment, we estimated the baselines accurately with 50 × 50 points selected uniformly in each interferogram. For baseline refinement based on the flat-earth phase, in this study, we firstly retrieved the flat-earth phase

Baseline Refinement of Simulation Data
On the basis of the experiments in Section 4.1.1, we first added the random errors with the range of [0,16] meters to the "true" DEM. Then, we simulated the DEM and the topographic phase by 5000 groups of the above-mentioned baselines with and without errors, to obtain the interferogram with baseline errors and DEM errors, as well as the "true" interferogram without any errors. Next, we differed the "true" interferogram and the one with errors to get the initial differential interferograms. To make the experimental results more reliable, we added the decorrelation error, atmospheric delay error, and random phase noise to the initial interferogram [34], so as to generate 5000 scenes of differential interferograms.
In the experiment, we refined the baselines of 5000 scenes of interferograms using the LS method. In order to reduce the effect of random phase noise, we filtered the differential interferograms by the Goldstein filter (filtering window is 32 × 32 pixels, filtering factor α = 0.5, and the size of overlap is 14 pixels) [35], and then unwrapped the filtered interferograms by minimum cost flow (MCF) [13] to obtain the differential interferograms. In the experiment, we estimated the baselines accurately with 50 × 50 points selected uniformly in each interferogram. For baseline refinement based on the flat-earth phase, in this study, we firstly retrieved the flat-earth phase φ f lat,res by quadratic polynomials from the unwrapped differential interferogram and calculated φ f lat,init with initial baselines, the sum of which is the "true" unwrapped flat-earth phase φ f lat,unw . Then, we refined the baselines based on the theories presented in Sections 3.1 and 3.2. Similarly, the accurate baseline estimation by the LS method is also based on the above-mentioned 50 × 50 GCPs of which the elevations are without DEM errors.

Real Dataset Experiment
On August 7, 2010, a sudden rainstorm occurred in the northeast mountains in Zhouqu County in Gansu Province, causing extraordinarily serious mountain torrents and geological disasters. The resulting debris flow was about 5 km long and 300 m wide, the flow region of which was wiped out. Two scenes of PALSAR data covering the interest area were used to verify the interferometric baseline refinement method mentioned above, and the details are provided in Table 1. Given the lack of GCPs in such an area, the conventional LS method cannot be applied here. Therefore, the new method proposed in this paper, refining the baselines without GCP, is a better choice. GAMMA software [13] was used for data preprocessing. Firstly, we transformed the Raw data into single look complex images (SLC) with the modular SAR processor (MSP) module. Then, we generated an interferogram from SLCs according to differential InSAR procedures, estimated the initial baselines based on orbit state vectors, and removed the flat-earth phase with SRTM DEM [36,37] by using the interferometric SAR processor (ISP) module. The generated interferogram is shown as Figure 4a.  GAMMA software [13] was used for data preprocessing. Firstly, we transformed the Raw data into single look complex images (SLC) with the modular SAR processor (MSP) module. Then, we generated an interferogram from SLCs according to differential InSAR procedures, estimated the initial baselines based on orbit state vectors, and removed the flat-earth phase with SRTM DEM [36,37] by using the interferometric SAR processor (ISP) module. The generated interferogram is shown as Figure 4a.

The Discussion about Polynomial Fitting
As we can see from Figure 3, if the length of the baseline is less than 2500 m, the quadratic polynomial can fit the flat-earth phase residuals well due to the baseline errors, of which the RMSEs are as small as a 10 −3 rad magnitude. With the increase of the lengths of baselines, fitting RMSEs will gradually converge to a small range. Even the RMSEs in such cases will present an increasing trend, the trend is still small and the maximum RMSE is less than 4.2 × 10 −3 rad. Therefore, it is proved that the quadratic polynomial can fit the flat-earth phase residual φ f lat,res well, where the fitting error is too small to be counted. In this way, we can accurately obtain the flat-earth phase by using the quadratic polynomial to fit its residuals with the initial baseline estimation, so that we can refine the baseline strictly based on the theories presented in Sections 3.1 and 3.2. c,0 , B 0 n,0 , α c , and α n , respectively. If the gaps between two kinds of errors are small, the scattered points in the figure will regress to the red (diagonal) line. By comparing the subgraphs (a) in Figures 5 and 6, we can see that the baseline parameter B c,0 obtained by our method regresses more than LS, which testifies that B c,0 estimated by our method is more accurate. By comparing subgraphs (b), the accuracy of B n,0 obtained by our method is a little poorer than that of LS. From the comparison between subgraphs (c) and (d) in Figures 5 and 6, we can see that our method has great advantages in α c and α n (baseline change rate), as the estimation results almost exclusively lie on the red (diagonal) line, while large dispersion appears for the conventional LS method.     From the quantitative evaluation results, we can see that parameters B c,0 and B n,0 are almost the same for the two methods. For a deeper evaluation, we calculated the differences between simulated errors and calculated errors for the two methods and statistically analyzed them in histograms (see Figures 7 and 8). The figures show that the estimation errors of B c,0 and B n,0 are more concentrated at 0, most of which are in the range of [−0.05, 0.05] m, while the errors of the LS method are more dispersed. For parameters α c and α n , the estimations are almost all concentrated at 0 for the new method, and those for the LS method are scattered, for which quite a number of the estimation errors are greater than 10 mm/s and some errors are even greater than 20 mm/s. To conclude, the new baseline estimation method is better than the conventional LS method, especially for estimation of the baseline change rate.      For a deeper quantitative evaluation, we calculated the RMSE of the differences between simulated errors and calculated errors for the two methods, respectively. They are shown in Table 2. For baseline parameter B c,0 , the accuracy of the LS method is 6.8 cm, and only 53.3% of the estimation results are less than 5 cm in 5000 groups of experiments, while the accuracy reaches 3.3 cm of the new method, which exhibits a 51.5% improvement in estimating over the LS method. In addition, 87.1% of the results are less than 5 cm, and 33.8% are higher than those of the LS method. For parameter B n,0 , it can be considered that the accuracies of the two methods are almost the same, since the differences between RMSE and the corresponding proportion are not so large. Considering baseline parameters B c,0 and B n,0 comprehensively, it can be found that the accuracy of B n,0 is higher than that of B c,0 for the LS method, while the trend for the proposed method is the opposite, but the accuracy is superior to that of the LS method. Further, the method proposed in this paper estimates baselines based on the flat-earth phase, independent of the external elevation, so the factor affected by elevation is invalid, leading to a higher accuracy. For baseline change rates α c and α n , the accuracies of the two methods are obviously different. The RMSEs are 7.4 and 5.1 mm/s, respectively, for the LS method, where the proportions of the error range less than 0.5 mm/s are only 5.4% and 8.2%. However, the RMSEs of the new method are 0.1 and 0.1 mm/s, respectively, and the proportions of the error range less than 0.5 mm/s reach 99.7% and 100%. With respect to a standard SAR image with a moderate resolution, the interval of the acquisition is generally around 15 s, for example, it is 14 s for PALSAR. Here, we select the acquisition time of the SAR scene center as the reference time of the baseline model, and the maximum relative time is 7 s. For the error of baseline change rate mentioned above, the maximum baseline error caused by the baseline change rate is only 0.7 mm for the new method, which is almost negligible. Conversely, the baseline errors can be up to 5.2 and 3.6 cm for the LS method. Furthermore, adding errors of B c,0 and B n,0 to them, we can get baseline errors of more than 0.1 m. It can be seen that the accuracy of interferometric baseline constant estimation in the C direction for the new method is obviously improved compared to the conventional LS method, the improvement of which is up to 51.5% (the ratio of reduced RMSE to the RMSE of the LS method). It proves the effectiveness of the proposed baseline estimation method.

The Discussion about Baseline Refinement
In the experiment, we selected 50 sets of baselines (the length of the baseline is in the range of 50-2500 m) with 100 experiments. Therefore, the deeper evaluation of the effect of baseline estimation due to baseline length could be analyzed from 100 repeated experiments for each set of baselines. The mean value and the standard deviation were calculated by baseline errors resulting from the repeated experiments and are shown in Figure 9. From Figure 9, the four parameters basically remain constant with the growth of baselines, and the dispersions also remain constant. This indicates that the change of baseline length does not affect the accuracy of the baseline estimations for both methods. However, the accuracy of the new method is higher than that of the LS method in the case of the same baseline length, especially for the estimation of the baseline change rate. In view of the above analysis, we can see that the new method is more effective and reliable than the LS method.

Discussion about Real Dataset Experiment
Since the area of debris flow is small (see the area within the dashed ellipse in 4d), and there is no obvious large-scale surface deformation during the interval of the two acquisitions, the fringes with longer wavelengths are caused by not completely removing the flat-earth phase caused by baseline errors, which is called the flat-earth phase residual. To obtain the flat-earth phase residual, we used MCF to unwrap the differential interferogram (Figure 4a), and used the quadratic polynomial introduced in [23] to fit the flat-earth phase residual , flat res φ . We wrapped the flat-earth phase into the range [−π, π] to display the data conveniently (see Figure 4b). Next, we calculated  figure) is well-shown. In addition, the interferogram is very flat, except for the deformation area, where there is almost no excessive fringe change rate. For a further illustration of the effect of baseline errors, the difference between Figure 4a,b is given and shown as Figure 4c. Figure 4a shows the differential interferogram with baseline errors, which will also affect

Discussion about Real Dataset Experiment
Since the area of debris flow is small (see the area within the dashed ellipse in 4d), and there is no obvious large-scale surface deformation during the interval of the two acquisitions, the fringes with longer wavelengths are caused by not completely removing the flat-earth phase caused by baseline errors, which is called the flat-earth phase residual. To obtain the flat-earth phase residual, we used MCF to unwrap the differential interferogram (Figure 4a), and used the quadratic polynomial introduced in [23] to fit the flat-earth phase residual φ f lat,res . We wrapped the flat-earth phase into the range [−π, π] to display the data conveniently (see Figure 4b). Next, we calculated φ f lat,init with the initial baseline and added it to φ f lat,res , to get the unwrapped flat-earth phase φ f lat,unw , and then refined the baseline estimation based on the theories presented in Sections 3.1 and 3.2. In order to visually show the accuracy of baseline refinement, we generated a new interferogram with refined baselines, shown as Figure 4d. From the figure, we can see that the flat-earth phase with longer wavelengths caused by baseline errors has been totally removed, and the main deformation area (the dashed ellipse in the figure) is well-shown. In addition, the interferogram is very flat, except for the deformation area, where there is almost no excessive fringe change rate. For a further illustration of the effect of baseline errors, the difference between Figure 4a,b is given and shown as Figure 4c. Figure 4a shows the differential interferogram with baseline errors, which will also affect the removal of the flat-earth phase and topographic phase. By subtracting Figure 4b from Figure 4a, Figure 4c obviously shows the phase residuals caused by baseline errors (e.g., arrows A/B pointing areas), while there are no such errors in the differential interferogram with refined baselines. This further demonstrates the new method is effective for baseline refinement.

Conclusions
Interferometric baseline estimation is a key procedure of InSAR data processing. The error of the interferometric baseline affects not only the removal of the flat-earth phase, but also the transformation coefficient between the topographic phase and elevation, which will affect the topographic phase removal for D-InSAR and the accuracy of the finally generated DEM product for InSAR. To obtain a highly accurate interferometric baseline, this paper firstly investigated the geometry of InSAR imaging, and established a rigorous relationship between the interferometric baseline and the flat-earth phase. Then, a baseline refinement method without a ground control point (GCP) was proposed, where a relevant theoretical model and resolving method were developed. Synthetic and real SAR datasets were used in the experiments, and comparisons with the conventional least-square (LS) baseline refinement method were made. The results demonstrate that the proposed method exhibits an obvious improvement over the conventional LS method, with percentages of up to 51.5% in the cross-track direction. Therefore, the proposed method is effective and advantageous.
In general, the application of the conventional LS method is limited due to the lack of available GCPs in the study area. The baseline refinement based on the LS method processes the observation (unwrapped interferometric phase) with equal weight by using the unwrapped interferometric phase and the rigorous geometric relationships of baselines when calculating the equations. Due to the errors, such as atmospheric delay error, the equal weight process will result in deviations of the results, which is obvious in the interferogram, which displays an apparent fluctuation of atmospheric error.
The interferometric baseline refinement proposed in this paper takes advantage of the unwrapped flat-earth phase and rigorous relationship of baselines, which is independent of GCP. During modeling, since the observation is the flat-earth phase without atmospheric error, the effect of atmospheric error on baseline estimation has been avoided. However, it should be noted that the flat-earth phase residual is obtained by using global polynomial fitting for the interferogram, so that the flat-earth phase may absorb some other signals, such as the deformation signal with longer wavelengths distributing globally. For this insufficiency, we can parameterize the flat-earth phase residual by using time series InSAR, where the flat-earth residual is calculated at the same time as the deformation parameters. Relevant thoughts have been studied by researchers, and will be introduced into the baseline refinement method proposed in this paper in the future to make it more general.