Bayesian Filtering Multi-Baseline Phase Unwrapping Method Based on a Two-Stage Programming Approach

Featured Application: This paper proposed a novel TSPA-based Bayesian filtering multibaseline phase unwrapping (MBPU) framework. This framework can provide a theoretical basis for future MBPU studies. The proposed framework can provide an effective theoretical basis for high-precision DEM acquisition. Abstract: Phase unwrapping (PU) has been a key step in the processing of interferometric synthetic aperture radar (InSAR) data, and its processing accuracy will directly affect the reconstruction results of digital elevation models (DEMs). The traditional single-baseline (SB) PU must be calculated under continuity assumptions. However, multi-baseline (MB) PU can get rid of the limitation of continuity assumption, so reasonable results can be obtained in regions with large gradient changes. However, the poor noise robustness of MBPU has always been a key problem. To address this issue, we transplant three Bayesian filtering methods with a two-stage programming approach (TSPA), and propose corresponding MBPU models. First, we propose a gradient-estimation method based on the first step of TSPA, and then the corresponding PU model is determined according to different Bayesian filtering. Finally, the wrapped phase can be obtained by unwrapping, one by one, using an effective quality map based on heapsort. These methods can improve the robustness of the MBPU methods. More significantly, this paper establishes a novel TSPA-based Bayesian filtering MBPU framework for the first time. This is of great significance for broadening the research of MBPU. The proposed methods experiments on simulated and real MB InSAR datasets. From the results, we can see that the TSPA-based Bayesian filtering MBPU framework can significantly improve the robustness of the MBPU method.


