Abstract
Synthetic aperture radar interferometry can obtain information on the surface of the earth. In fact, the accuracy of phase unwrapping directly affects the accuracy of the digital elevation model, but the problem that has been puzzling people is the low accuracy of unwrapping in high-noise regions. To solve this problem, we propose a new approach based on an improved branch-cut method, which classifies the branch cut. First, connect the short branch cuts and then use the network-flow method to connect the long branch cuts. By categorized branch cuts, each type of branch cut processes part of the residues and divides the region with high residue density into two regions with low residue density for unwrapped, which solves the problem of low phase unwrapping accuracy in high-noise regions. The experimental results of simulation and measured data show that RCBC method shows high unwrapping accuracy in high-noise regions, very short total length of connected branch cuts, and strong robustness.
1. Introduction
Every year, the casualties and economic losses caused by natural disasters such as earthquakes and landslides are huge. Therefore, the requirements for early warning of natural disasters and emergency responses after disasters are more stringent. Synthetic aperture radar interferometry (InSAR) [1–3] is a powerful and well-established remote-sensing technique used to measure many important geophysical parameters, e.g., surface height of topography [4] or ground deformation [5–8], which plays an important role in early warning of natural disasters and emergency responses after disasters. Phase unwrapping (PU) is one of the key steps of InSAR [9]. People get the wrapped phase with a value range of [10], and the information that people require is the true phase. To obtain the true phase, it is necessary to calculate the ambiguity number first, and then add the wrapped phase. The PU can be expressed by where is the true phase of the sth pixel, is the wrapped phase of the sth pixel, and is the integer ambiguity number of the sth pixel.
Equation (1) shows that the main difficulty of PU is that it is an ill-posed inverse problem, and there are two unknowns and in Equation (1) [11]. So, the PU algorithm depends on the assumptions of consistency and phase continuity. However, phase continuity is destroyed by noise and other factors, so it is necessary to eliminate the influence of discontinuity on PU. According to the method of elimination, the PU algorithm is divided into three main directions: path-following method, minimum-norm method, and other methods. The branch-cut (BC) method [12] is the representative path-following method, and the minimum-cost flow (MCF) method [13] is the representative of the minimum-norm method. The BC method reduces the effect of discontinuity on PU by residues and branch cuts. The BC method considers that the total length of the connected branch cuts should be the shortest, and the total polarity of the branch cuts should be zero. The quality-guided (QG) method [14] and other methods [15, 16] has been proposed successively. The MCF method calculates the minimum error between the unwrapped phase gradient and the wrapped phase gradient to reduce the interference of inconsistency and improves algorithms such as the statistical-cost network-flow method [17] and the minimum discontinuity method [18]. Other methods include the filtering algorithm led by the Kalman Filtering methods [19–21], the multibaseline PU methods [22, 23], and the PU algorithm such as the deep learning technique [24–27]. The PU algorithms have lower unwrapping accuracy in high-noise regions. In addition, the efficiency requirements of the PU algorithms for larger-scale InSAR data will be correspondingly higher. Therefore, PU has broad prospects for development.
The structure of the paper is as follows: The second part introduces the algorithm backgrounds of the BC method and the MCF method. The third part describes the principle of connecting branch cuts through the MCF method and the process and steps of our proposed branch-cut method based on residue classification (RCBC). The fourth part verifies the performance of the RCBC method through simulated data and measured data. The fifth part is the conclusion.
2. Backgrounds of the BC Method and the MCF Method
The BC method [12] is to convert the pixel phase discontinuity into the residues with polarity and connect the positive and negative residues with branch cuts, so that the polarity of the branch cuts is zero, which is used to mark the path that the integration cannot pass; then, obtain the unwrapped phase through flood-fill integration. The MCF method [13] is to calculate the phase deviations in the horizontal and vertical directions and then select the horizontal or vertical integration. The study finds that the MCF method and the BC method had certain commonalities, so we propose a new PU algorithm, establish the objective function of the shortest path based on the MCF method, and then use the optimization algorithm to connect the long branch cuts. The following are the backgrounds of the BC method and the MCF method.
2.1. Background of the BC Method
The BC method [12] assumes that the wrapped phase gradient is equal to the unwrapped phase gradient, and the unwrapped phase is obtained by the wrapped phase gradient for integration, and then adds the reference phase of the starting pixe. However, the discontinuity will cause noise transmission, so it is necessary to connect the branch cuts between the residues to reduce the influence of the discontinuity on the PU. The integration path including residues has noise transmission, but when the path polarity sum is zero, the integration has nothing to do with the path, and the error is only transmitted within the branch cuts. The branch cuts contain nonresidues, which is also a kind of noise transmission. Therefore, the shortest branch cuts contain all residues when the path polarity is zero, and there are as few nonresidues as possible in the branch cuts. The objective function of the shortest path for the image can be expressed by where . The values of and are 0 or 1. When it is 0, it means that there is no branch cut here, and when it is 1, it means that there is branch cut here. The algorithm steps of the BC method are as follows: (1)Calculate the residues by the wrapped phase and mark it with ±1(2)Create a window with the searched residues as the center and then connect the positive and negative residues in the window. When the sum of the polarity of the path is not zero or there is no residue in the window, expand the radius of the window and continue to connect the residues in the window(3)Connect the residues to the boundary when the preset maximum radius is exceeded. Connect all the residues as described above(4)The flood-fill integral is used to obtain the unwrapped phase(5)The unwrapped phase is obtained by the above steps, and then the estimated phase of the study area can be obtained by adding it to the reference phase of the starting pixel
2.2. Background of the MCF Method
The MCF method [13] assumes that there is a difference between the wrapped phase gradient and unwrapped phase gradient. The difference comes from phase-jump errors. In order to reduce the error of PU, the sum of the difference is required to be the smallest. The MCF method solves the problem of oil transportation by constructing a minimum objective function and uses an optimization algorithm to obtain path planning at the minimum cost. The MCF method uses this principle to transform the PU problem of minimum gradient error into the target function of minimum expenses.
The wrapped phase gradient and the unwrapped phase gradient are obtained by the Equations (3) and (4):
where and are the wrapped phase gradients in the vertical and horizontal directions, respectively. Where and are the unwrapped phase gradients in the vertical and horizontal directions.
There is a phase difference between the wrapped phase gradient and the unwrapped phase gradient, so the unwrapped phase can be obtained by
In order to obtain the smallest error between the wrapped phase gradient and the unwrapped phase gradient, calculate and, as shown in the formulas (6) and (7):
The constraints are as follows: whereandrepresent the reliability of the pixel, and their values are between [0,1]. Where and are the phase deviations between the unwrapped phase gradient and the wrapped phase gradient. The algorithm of the MCF is as follows: (1)Establish the minimum objective function(2)Solve the value of through the optimization algorithm(3)Add to the wrapped phase gradient at the phase jump and then select vertical or horizontal accumulation to obtain the unwrapped phase
The BC method has many advantages, such as simple structure, fast speed, and high accuracy in low-noise regions. However, the unwrapping accuracy is low in high-noise regions. The MCF method has the advantages of high unwrapping accuracy and strong practicability under low-noise regions, but it has the characteristics of low unwrapping accuracy and poor robustness in high-noise regions.
3. Improved the Branch-Cut Method
Based on extensive experiments, we found that the distribution of residue shows certain characteristics. There are more positive and negative residues located in four adjacent and diagonally adjacent positions, so the branch cuts are divided into two categories: short branch-cut and long branch-cut. The short branch cuts are the branch cuts that connect all the positive and negative residue pairs between the four adjacent and diagonally adjacent pixels, and the remaining residues are connected by the long branch-cuts, which is connected by the MCF method. The classification of branch cuts greatly reduces the density of residues that need to be processed when connecting each type of branch cut. The short branch cuts consume a large number of residues, and then the MCF method is used to connect the long branch cuts under the condition of low-noise density; the total length of the branch cuts is very short, so the method has high unwrapping accuracy in the high-noise regions.
3.1. The Principle of Short Branch-Cut Connection
In order to avoid branch cuts crossing each other cuts, the distribution of residues in the short branch cuts only includes the situation shown in Figure 1: four-adjacent distributions and diagonal-adjacent distributions.
(a)
(b)
In the two windows in Figure 1, the gray position is the central residue. Figure 1(a) is the position of the residues of four-adjacent distribution. The four-adjacent distribution means that one of the four positions 1, 2, 3, and 4 in Figure 1(a) has residue with opposite polarity to the gray area. Figure 1(b) is a diagram of the position of the residues of the diagonal-adjacent distribution. The diagonal-adjacent distribution means that one of the four positions 1, 2, 3, and 4 in Figure 1(b) has residue with opposite polarity to the gray area. The short branch cuts are classified twice: type a subbranch cuts and type b subbranch cuts. The type a subbranch cuts are formed by branch cuts that connect the positive and negative residue pairs shown in Figure 1(a). The type b subbranch cuts are formed by branch cuts that connect the positive and negative residue pairs shown in Figure 1(b). Here are the steps to establish the short branch cuts: (1)Search for residues in order from top to bottom and from left to right and create a window as shown in Figure 1(a) with the searched residues as the center(2)Search whether there is a residue with the opposite polarity to the central residue in the order of 1, 2, 3, and 4. If there is a residue with the opposite polarity, connect the central residue with the first found residue and then return to step 1(3)If there is no residue with opposite polarity, ignore this residue and return to step 1 until all the short branch cuts are connected
The short branch cuts are established first, and then, the long branch cuts are established through the MCF method. The process of the MCF method is as follows: establish a network model, construct an objective function, and solve the optimization algorithm. The following will introduce the principle of establishing long branch cuts by the MCF method.
3.2. The Principle of Long Branch-Cut Connection
Establish the network model, as shown in Figure 2. The area surrounded by the four bidirectional arrows in the figure is the wrapped phase pixel, the value represents the wrapped phase. The circle around each pixel indicates the polarity of the surrounding four pixels, of which ±1 indicates positive and negative residues, respectively, and the remaining blank places indicate no polarity. The nodes are connected by arrows, and the different nodes make the arrow have flows of different sizes and directions. When the flow on the arrow is not equal to zero, the arrows are the branch cuts, and the change in the polarity of the branch cuts is consistent with the change of the flow, and the direction of the flow is the connection direction of the branch cut. The network model closely connects the wrapped phases, the residues and the branch cuts through the magnitude and direction of the arrows, and the upward flow of the arrows.
The MCF method solves the multistart and multiend freight minimum problem. When the unit-price of freight is 1, the length of the path is equal to the freight. Let in formula (6), and the value of is 0 or 1, and the formula (3) is incorporated into the formula (7); then, the formulas (8) and (9) are obtained:
The constraints are as follows:
The right side of the constraint equation is the method of residue calculation, and the variable of the left side of the equation indicates that the branch cut is established in one of the four-adjacent directions of the residue. Comparing Equations (2), (6), (7), (8), and (9), it is found that the objective function has changed from the smallest error to the shortest path, so the MCF method can be used to connect branch cuts. After the branch cuts are classified, only part of the residues is processed when connecting each type of branch cuts. That is, the density of residues is greatly reduced, which is equivalent to dividing a wrapped phase image with high residue density into two low-residue density images and then separately establish branch cuts and superimpose the branch cuts. In addition, the total length of the short branch cuts is very short, and the total length of the long branch cuts connected by the MCF method is also very short, so the large errors are limited to a few pixels. Due to the special structure of the short branch cuts and the weight constraint when the MCF method connects the long branch cuts, the branch cuts will not cross other branch cuts when the branch cuts are superimposed. Therefore, the algorithm has higher PU accuracy in high-noise regions. After the branch cut connection is completed, the unwrapped phase of most areas is solved by integration, and then, the following method is used to solve the unwrapped phase of the pixels on the branch cut.
3.3. Unwrapped Phase of Pixels on Branch Cuts
There are errors in the pixels on the branch cuts, so a special method is needed to solve the unwrapped phase of these pixels. The position distribution of branch cut and unwrapped phase is shown in Figure 3.
The gray areas in Figure 3 are the pixels of the branch cut, and the others are the unwrapped phase pixels. The position distribution between branch cut and the unwrapped phase pixels basically conform to the four situations in Figure 3, and others are the combination or deformation of the four situations. “Islands” generally cannot be used as seed pixels when other pixels are PU (islands are the pixels that are not on a branch cut but are surrounded by branch cuts). The numbers in Figure 3, respectively, represent the pixels with corresponding number of branch cuts in the four-adjacent positions of the number pixel. Because on the assumption of the continuity of PU, the unwrapped phase is only related to the unwrapped phase of four-adjacent positions, so the unwrapped phase at the pixels on the branch cut is solved according to the following steps: search for branch cuts in order from top to bottom and from left to right. Search for branch cuts, center on the pixel of the searched branch cuts, and detect whether there are branch cuts at four-adjacent positions. If the position is shown as 1 in Figure 3, the average value of the sum of the unwrapped phases of the three pixels unwrapped phases of a, b, and h. If the position is shown as 2 in Figure 3, the unwrapped phase is equal to the average value of the sum of the two unwrapped phases at d and e. If the position is shown as 3 in Figure 3, the unwrapped phase value is equal to the unwrapped phase at position b plus the wrapped phase gradient between 3 and b. If the position is shown as 4 in Figure 3, the unwrapped phase is equal to the average value of the sum of the unwrapped phases of the four branch-cut pixels at the nearest position. And the unwrapped phases of all the branch-cut pixels are solved in the above manner.
3.4. Postfiltering
Since the method is suitable for PU in high-noise regions, postfiltering is added to improve the robustness of the method. The RCBC method did not use the window function of mean filtering as shown in Figure 4(a) but designed the window function shown in Figure 4(b). When the window coefficients shown in Figure 4(a) are used, the window is excessively smooth, resulting in unwrapped phase distortion. Combined the characteristics of the mean filter function and the Gaussian filter function, the RCBC method designed a new filter function.
(a)
(b)
3.5. Algorithm Steps
Figure 5 is the flowchart of RCBC PU. (1)Calculate the residue by wrapped phase(2)Preferentially connect all the a-type subbranch cuts(3)Then, connect all the b-type subbranch cuts between the remained residues(4)Then, the objective function of the shortest path is established, and the remaining residues are input into the network model(5)Use the optimization algorithm to solve the objective function, and the position where the value is 1 is recorded as the position of the long branch cuts(6)Bypassing the branch cuts to perform flood-fill integration(7)Solve the unwrapped phase of the branch-cut pixels(8)Solve the unwrapped phase of the islands(9)Postfilter processing(10)Obtain the unwrapped phase of the study area through the above nine steps
4. Algorithm Performance Analysis
The experiments are all operated on a computer with i5-1035G, 1GHZ, and 16G RAM configuration. To show the performance of the RCBC method, this method will be compared with the BC method [12], the QG method [14], the iterative least-square (ILS) method [28], the MCF method [13], the unscented Kalman filtering (UKF) method [19], and the iterated unscented Kalman filtering (IUKF) algorithm [20].
4.1. Simulation Data
Algorithms used the peak function with a dimension of and a signal-to-noise ratio (SNR) of 0.088 dB to generate simulation data [29, 30], and the vertical and horizontal coordinates of the simulated data and the measured data indicate the number of pixels. The unit of unwrapping phase value is .
where is the effective value of the wrapped phase signal and is the power of the noise signal.
None of the algorithms in the PU experiment of the simulated data used the prefiltering algorithm. Except for the UKF and IUKF algorithms, other algorithms used the postfilter processing. The following are the unwrapped images and error distribution images of the simulated data of each algorithm. The error distribution images are obtained by subtracting the unwrapped phase from the reference phase.
Figures 6–12, (a) are the unwrapped phase images of each algorithm, and (b) are histogram of the between the unwrapped phase and the reference phase. Figures 6–12 are data images of PU using the BC, QG, ILS, MCF, UKF, IUKF, and RCBC methods. Comparing the error range distribution diagrams of each algorithm, it can be found that most of the errors of the RCBC method are concentrated near 0, and the absolute value of the error is less than 6 rad. It can be seen that compared with the four methods of BC, QG, ILS, and MCF, the RCBC method reduces the PU error and improves the PU accuracy. Compared with the two methods of UKF and IUKF, the RCBC method of accuracy has been slightly reduced. Because the UKF and IUKF methods merge the PU and filtering process, their accuracy of PU is high, but the methods are inefficient. Figure 12(c) is the branch-cut image connected by the RCBC method. Figure 12(d) is the wrapped phase image. From the residue density of Figure 12(c) and the wrapped phase image of Figure 12(d), it can be seen that the wrapped phase contains high noise.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
Table 1 shows the parameter comparison between BC method and the RCBC method in the simulated data experiment. Where “residue” is the number of residues, “length” is the total length of the branch cuts, and “ratio” is the percentage of the number of residues in the total number of the branch-cut pixels. It can be seen from Table 1 and Figure 12(c) that the branch cuts connected by the RCBC method are very short, and the proportion of nonresidues in the branch cuts is low.
Experiments have proved that the algorithm time-consuming changes with the SNR is small, and the general change range is less than 10% of the total time. Therefore, we only give the time of each algorithm in the two dimensions when the SNR is 0.088 dB. Table 2 shows the root-mean-square (RMS) error between the unwrapped phase and the reference phase of each algorithm under different SNR conditions, and the unit is rad. The time used by each algorithm is displayed in Table 3.
It can be seen from Tables 2 and 3 that the RCBC method takes longer than the BC, ILS, and MCF methods, but the accuracy is greatly improved. Compared with the QG method, it has advantages in time and accuracy. Compared with the UKF and IUKF methods, the time is greatly reduced when the accuracy is slightly reduced. Through the comparison between the above unwrapped images and the experimental data, it can be found that the RCBC method has high unwrapping accuracy in high-noise regions. As the noise increases, the RMS error of unwrapping slowly increases without sudden changes, so the method has strong robustness.
4.2. Measured Data
The measured data used in the paper is the interference image of Italian volcano Etna in Italy, provided by SIR-C/X SAR [19]. The following is the unwrapped phase images of each algorithm.
Figure 13(a) shows the wrapped phase image, and Figure 13(b) is the coherence map image. Figures 13(c)–13(i) are the unwrapped phase images of the BC, QG, ILS, MCF, UKF, IUKF, and RCBC method, respectively. The time used by each algorithm is presented in Table 4. Comparing Figures 13(c)–13(i), it can be seen that the two red wire-frames in the figures are the areas where the PU of each algorithm has a large error. From the unwrapped phase image in the red wire-frame in the lower left corner of Figure 13(i) is very smooth, and there is a shorter line with larger color difference in the red wire-frame on the right. Compared with other algorithms, the RCBC method has a smaller area of poor unwrapped effect. Observing Table 3, the RCBC method has a slight increase in time than the BC, MCF, and ILS method, but the accuracy is greatly improved. Compared with the QG method, time and accuracy have absolute advantages. Compared with the UKF and IUKF methods, the RCBC method has smaller unwrapped error area and higher efficiency. It can be seen from the above that the RCBC method also has high unwrapped accuracy in high-noise areas, requires less time for PU, has a certain practicability, and has strong robustness.
(a) Wrapped phase
(b) Coherence map
(c) BC method
(d) QG method
(e) ILS method
(f) MCF method
(g) UKF method
(h) IUKF method
(i) RCBC method
Sentinel-1 data is used for the second measurement data, and the SNR is 0.088 dB. The following Figure 14(e) is the unwrapped phase image of the RCBC method.
(a)
(b)
(c)
(d)
(e)
Where (a) is the wrapped phase image with noise, (b) is an enlarged view of the red wire-frames mark in (a), (c) is branch cuts, (d) is the true phase image, and (e) is the unwrapped phase image. This RMS error is 0.512 rad and takes 178 s.
5. Conclusion
The RCBC method classify the branch cuts and arranges various types of branch cuts through different methods, which greatly reduces the density of residues that need to be processed when connecting each type of branch cuts. In addition, the total length of the long branch cuts connected by the MCF method is also very short. Therefore, the RCBC method not only reduces the “island” but also shortens the length of the branch cut, and also reduces the error transmission, and enhances the antinoise performance of the method, so that the method maintains a high phase unwrapping accuracy in the high-noise regions. Then, the robustness of the method is enhanced by the postfiltering. In addition, a large number of residues are consumed when connecting the short branch cuts, which reduces the number of residues in the input of the MCF method, reduces the complexity of the method, and saves computing resources. The algorithm step before integration takes the time 2/3 of the MCF method. However, the method still has shortcomings, such as in the MCF method, regular block and fusion are used, and the repeated area is 15%, which reduces the accuracy of the method. Therefore, the next research direction is to find new block fusion methods and new integration methods, so as to reduce the time consumption of the method without loss of accuracy.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflict of interest.
Acknowledgments
This work was supported by the Dean Project of Guangxi Key Laboratory of Wireless Broadband Communication and Signal Processing and Guangxi Special Fund Project for Innovation-Driven Development (No. GuikeAA21077008).