The solution of the ambiguity number can be divided into basic model and auxiliary model according to the relationship of between a 1 and   a 2 . Figure 4A shows the basic model with positive intercept of constraint condition and horizontal axis when a 2 > a 1 . Figure 4B shows the auxiliary model with positive intercept of constraint condition and vertical axis when a 2 < a 1 .

The equation that k 1 as the independent variable and k 2 as the dependent variable can be set up, and the constraint condition is: k 2 = k 1 m 1 m 2 + a 1 − a 2 m 2 . (8)

The optimal solution of the basic and auxiliary models can be solved, respectively: b a s i c   m o d e l { k 1 ∗ = a 2 − a 1 m 1 , k 2 ∗ = 0 ,          a u x i l i a r y   m o d e l { k 2 ∗ = a 1 − a 2 m 2 , k 1 ∗ = 0. (9)

According to the intersection of the constraint condition and the k i axis, we branch the ambiguity number variable. For the basic and auxiliary models, k i or k 2 component is branched, respectively, denoted as L P 1 and L P 2 models.   b a s i c   m o d e l { ( L P 1 ) { ( L P ) k 1 ≤ [ a 2 − a 1 m 1 ] ( L P 2 ) { ( L P ) k 1 ≥ [ a 2 − a 1 m 1 ] + 1        a u x i l i a r y   m o d e l { ( L P 1 ) { ( L P ) k 2 ≤ [ a 1 − a 2 m 2 ] ( L P 2 ) { ( L P ) k 2 ≥ [ a 1 − a 2 m 2 ] + 1                                                    (10) where [] is the rounding function. As shown in Figure 4, whether it is the basic model or the auxiliary model, if only the traversal of the ambiguity number to positive direction of the coordinate axis is considered, only the L P 2 models need to be considered.

According to Eq. 5, the objective function expression can be changed: k 2 = − k 1 m 1 m 2 + Y m 2 . (11)

The initial values of the ambiguity numbers k 1 and k 2 for the two models can be set: b a s i c   m o d e l { k 1 = [ a 2 − a 1 m 1 ] + 1 , k 2 = 0 ,             a u x i l i a r y   m o d e l { k 2 = [ a 1 − a 2 m 2 ] + 1 , k 1 = 0. (12)

The traversal path is determined by the slope. If m 1 / m 2 > 1 , then k 1 is traversed. If m 1 / m 2 < 1 , then k 2 is traversed.

In case of m 1 / m 2 > 1 , as shown in Figure 5, k 1 is incremented by one unit length in turn, and then it is determined whether k 2 satisfies an integer according to Eq. 13, until the integer solution corresponding to the minimum objective function value under this following condition is found: k 2 = k 1 m 1 m 2 + a 1 − a 2 m 2 . (13)

For the basic model, because of m 1 / m 2 > 1 , the initial value has the following restrictions: k 1 is traversed sequentially from 1. For the auxiliary model, k 1 is traversed directly from 0. 0 < a 2 − a 1 m 1 < a 2 − a 1 m 2 < 1. (14)

In case of m 1 / m 2 < 1 , as shown in Figure 6, k 2 is incremented by one unit length in turn, and then it is determined whether k 1 satisfies an integer according to Eq. 15, until the integer solution corresponding to the minimum objective function value under this following condition is found. k 1 = k 2 m 2 m 1 + a 2 − a 1 m 1 . (15)

For the basic model, k 2 is traversed directly from 0. For the auxiliary model, because of m 1 / m 2 < 1 , the initial value has the following restrictions: k 2 is traversed sequentially from 1. 0 < a 1 − a 2 m 2 < a 1 − a 2 m 1 < 1. (16)

So far, the positive traversal of determining the ambiguity number k i based on the relationship between a 1 and   a 2 has ended. The negative traversal can solve the corresponding negative ambiguity number according to the aforementioned theories.

In case of a 1 =   a 2 , as shown in Figure 7, the intercept between the constraint condition and horizontal axis becomes 0, and the optimal solution can start from the point (0,0). In this situation, the ambiguity numbers have all satisfied the integer condition, so there is no need to traverse. The corresponding optimal integer solution are as follows: k 1 = k 2 = 0. (17)

It can be seen from the aforementioned analysis that after transforming the dual-baseline InSAR phase unwrapping problem into the PIP problem, we only need to find the corresponding algorithm to solve the integer K, and there is no requirement for the adjacent interferometric phase difference to be less than half a period. Therefore, the proposed method can extend the non-ambiguity interval to [ − m π , m π ) and also has a better unwrapping ability in the phase under-sampling region. In addition, the ratio of m 1 and m 2 determines the slope of the straight line, and the pure integer programming phase unwrapping method also weakens the baseline requirement of interferometric pairs. As long as the two baselines have different lengths, the unwrapping phase can be effectively solved.

In order to further improve unwrapping efficiency, the axis symmetry theory is also introduced in this proposed method. As shown in Figure 8, if a 1 and a 2 and m 1 and m 2 in the optimal solution parameters are interchanged, the auxiliary model can be transformed into the basic model. Just as that the black line represents the basic model, and the blue line represents the auxiliary model; the two models are symmetric about the line k 1 = k 2 . Therefore, in this case, the corresponding ambiguity number can be solved by repeating the aforementioned steps. The axis symmetry theory can optimize the four traversing approaches to be considered into two traversing approaches (Figure 9), which reduces the storage space of traversing variables and improves operating efficiency.

In order to verify the feasibility of the branch and bound PIP–PU algorithm for dual-baseline InSAR, the first experiment is performed on the simulated data using the sinc function. Figures 10A,B show the three-dimensional (3D) and two-dimensional (2D) digital elevation model (DEM), respectively. Figures 10C,D, respectively, show the interferograms of the short and long baselines.

After unwrapping by the branch and bound PIP–PU method, the 3D map of the solved ambiguity numbers, unwrapping results, and phase errors are shown in Figure 11. The solved ambiguity numbers are similar to the integer layer of 3D DEM. The longer the baseline, the denser the solved ambiguity numbers, and the closer the 3D map is to the 3D DEM. This is because the baseline length is proportional to the absolute phase, and the larger the absolute phase, the larger the solved ambiguity numbers. The unwrapping results show that the overall phases have good continuity; the phase errors are not only concentrated around 0, but also the order of magnitude can reach 10^–15, which prove that this proposed algorithm based on PIP mathematical ideas can effectively guarantee the accuracy of unwrapping.

In the second experiment, we select Isolation Peak National Park data to verify the effectiveness and noise robustness performance of this proposed method. Figures 12A,B show 2D and 3D DEM. Figures 12C,D show simulated noisy interferograms with 105 m baseline and 189 m baseline, respectively, in which we added a random phase noise with a mean value of 0 and a variance of 0.1 rad². The main parameters are shown in Table 1.

There exists phase under-sampling area in the red frame of the aforementioned interferograms. The unwrapping results of 105 m and 189 m baseline interferograms using the branch and bound PIP–PU method are shown in Figures 13A,C, respectively, which are in good agreement with Figure 12A. By further making the difference between the unwrapping result and the original reference phase, the phase error distribution map is as shown in Figures 13B,D. Better unwrapping results can be obtained even in phase under-sampling areas and abrupt topographic change areas. Compared with Figures 11C,F, although the phase error is distributed around 0, the order of magnitude becomes significantly large. It shows that this proposed algorithm is relatively sensitive to the noise. The next step should be to improve the noise robustness ability of this algorithm. But even so, the unwrapping accuracy of the proposed algorithm is still optimal, which will be proved in the following quantitative analysis.

To further verify the effectiveness of the proposed method, the branch-cut method, quality-guided method, LS method, and MCF method are used to unwrap the 105 m baseline interferogram. Figures 13E,F show the unwrapping results of the branch-cut method. There exist obvious phase errors in the phase under-sampling area of the red frame, and the unwrapping effect is also poor in the abrupt topographic change area. Figures 13G,H show the unwrapping results of the quality-guided method. There is no excessive deviation, but the unwrapping effect is also poor in the phase under-sampling area and the abrupt topographic change area. Figures 13I,J show the results of the LS method. From the overall point of view, the phase continuity is guaranteed. However, the unwrapping result has too much phase deviation. Most of the data are almost 10 rad different from the original phase, and some even reach 20 rad. The main reason is that the discrete phase gradient estimation cannot reflect the true phase gradient, which introduces the large error. Figures 13K,L show the results of the MCF method. The number of error points is greatly reduced, especially in the left area. The PU effect has been greatly improved but is still not ideal in phase under-sampling areas and abrupt topographic change areas.

Taking into account the randomness of the noise generated by the simulation experiment, five groups of repeated experiments were performed by adding the same level of noise to the interferograms. Figures 14, 15 show the scatter plots of the mean and standard deviation of phase errors. Different algorithms focus around a different number. Both the error indicators of this proposed algorithm are always the smallest, followed by the MCF method, quality-guided method, branch-cut method, and LS method. Table 2 gives the average values of the five groups of repeated experiments. The branch-cut method has a large phase error mean and standard deviation, so the unwrapping effect is poor, but the unwrapping efficiency is higher. Compared with it, the quality-guided method has a lower mean and standard deviation, but the unwrapping efficiency is the worst. The LS method has the highest unwrapping efficiency, but the mean and standard deviation are the largest, and the unwrapping effect is also the worst. The error indicators of the MCF method are obviously smaller, and the unwrapping effect is better, but the diversity of network planning results in its low efficiency. For the PIP–PU method, the mean and standard deviation of phase errors are greatly reduced, and the unwrapping efficiency is sub-optimal. The reason is that the PIP–PU method uses rigorous mathematical expressions to accurately solve the ambiguity number, which is equivalent to realize the point-by-point unwrapping of the four neighborhoods of each reference point. Even in the phase abrupt change area, the ambiguity number can also be accurately solved, so the unwrapping accuracy is the highest. In addition, the traversal method adopts four-neighborhood expansion, which is similar to the branch-cut traversal method, so the time cost is not much different from the branch-cut method. Therefore, according to the overall performance evaluation, the unwrapping effect is relatively better, which verifies the effectiveness of the proposed method.

In the third experiment, we take the real Etna ERS-1/2 interferogram data (5,186 × 1998 pixels) to evaluate the practicality of this proposed method, which are shown in Figure 16.

Figures 17A,C show the 3D maps of the solved ambiguity numbers corresponding to different baseline interferograms. It can be seen that the terrain has an extreme value of the ambiguity number in the left area, so there may exist a towering terrain. Compared with the unwrapping results in Figures 17B,D, the phase distribution regularity is also consistent with that of the solved ambiguity number. Figure 17E shows the re-wrapped result of the long-baseline unwrapping phase, which is highly consistent with the original interferogram (Figure 16B), then the difference between the two can get the phase error distribution map as shown in Figure 17F, and the errors are almost 0. Finally, the inverted DEM using the long-baseline unwrapping result is shown in Figure 17G. Compared with the 3D maps of the solved ambiguity numbers (Figures 17A,C), the distribution regularity of the three are also highly consistent, which is sufficient to prove the effectiveness and practicality of the proposed algorithm.

The representative algorithms were selected from the path-following type, the minimum norm type, and the optimal estimation type for this long-baseline unwrapping experiment. The unwrapping results of the branch-cut method, LS method, and MCF method are shown in Figures 18A–C, respectively. Since the integral path fails to bypass the residual point, the error propagation leads to a large number of wrong unwrapped points in the whole row of unwrapped results of the branch-cut method. The results of the LS method are similar to the unwrapping effect of the Isolation Peak National Park in the second experiment. Although the phase continuity is guaranteed, the phase errors are relatively large. The results of the MCF method are particularly good visually.

Figures 18D–F show the re-wrapped results of long-baseline unwrapping phase, and Figures 18G–I show the phase error between re-wrapped results and original interferogram. Comparing Figures 17E,F, the re-wrapped results of the branch-cut method still show a large number of wrong unwrapped points in the whole row, and phase errors from the original interferogram exhibit a ring-shaped distribution. The re-wrapped results of the LS method appear as sparser fringes, missing details of the original interferogram, so phase errors are so large that they also exhibit a more obvious ring-shaped distribution. The re-wrapped results and phase error of the MCF method are highly consistent with the original interferogram (Figure 16B), and there are some differences only in the right area. Therefore, from the visual effect of unwrapping results for real data, the proposed algorithm in this article also has the highest accuracy.