Introduction
Interferometric synthetic aperture radar (InSAR) has become a valuable and important method for digital elevation models (DEMs) production [1][2][3][4]. Phase unwrapping (PU), as a difficulty in the research of InSAR data processing, has received widespread attention in recent decades [5][6][7]. Traditional single-baseline PU (SBPU) obtains the unwrapped results based on a phase continuity assumption (PCA). As shown in Equation (1), the gradients of adjacent pixels are limited to  . Therefore, due to the limitation of the PCA, SBPU cannot obtain better results in regions with large gradient changes [8,9]. Compared to SBPU, the multi-baseline (MB) PU method will increase the gradient ambiguity coefficient from  [ 1,1] to  [ , ] n n , so that the MBPU method can get rid of the PCA limitation [10,11]. Therefore, the method can obtain ideal results, even in areas with terrain exhibiting phase layover and phase undersampling, which has caused MBPU to receive increasing attention in recent years [12][13][14].
where    ( , 1) k k is the phase gradient between neighboring pixels. In recent years, more and more MBPU methods have been proposed. These can be divided into two categories. One of methods requires parameter phase estimation, which mainly includes maximum likelihood (ML) estimation [15] and maximum a posterior estimation (MAP) [16]. This type of method solves MB interferograms according to the probability density function, and can obtain the final DEM without PU. However, the noise robustness of the ML method is poor, and the method is extremely susceptible to noise. A small amount of noise will cause an obvious "burr" phenomenon in the unwrapped results, and it requires a sufficient number of high-quality interferograms to obtain satisfactory noise robustness. The MAP method is obtained by establishing a Markov random field, based on the ML method. The MAP method can obtain more accurate DEM reconstruction results, but its reconstruction time is longer, and the optimization steps are more complicated. These shortcomings have limited the application of the ML and MAP methods, and improved methods have been proposed [17,18]. For example, Y. Dong [19] proposed an efficient ML estimation approach of MB SAR interferometry for refined topographic mapping in mountainous areas. The Chinese remainder theorem (CRT) is a representative of the second category of MBPU method [20]. The second category of method is different from the first type, in that it does not require parameter phase estimation. The CRT method can accurately obtain the MBPU results in the absence of noise, but noise is inevitable in real data. To solve this issue, a closed-form CRT-based algorithm was proposed [21]. The MBPU problem was transformed into a cluster-analysis problem, which reduced the shortcomings of the MB unwrapping algorithm with poor noise robustness [10,22]. Thompson [13] proposed a recursive MBPU method that analyzes the relationship between the long-baseline and the short-baseline interferograms, and then uses the short-baseline to obtain the PU results of the long-baseline. B. Jin [23] proposed an MB InSAR PU method based on a mixed integer optimization model. H. Yu [12] proposed a groundbreaking MBPU framework that can introduce an SBPU model to an MBPU method, which has important implications for MBPU research. Y. Lan [24] proposed a refined two-stage programming-based (TSPA) MB approach using a local plane model. In summary, MBPU can obtain better unwraped results, even in areas with complex terrain, but its robustness and efficiency still need further improvement. The Kalman filter is a better filtering method based on a Bayesian framework. Kim and Griffiths [25] proposed Kalman filtering to solve the problem of PU for the first time. The Kalman filter PU method has attracted wide attention because of its better noise robustness. Loffeld et al. [26] proposed an extended Kalman filter (EKF) PU method. Juan J. Martinez-Espl [27] introduced a particle filter (PF) to the PU method. W. Liu [28] proposed a cubature Kalman filter (CKF) PU method. The Bayesian filtering PU methods are not only applied in the SBPU method but also in MBPU methods. For example, Kim [25] presented an MBPU method based on 3-D Kalman filtering. X. Xie [29] proposed an extended particle filtering (EPF) MBPU algorithm based on an amended matrix pencil model and quantized path-following strategy. However, these Bayesian filtering MBPU methods start at the shortest baseline and proceed to longer and longer baselines. Strictly speaking, these methods do not completely eliminate the limitation of the PCA. Y. Gao [30] proposed a refined TSPA of PU for MB InSAR data using the unscented Kalman filter (TSPAUKF). This method uses the first step of TSPA for gradient estimation, so that TSPAUKF PU can completely get rid of the limitations of the PCA. In this paper, three MBPU methods base on TSPA with better noise robustness are proposed. This paper transplants EKF, CKF, and an unscented information filter (UIF) to the MBPU domain. More meaningfully, based on these methods, we propose a framework of Bayesian filtering MBPU methods based on TSPA. The framework of the MBPU method proposes a gradient estimation method based on TSPA, and uses median filtering method to slightly filter the estimated phase gradients. Finally, the Bayesian filtering method unwraps the pixels, one by one, using an effective quality map based on heapsort. In this paper, we use several sets of simulated data and real data experimental analysis. The results show that the TSPA-based Bayesian filtering methods have better noise robustness, and the framework provides a basis for the introduction of other Bayesian filtering methods to the MBPU method.
The rest of this paper is organized as follows. In Section 2, the MBPU theory and phase gradient estimation method are provided. In Section 3, several TSPA-based Bayesian filtering methods are introduced. In Section 4, we compare and analyze TSPA and TSPA-based Bayesian filtering methods through experiments. A summary and conclusions are provided in Section 5.

Mathematical Foundation of MBPU
InSAR technology reconstructs DEM, based on the geometric relationship of its observation. According to the characteristics of the InSAR receiver, the received interferometric phase is , so the wrapped interferometric phase must be unwrapped. The relationship between the wrapped phase and absolute phase can be expressed as: where  ( ) k , ( ) k , and ( ) n k are the absolute phase, wrapped phase, and ambiguity number of the k th pixel, respectively. From Equation (2), when the interferometric phase does not contain noise, the absolute phase  ( ) k can be obtained by finding the ambiguity number ( ) n k . According to the observation geometry, we can obtain the following formula: where ( ) h k represents the terrain height,  represents the wavelength, ( ) r k represents the slant range,  represents the incidence angle, B represents the normal baseline, and  is 2 for monostatic mode and 1 for bistatic mode. For the same study area, the terrain heights corresponding to different baseline lengths are the same. In this paper, we take the dual-baseline as an example, the equation of the baseline length is: where represents the baselines of different lengths ( i represents different interferograms), represents the ambiguity number, and  ( ) i k represents the wrapped interferometric phase. Equation (4) is the basic equation for the MBPU method. We cannot solve for two unknowns through one equation. However, because the ambiguity numbers are integers, some methods can enable us to obtain 1 ( ) n k and 2 ( ) n k using only Equation (4).

Phase Gradient Estimation
Phase gradient estimation is a key step of PU. Whether we use an SBPU or MBPU method, accurate phase gradient estimation method is necessary to obtain a high-precision wrapped phase. According to Equation (4), we can obtain similar equations of a neighboring pixel  1 k of pixel k , so we will have: According to Equation (5) and Equation (4), we can obtain as follows: We can further rewrite Equation (6) as: where     are the gradient information between neighboring pixels, and they are both integers. So, we can obtain them by the CRT-based optimization model: . According to the above equations, we can obtain the ambiguity numbers of the gradient, which is the first step of TSPA. We assume that an interferometric phase is  M N pixels. According to Equation (8), we can obtain the ambiguity numbers of these gradients. The phase gradients can be obtained based on the relationship between a point and its neighbors. Therefore, the estimated gradients can be obtained by the following formula:  represent the gradient estimation errors. As can be seen from Equation (9), the gradients estimation method proposed in this paper is not limited by the PCA.

Bayesian Filtering MBPU method
The Bayesian filtering method has been widely used in various fields. The first step of the Bayesian filtering method is to obtain the corresponding state equation and observation equation, according to the specific problem. Then, the state value and the variance of the state equation are predicted, according to the state equation. Finally, we update the data according to the observation equation. As the Bayesian filtering method involves obtaining the unwrapped phase through the state estimation method, not a simple phase gradient summation method, it has also been widely used in PU. To address the MBPU problem, we can obtain the corresponding state equation and observation equation as follows:

Multi-baseline EKFPU Algorithm
EKF further applies the Kalman filter theory to the nonlinear field, and it has been widely used in satellite positioning and signal processing. EKF expands the nonlinear function into a Taylor series and omits higher-order terms. Generally applicable to weakly nonlinear systems, the filtering performance is poor in the case of strong nonlinearity. Therefore, TSPAEKF method is more suitable for obtaining unwrapped phase with uniform noise distribution. Based on Equation (10), the proposed TSPAEKF method one-step prediction equation is: where    ( ) i k represents the predicted value of the unwrapped phase, and where ( ) i H k is the parameter of the measurement system, and  ( ) i k is the gain matrix. Equations (11) and (12) give TSPAEKF models.

Multi-baseline CKFPU Algorithm
The CKF filter uses the spherical-radial criterion to approximate the state posterior mean and covariance, and uses 2n cubature points of the same weight to obtain new state points through the system equation. Then, a prediction point of the state at the next moment is obtained. This process does not require the calculation of the complex Jacobian matrix, and it has a good filtering effect. However, the weights of the Cubature points of CKF are equal and appear symmetrically. Therefore, the accuracy of TSPACKF is lower than the other two proposed methods. TSPACKF is more suitable for fast PU with less noise. According to CKF, the cubature points and weights of TSPACKF can be expressed as: jth propagated cubature point. Therefore, we can obtain the equation as: where ] [ Chol is the Cholesky decomposition, and ( ) j v k is the observation variance. According to CKF theory, we can obtain the equation as: where   ( ) Therefore, we can obtain the unwrapped phase and its corresponding error covariance matrix as the following: Equations (10), (11), and (13)- (16) are the equations of the TSPACKF method.

Multi-baseline UIFPU Algorithm
UIF is a newly proposed nonlinear filtering method based on the UKF. UIF improves the robustness of the UKF model by adding the information matrix and its corresponding information state vector to the UKF model. UIF performs deterministic sampling based on UT transformation, and uses a series of specified samples to approximate the posterior probability density of the system state. Therefore, TSPAUIF is more suitable to obtain the phase of complex noise regions. The sigma points can be obtained by the following formula: we can obtain the one-step prediction equation as: We let ( ) k and  ( ) k be the information matrix and its corresponding information state vector, respectively, at pixel k . We can obtain ( ) k and  ( ) k by the following formula: We can obtain   ( ) i k and   ( ) i P k using the following formula: Through the above analysis, we can obtain  ( ) i k and  ( ) i P k as follows: The above three models are MBPU methods, based on the Bayesian filtering method proposed in this paper. All three methods can simultaneously perform PU and filtering. For the convenience of expression, the above equations are all expressed in a one-dimensional shape. We can replace k in the formula with ( , ) m n to represent the two-dimensional method.
We can use Figure 1 to represent the two-dimensional Bayesian filtering PU method. In the figure, the red dots represent the pixel to be unwrapped. First, we find the eight pixels adjacent to the pixel to be unwrapped and find the unwrapped phase. The green rectangle in Figure 1 represents the unwrapped pixels. Then, we obtain the error map according to the coherence coefficient, and obtain the weight values according to the error of the corresponding pixels. Finally, we weight the red dots pixel with green rectangular pixels to obtain the final unwrapped pixels. The red rectangle in Figure 1 represents the final unwrapped pixels.

Framework of TSPA-based Bayesian filtering MBPU
In summary, this paper proposes the TSPAEKF, TSPACKF, and TSPAUIF MBPU methods, which are all Bayesian filtering MBPU methods. The MBPU methods with filtering functions proposed in this paper are essentially different from methods proposed in previous research. This paper proposes MBPU methods with good noise robustness, and it is more meaningful that this paper proposes a TSPA-based Bayesian filtering MBPU method framework, whose schematic representation is shown as Figure 2. First, we propose a method based on the first step of TSPA to estimate the distance and azimuth gradients (Equation (9)). The gradients estimated by this method can eliminate the continuity assumption. Then, the equations are obtained that correspond to different Bayesian filtering models. Finally, the unwrapped phase is obtained by unwrapping one pixel at a time according to the heapsort strategy. Through this framework, some other Bayesian filtering methods can also be transplanted to MBPU methods, such as particle filtering, unscented particle filtering, and adaptive filtering. The proposed framework has great significance for broadening the research of MBPU methods.

Results and Discussion
In this section, we use two simulated InSAR datasets and a realistic dataset to verify the performance of the proposed TSPA-based Bayesian filtering MBPU methods. The second group of simulation data is used to analyze the noise robustness of different methods. To quantitatively evaluate the performance, the root mean-square error (RMSE) is used to evaluate the results (as Equation (23)): where  is the unwrapped phase collected from the different MBPU methods,  ture is the reference unwrapped phase, and  2 is the quadratic norm. The high-precision ICESat data points are used to verify the accuracy of the realistic dataset.

Experiment 1
In the first experiment, the performance of the TSPAEKF, TSPACKF, and TSPAUIF methods is verified with a simulated InSAR dataset, and compared with traditional TSPA. Figure 3a is the reference elevation used in this experiment. We can see that the simulated data has complicated terrain fluctuations. Figure 3b shows the reference unwrapped phases with a baseline length of 389.20 m, whereas Figure 3c shows the reference unwrapped phases with a baseline length of 112.10 m. Figure 3d shows the simulated wrapped phases of Figure 3b, and Figure 3e shows the simulated wrapped phases of Figure 3c. From Table 1, we can see the parameters of simulated InSAR data. Figure 3f is the PU result of Figure 3d by TSPA, while Figure 3g-i illustrate the unwrapped phase of Figure 3d obtained by the TSPAEKF, TSPACKF, and TSPAUIF methods, respectively. Figure 3jm illustrates the errors between Figure 3f-i and Figure 3b. It can be seen from the figure that TSPA produces a significant error transmission phenomenon, due to the influence of noise. Although TSPAEKF, TSPACKF, and TSPAUIF also produce error transmission phenomena, their error transmission areas are significantly smaller than that of TSPA. We can see that different methods produce errors in areas with large gradient changes. This is due to the low interference phase quality in areas with large gradient changes. Although they limit the range of error transmission to a small range, they produce obvious over-smoothing in areas with large gradient changes. Figure 3n Figure 3c. From the results of the short baseline, it can be seen that these interferograms can obtain better PU results than long baseline interferograms, but that there is still an error transmission phenomenon in the areas of dense noise. The MBPU method of this Bayesian filtering method has better accuracy. In comparison, several Bayesian filtering MBPU methods can obtain more accurate results. From Table 2, we can see that the long baseline results obtained by TSPA are 3.2803 radians, while the results of TSPAEKF, TSPACKF, and TSPAUIF are 1.6736 radians, 1.6856 radians, and 1.6530 radians, respectively. The short baseline results are similar to the long baseline results. However, the short baseline interferogram is less difficult to unwrap, so the accuracy of these PU methods is very close. Nevertheless, the accuracy of the proposed MBPU methods is still better than that of traditional TSPA. It can be seen from the results that the TSPAbased Bayesian filtering MBPU methods proposed in this paper have a better noise robustness than TSPA. In this experiment, the TSPA consumes 40.55 s, while the TSPAEKF, TSPACKF, and TSPAUIF consume 13.84 s, 11.67 s, and 8.38 s, respectively. The efficiency of proposed TSPA-based Bayesian filtering MBPU methods are better than TSPA.     Figure 3(j-m,r-u) (unit: radians).

Experiment 2
In this experiment, the performance of TSPA and TSPA-based Bayesian filtering MBPU methods was verified using the simulation data in Figure 4d,e. Figure 4a is the employed digital elevation model in the experiment. The data are obtained by MATLAB simulation. Figure 4b,c show the reference unwrapped phases with different normal baseline lengths. Figure 4d shows the simulated wrapped phases of Figure 4b, with added phase noise using 0.95 mean correlation coefficient, and Figure 4e shows the simulated wrapped phases of Figure 4c, with added phase noise using 0.9 mean correlation coefficient. Figure 4f illustrates the unwrapped phase of Figure 4d  However, due to the influence of noise, an obvious PU error has been generated in the upper right corner. From Figure 4g-I,k-m, it can be seen that the proposed TSPA-based Bayesian filtering MBPU methods can obtain better results than traditional TSPA, and obviously have better noise robustness. However, there is still a small amount of over-smoothing in the gradient-change regions. Figure 4n illustrates the unwrapped phase of Figure Figure 4c. It can be seen from the results that, although the short baseline interferogram is easier to unwrap, the result of the traditional TSPA is still not ideal, due to the influence of noise. However, the TSPA-based Bayesian filtering MBPU methods not only obtain great PU results, but have better noise robustness. As can be seen from Table 3, the long and short baseline PU accuracy of traditional TSPA are 2.5464 radians and 1.7760 radians, respectively. Due to the influence of noise, the accuracy of traditional TSPA is not very good. The TSPA-based Bayesian filtering MBPU methods proposed in this study, especially the TSPAEKF and TSPAUIF methods, can obtain better results. The PU results of TSPACKF are relatively ineffective, but have better noise robustness. The efficiency of proposed TSPA-based Bayesian filtering MBPU methods are 7.97 s, 7.13 s, and 5.28 s, respectively. The efficiency of TSPA is 8.23 s, and it can be seen from the experiment that the proposed methods not only have better robustness, but also have better efficiency.

Experiment 3
In this experiment, we validate the proposed TSPA-based Bayesian filtering MBPU method with a real MB InSAR datasets. Figure 5a shows the real surface conditions of the study area. From the figure, we can see that the study area contains a large number of high mountains. In Table 4, we can see the major parameters of the real dataset. Figure 5b is an interferogram with a baseline length of −370.46 m, and Figure 5c is an interferogram with a baseline length of −127.79 m. From Figure 5b,c, the white lines are composed of the ICESat points. The accuracy of these ICESat points is about 14cm. It can be seen from Figure 5a that the research area has obvious mountain undulations, and it was proved [30] that the SBPU method cannot obtain high-precision PU results. Figure 5d-g illustrates the unwrapped phase of Figure 5b obtained by TSPA, and the proposed TSPA-based Bayesian filtering MBPU methods, respectively. Figure 5h-k illustrates the unwrapped phase of Figure 5c obtained by TSPA, and the proposed TSPA-based Bayesian filtering MBPU method, respectively. It can be seen at a glance that the results obtained by different MBPU methods are very similar. This is because the experimental data have less noise and a larger study area. To further compare the performance of different MBPU methods, we provide the DEMs generated from the data of different MBPU results, and use ICESat data points to verify the accuracy of the results of different methods. As can be seen from Table 5, the standard deviations obtained by traditional TSPA method are 3.25 m for Figure 5d  .37 s, respectively. They are much faster than the efficiency of the TSPA, which required 26,156.58 s. It can be seen from the experimental results that the TSPA-based Bayesian filtering MBPU methods proposed in this paper can obtain better results than traditional TSPA methods, and their efficiency is better.

Experiment 4
To verify the noise robustness of the proposed TSPA-based Bayesian filtering MBPU methods, this study uses Figure 4b,c as a reference unwrapped phase, to add phase noise with different mean correlation coefficients to the data, and we obtained the wrapped phases. As shown in Figure 6, we added phase noise to the simulated data with phase correlation coefficients of 1, 0.99, 0.98, 0.96, 0.95, 0.93, 0.90, 0.88, 0.85, 0.88, 0.75, and 0.7. It can be seen from the results that when the wrapped phases are without noise, different MBPU methods can accurately obtain the PU results. However, with the decrease of the phase correlation coefficient, the results of traditional TSPA become significantly worse, especially when the correlation coefficient is less than 0.90, as the unwrap phase produces a serious unwrap error phenomenon. However, the proposed TSPA-based Bayesian filtering MBPU methods can still obtain results when the phase correlation coefficients are 0.8, which further validates that the methods are better than the traditional TSPA. However, there are also differences among the TSPA-based Bayesian filtering MBPU methods. For example, when the phase correlation coefficient is 0.93, the long baseline result of TSPACKF produces a phenomenon of unwrap error in the same areas, but this does not occur in the other two methods. This may be caused by different methods of variance estimation. This provides a research point for the more accurate TSPA-based Bayesian filtering MBPU methods proposed in the future.

Conclusions
Three TSPA-based Bayesian filtering MBPU methods have been proposed. These proposed methods transplant the EKF, CKF, and UIF to the MBPU methods, respectively. More meaningfully, this paper proposed a TSPA-based Bayesian filtering MBPU framework. We proposed a gradient estimation method based on the first step of TSPA, which can eliminate the limitation of phase gradients by the PCA. Then, the corresponding PU model is determined according to different Bayesian filtering methods. Finally, the TSPA-based Bayesian filtering MBPU method unwraps the pixels, one by one, using an effective quality map based on heapsort. This proposed MBPU framework provides a basis for the research of MBPU methods with filtering functions.
From the results of the simulated data, the TSPA-based Bayesian filtering MBPU methods not only ignore the PCA, but also have good noise robustness. It is further seen through the robustness analysis experiments that TSPA is easily affected by noise. However, TSPA-based Bayesian filtering MBPU methods have great noise robustness. Experiments with a set of TanDEM-X data with topography exhibit severe undulations. In addition, ICESat data points are used to evaluate the accuracy of the DEMs generated by the different methods. The results of the experiments demonstrate that the proposed TSPA-based Bayesian filtering MBPU methods can obtain better results, and are more efficient than the traditional TSPA.