Phase Cost Indicator-Guided Central Difference Information Filter 2-D Phase Unwrapping Approach

Phase unwrapping (PU) as an ill-posed inverse scientific problem has invariably been the significant difficulty in interferometric synthetic aperture radar technique for topographic mapping. Conventional PU approaches may produce a local or global unwrapping bias in severe phase changes caused by noise and large gradient terrain because the phase continuity assumption is not satisfied. To address this issue, we have created a novel phase quality-guided integrated filtering and unwrapping process for PU. In particular, the PU path is guided by an innovative phase cost indicator, which can automatically coordinate contributions from several phase quality estimation methods so that the unwrapping path can circumvent phase discontinuity regions and delay unwrapping these pixels. Since the central difference transform is another extra effective nonlinear transformation mode based on sigma points and only has one regulation parameter, we have employed the central difference information filter (CDIF) for the first time to realize the dynamic estimation of the unwrapped phase in the PU process. We tested the proposed approach using two real datasets with different mountainous areas. The results demonstrate that the PU accuracy of the proposed approach based on these two datasets is improved by 25.2% and 40.4% compared to the unscented information filtering (UIF) PU method, and the anti-noise capacity and unwrapping efficiency of the CDIF are both slightly superior to the UIF.

Phase Cost Indicator-Guided Central Difference Information Filter 2-D Phase Unwrapping Approach Yansuo Zhang , Yandong Gao , Shijin Li , and Shubi Zhang Abstract-Phase unwrapping (PU) as an ill-posed inverse scientific problem has invariably been the significant difficulty in interferometric synthetic aperture radar technique for topographic mapping. Conventional PU approaches may produce a local or global unwrapping bias in severe phase changes caused by noise and large gradient terrain because the phase continuity assumption is not satisfied. To address this issue, we have created a novel phase quality-guided integrated filtering and unwrapping process for PU. In particular, the PU path is guided by an innovative phase cost indicator, which can automatically coordinate contributions from several phase quality estimation methods so that the unwrapping path can circumvent phase discontinuity regions and delay unwrapping these pixels. Since the central difference transform is another extra effective nonlinear transformation mode based on sigma points and only has one regulation parameter, we have employed the central difference information filter (CDIF) for the first time to realize the dynamic estimation of the unwrapped phase in the PU process. We tested the proposed approach using two real datasets with different mountainous areas. The results demonstrate that the PU accuracy of the proposed approach based on these two datasets is improved by 25.2% and 40.4% compared to the unscented information filtering (UIF) PU method, and the anti-noise capacity and unwrapping efficiency of the CDIF are both slightly superior to the UIF.

I. INTRODUCTION
T WO-DIMENSIONAL phase unwrapping (PU) is a crucial procedure in the domains of interferometric synthetic aperture radar (InSAR) processors [1]. InSAR processor is commonly used to produce surface digital elevation model (DEM) [2], [3], [4], [5], which is typically used for three-dimensional modeling, contour extraction, terrain analysis, and other purposes. The interferometric phase measurement is mainly within the interval [−π, π), which is obtained by wrapping the actual phase in a nonlinear manner. Therefore, the measured interferometric phase is known as the wrapped phase, which differs from the true phase by 2kπ (k is an integer). The classical PU approaches often use local or global ways to get the absolute phase starting from the seed pixel according to the phase gradient. The premise of realizing these approaches is that the absolute value of the actual phase difference of any two adjacent pixels is less than π (i.e., phase continuity assumption) [1], [6]. However, large gradient terrain phase, noise and other factors will directly destroy the hypothesis, leading to a local or global unwrapping bias, which is the challenge that conventional PU approaches have always dealt with and is a major scientific drawback single-baseline PU approaches currently need to address. The single-baseline PU approach targets an interferogram, and its unwrapping efficiency is generally high, but it requires spatial continuity in the observation scene. However, the spatial variation pattern of the discovered object is sometimes uncontrollable, i.e., spatial heterogeneity. With the increase of interferometric data, the diversification of external information, and the promotion of artificial intelligence techniques, numerous PU approaches have flourished [7], [8], [9]. Commonly, the unwrapping accuracy of PU approaches in high noise and large gradient terrain regions have perpetually been the focus of PU researchers' efforts. For example, Yu et al. [10] proposed the two-stage programming approach (TSPA) multibaseline PU method. Under the framework of the TSPA, we can extend many representative single-baseline PU methods to the multibaseline domain, significantly increasing the unwrapping accuracy of traditional single-baseline PU approaches. Spoorthi et al. [11] proposed a novel framework for PU using the deep convolutional neural network (DCNN) known as PhaseNet. The proposed framework performs better than the quality-guided PU (QGPU) method, even under severe noise conditions with less computational time. Zhou et al. [12] also introduced deep learning into the PU field. They presented the PGNet methodology, which uses DCNN to predict phase gradient information, making it possible for the single-baseline PU method to directly get rid of the limitations of the phase continuity assumption. In the previous study [8], we proposed a PU method with external low-precision DEM-assisted phase slicing, which dramatically improved the quality of high-precision DEM products in large gradient regions produced by the InSAR processor. Wu et al. [13] proposed a Stereo-Radargrammetry assisted PU method to get the differential phase by subtracting the absolute phase converted from the pixel parallax information in Stereo-Radargrammetry from the InSAR-wrapped phase. Then they used the optimization-based global model to unwrap it, which improves the DEM reconstruction accuracy in steep regions. Li et al. [14] used gradient This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ information fusion to reduce the training and prediction costs of the deep learning unwrapping network while reducing the accumulation of gradient errors and improving the robustness of the PU method. The above studies provide examples of upgrading single-baseline PU methods using multibaseline, phase gradient estimation, external information assistance, and multitechnology fusion. These approaches do not necessarily require that the spatial phase of the observation scenario be continuous and can provide a reliable unwrapped phase. However, only a few PU methods can guarantee ideal unwrapping results when faced with complex interference fringes such as high uncorrelated noise or large gradient terrain.
Intriguingly, our recent study demonstrates that the QGPU method can circumvent instead of directly crossing the phase discontinuity region by establishing an appropriate unwrapping path [15]. It implies that the QGPU algorithm has the power to be not entirely forced by the phase continuity assumption by changing the integration order, which makes it possible for traditional single-baseline PU methods to avoid unwrapping error transmission in large gradient regions. A precise and dependable phase quality map is crucial to the QGPU algorithm. Currently, phase quality estimation primarily concentrates on residual, coherence, and phase derivative. The most frequently used phase quality estimation methods include the correlation coefficient (CC), the pseudo-correlation coefficient, the phase derivative variance (PDV), and the maximum phase gradient [6]. Osmanoglu et al. [16] examined the effects of unwrapping paths guided by five different phase quality functions on PU results and proposed a phase quality map based on Fisher information theory, which is more reliable than CC and second-order phase derivative. To more accurately reflect the texture features of interferograms, Liu et al. [17] proposed a phase quality map based on the gray level co-occurrence matrix, which has a subtle discriminatory ability with the CC map. Miguel et al. [18] proposed a phase quality estimation method based on residual theory and integrated it into the QGPU algorithm, primarily associated with the likelihood of 2 × 2 neighboring pixels becoming residuals once adding random tilt. It is not only computationally effective but also unaffected by the carrier signal of the initial phase. Sica et al. [19] introduced the interferogram and interferometric coherence as input features to the convolutional neural network models for PU. Using the coherence feature results in overall performance improvement on the unwrapping result. Similarly, Wang et al. [20] suggested using any of the most common quality maps as an input channel for GAUNet, incorporating the quality map into the formulation of the loss function to increase the precision of the phase ambiguity gradient estimation, and choosing the best PU path based on the residual quality map. The PU results demonstrate that the GAUNet-based QGPU algorithm makes full use of the advantages of the existing quality maps and performs effectively in monitoring deformations caused by earthquakes. Additionally, combining various current quality maps is a decent means for enhancing the accuracy of phase quality measurement [21], [22]. The research mentioned above has significantly improved the precision of phase quality estimation, but the majority hide the shortcomings of the typical quality indicators, which are still present and have not been addressed directly and effectively. Therefore, addressing the drawbacks of the frequently employed fast and accurate quality indicators will ensure that the unwrapping strategy can set a better unwrapping path based on a more reasonable phase quality value, which makes it possible to perform the QGPU algorithm by changing the order of integration without being entirely constrained by the phase continuity assumption.
To minimize the local transmission of unwrapping error in the QGPU algorithm, particularly in the large gradient region, we propose a phase cost indicator-guided central difference information filter (CDIF) PU approach, abbreviated as CDIFPU. The CDIFPU approach belongs to the integrated filtering and unwrapping algorithm and the QGPU algorithm, which involves two primary upgrade procedures. On the one hand, we change the unwrapping path to deal with the local transmission of unwrapping error in large gradient regions. On the other hand, we introduce CDIF to modify the dynamic phase gradient in real time on the unwrapping path. Specifically, optimizing the unwrapping path is achieved by improving the reliability of the phase quality estimation. We apply the Roberts gradient operator for the first time to compensate for the theoretical imprecision of the PDV indicator. Additionally, we directly perform variance data statistics on the wrapped phase, which can reflect the texture characteristics of the interference fringes and quickly identify the potential site of the phase jump. We refer to this novel phase quality estimation method as a phase cost indicator, which maximizes the PDV properties without relying on subjective parameterization. What is more, to the best of our knowledge, the central difference transform is combined with the information filter-based PU for the first time in this article, which significantly simplifies the setting of parameters of the unscented transform and makes the information filter unwrapping process more stable and efficient. We performed the experimental tests using TanDEM-X/TerraSAR-X long-baseline and large-area interferometric data. The results demonstrate that the CDIFPU approach is much more accurate than the representative PU methods and can significantly reduce local unwrapping bias in large gradient regions.
The rest of this article is organized as follows. Section II describes the fundamental theories of QGPU and its limitations. Section III provides a demonstration of the phase cost indicator-guided CDIFPU approach. Then, Section IV analyzes the performance of CDIFPU on the long-baseline and large-area interferometric datasets. Section V presents a discussion of the advantages and limitations of the proposed approach. Finally, Section VI concludes this article.

II. PHASE UNWRAPPING PROBLEM ANALYSIS
In this section, we will go through the specifics of how the QGPU approach works and analyze the reasons for the local transmission of unwrapping errors. Then the advantages and limitations of the unscented information filter PU (UIFPU) [23] are discussed. Finally, the shortcomings of the PDV indicator, which is most frequently used in phase quality estimation, are highlighted.

A. Quality-Guided Phase Unwrapping
The QGPU algorithm [24] is a representation of path following-based PU [25], [26] approaches. Compared with the optimization-based PU [27], [28], [29] methods, the path following-based PU approaches have two advantages. First, it can limit the unwrapping error to a local area rather than transmitting the PU result error globally. The second is to trust the measurement fully. In other words, the unwrapped phase is a one-step truth value produced by integrating the measured phase gradients without further correction. The core of the QGPU algorithm is that the unwrapping process diffuses from high-quality pixels to low-quality pixels with the help of a phase quality map. The specific process is as follows.
Step 1. Evaluate the phase quality of each interferogram pixel.
Step 2. Choose the highest quality pixel as the starting pixel for PU and put its neighboring pixels in an adjacency list.
Step 3. Arrange the adjacency list according to the quality of the pixels, then unwrap the pixel with the highest quality, and add its neighboring pixels to the adjacency list. Finally, remove the unwrapped pixel from the adjacency list.
Step 4. Repeat Step 3 until all the waiting pixels have been unwrapped. Establishing a reliable unwrapping path is extraordinarily necessary for the QGPU algorithms. Otherwise, it will likely cause local transmission of PU result errors. According to the analysis of the path following-based PU methods principle, if there is a phase integer leap pixel on the unwrapping path due to the large gradient terrain, the error of 2kπ will be introduced into the unwrapped phase of the pixels behind the path. For geographic mapping, the causes of phase integer leap include steep terrain and high noise. As shown in Fig. 1, the dark blue dotted line represents that the phase of the region is discontinuous and might jump periodically. If the pixels within the area are correctly identified as low-quality pixels, the unwrapping procedure will be performed as shown in PU path 1. If the pixels within the region are mistakenly thought to be high-quality, the unwrapping procedure can be carried out as shown in PU path 2. Compared to PU path 1, PU path 2 will introduce an error to the following unwrapped pixel from the misjudged pixel, causing the unwrapped phase of the local pixel to deviate from the actual phase. Therefore, the effective quality estimation in the phase discontinuous region is helpful for setting the path of the QGPU algorithm and subsequently enhancing PU accuracy.
The PU steps of the UIFPU are the same as those of the QGPU methods, except that the unscented transformation and information filter are introduced while unwrapping the highest quality pixel in the adjacency list. Since the unwrapping and filtering procedures are carried out simultaneously, the UIFPU can lessen the impact of noise on the unwrapping process and has higher PU accuracy than the standard QGPU methods. Information filter is an information form of Kalman filter, which uses an information matrix and vector to explain Gaussian distribution, and the computational complexity of the state updating procedure is relatively low. It is also the reason why the unwrapping efficiency of the UIFPU exceeds that of the unscented Kalman filter (UKF) PU (UKFPU). Since the Kalman filter PU is a nonlinear problem, it is simpler to estimate the probability distribution via unscented transformation than to calculate the nonlinear transform directly. The selection of a set of predetermined sigma points for nonlinear transformation is the basis of the proportional unscented transformation. The values of sigma points and their weights are determined by a group of parameters (α, β, κ), where α and κ control the sigma point scattering range, and β is used to introduce the high-order moment information of the probability distribution of random variables. The settings of those parameters are often found in paper [30]. However, while processing actual interferometric data, we could obtain more accurate PU results using nonrecommended parameter values, with the setting of α having the most significant bearing on the UKFPU. It indirectly indicates that the recommended parameter values will ensure the steadiness of the algorithm. However, it is not always the best choice for parameter values. Therefore, optimizing the sigma point sampling technique or simplifying the parameter setting is helpful to further perfect the PU theory of UKFPU or UIFPU approach and enhance the PU accuracy.

B. Phase Quality Estimation
The unwrapping path of the QGPU algorithm depends on a phase quality map. Currently, the PDV is regarded as the most trustworthy and widely used quality map. The PDV of the wrapped pixel is calculated by the variance data statistics of the phase gradient (i.e., partial derivative) within the central pixel and its neighboring pixels. Since the unit and magnitude of the variance are not easily explicable from the standpoint of phase quality, the standard deviation is frequently employed to describe the degree of statistical data in practical scenarios, but it is still represented by variance. Generally speaking, the PDV computation consists of two steps. The first step is to calculate the phase gradient of each interferometric pixel in azimuth and range as follows: where ϕ represents the wrapped phase, wrap(·) represents the wrapping operator, c represents the differential step of the wrapped phase, usually taking c = 1, ∇ x m,n and ∇ y m,n represent the azimuth phase gradient and range phase gradient of the pixel with m row and n column in the interferogram, respectively. The standard deviation of the azimuth and range phase gradients within the local window is determined in the second stage. Then, the PDV of the interferometric pixel is generated by adding the statistical standard deviation values. The mathematical expression is as follows: (2) where q PDV m,n represents the PDV of the pixel with coordinate (m, n) in the interferogram, and∇ x m,n and∇ y m,n represent the average phase gradient in the azimuth and range of the (2c + 1) × (2c + 1) local window.
The PDV represents the degree of deviation between the phase gradient and its mathematical expectation. Usually, a large PDV value shows that the phase quality is poor. The first-order phase derivative is calculated by wrapped phase difference, and a topographic phase change or high noise may cause its variations. In actual processing, the application of the PDV will always produce better PU results than conventional phase quality maps because the PDV is sensitive to phase gradient changes, especially in areas with large topographic changes.
Although the PDV is more reliable in PU than other standard phase quality maps, there may be better ones. We may intuitively see from (1) that the phase gradient is the rewrapped result of the interferometric phase difference value and does not perfectly replicate the initial wrapped phase information. It can be seen from (2) that the PDV procedure treats the azimuth and range phase gradient as two separate random variables for variance summation. There are two issues with this treatment. One is that the PDV quality map is directional since the range and azimuth phase gradient are directional. Second, the direct addition of the azimuth phase gradient and the range phase gradient is not the actual gradient. The issues mentioned above in PDV will have an impact on PU. If we address these issues pertinently, the unwrapping path of QGPU will be more reasonable, which will help resolve the local PU error transmission in large gradient regions.

III. PROPOSED APPROACH
In response to the abovementioned issues, a more effective phase quality estimation technique is proposed in this section. Based on the UIFPU approach, the concept of central difference transformation is introduced, and the phase cost indicatorguided CDIFPU approach is proposed.

A. Phase Cost Indicator
We prefer to improve the most widely used and trustworthy PDV indicator and propose a more effective phase quality estimation method. First, it is necessary to solve the problem that PDV cannot wholly reflect the information of the initial wrapped phase. Therefore, we tend to calculate the variance of the initial wrapped phase directly where q PV m,n represents the standard deviation of the wrapped phase of the pixel with coordinate (m, n) in the interferogram, andφ m,n represents the average of the interferometric phase within the (2c + 1) × (2c + 1) local window.
Second, we would like to solve the problem of the directionality of the PDV quality map. If the final phase quality map is non-directional, it is essential to ensure the phase gradient is isotropic when computing the phase gradient. In general, the gradient is directional, whereas the modulus of the gradient is isotropic. For any pixel in the interferogram, we can define the first-order differential Then its modulus is as follows: where grad(m, n) represents the modulus with an isotropic phase gradient. The core of PDV is to calculate the variance of gradient grad(m, n) = ∇ x m,n + ∇ y m,n . Obviously, the two are entirely different. To be more precise, we now obtain the phase gradient using the Roberts operator with rotation invariance rather than the gradient computation technique presented in (1). The Roberts operator, based on cross-difference, is widely employed for edge detection in digital images [22]. In contrast to the horizontal and vertical phase gradients given in (1), the Roberts operator calculates the phase difference in the ±45 • direction within the 2 × 2 neighboring pixels, then (4) can be further inferred as follows: It is worth noting that before acquiring the Roberts gradient of the edge pixel of the interferogram, it is necessary to extend the edge of the interferogram, and the position and mode of the extended pixels are indicated in Fig. 2.
Furthermore, the standard deviation of the Roberts gradient variance is expressed as follows: where q RGV m,n represents the Roberts gradient standard deviation of the pixel with coordinate (m, n) in the interferogram, and grad(m, n) represents the average of the Roberts gradient modulus within the (2c + 1) × (2c + 1) local window. Thanks to the fact that the Roberts gradient operator is more precise in calculating the phase gradient, the theory of Roberts gradient variance is more rigorous than PDV. To improve the effectiveness of the phase quality estimation, we integrate the PDV, phase variance, and Roberts gradient variance, which entails giving the various phase quality indicators the proper weights, as shown in (8). Notably, we should normalize the PDV, phase variance, and Roberts gradient variance to guarantee the consistency of the magnitudes before the fusion process where q m,n represents the fused phase quality estimation result, which can be directly used for unwrapping path settings;q PDV m,n , q PV m,n , and q RGV m,n represent the values of the normalized PDV, phase variance, and Roberts gradient variance, respectively; and ω 1 , ω 2 , and ω 3 represent the weights assigned to the normalized PDV, phase variance, and Roberts gradient variance, respectively. The phase assessment standard that incorporates the results of multiple phase quality estimation methods is defined as the phase cost indicator. It should be noted that the low phase cost of an interferometric pixel means that its actual phase quality is high. However, the high phase cost does not mean its real phase quality is always poor, possibly that its actual phase quality is high. The higher the phase cost of the interferometric pixel, the later it appears on the unwrapping path. In general, the assignment of weights is set by subjective experience, and the proportion of weights should be different for various interferograms. To prevent the phase quality estimation from being impacted by irrational human sets of the weight parameters, we adopt an allocation model to realize the automatic setting of the three weights. The allocation model depends on the root-mean-square error of the PDV quality map, the phase variance quality map and the Roberts gradient variance quality map. Since there is no exceptionally accurate phase quality estimation method and the PDV quality map has shown an excellent path guidance ability in the actual data PU processing, the PDV quality map is taken as the standard to calculate the root-mean-square errors of the three quality maps. The equation of the allocation model is as follows: where Re PDV , Re PV , and Re RGV represent the root-meansquare error of the PDV map, the phase variance map, and the Roberts gradient variance map, respectively. Taking the root-mean-square error of the phase variance map as an example, the equation is as follows: where M and N represent the total interferogram's number of rows and columns. It is obvious that due to Re PDV = 0, the weight ω 1 = 0.5 of PDV will always be the same, guaranteeing that the PDV is a significant component of the phase cost indicator. The abovementioned issues with PDV are exactly resolved by the phase cost indicator, which integrates the PDV, phase variance and Roberts gradient variance. As a result, the phase quality estimation is more reasonable and effective.

B. Central Difference Information Filter Phase Unwrapping
We create the CDIFPU algorithm based on the generic form of CDIF. Similar to UKF, the CDIF is also one of the Sigma-Point Kalman Filter algorithms initially proposed by Ito and Norgaard [31], [32]. The CDIF relies mainly on the Stirling center difference theory and uses polynomials to approximate nonlinear function derivatives [33]. It is an efficient nonlinear transformation technique widely used in applications such as navigation, attitude estimation, and battery state of charge estimation. The core of the central difference transformation is to use the first-and second-order difference operators to replace the first-and second-order components of the Taylor expansion of the nonlinear function. Since the CDIF has a slightly higher theoretical accuracy than the UIF and is regulated by only one parameter (i.e., the central difference step), we use the central difference transformation and the information filtering procedure to get the effective upgrading of the UIFPU approach. The specific process is as follows.
The CDIFPU approach continues to adopt the same status and measurement equations described in paper [30]. First, a set of sigma points is constructed from the prior unwrapped phase and its error variance where σ represents the set of constructed sigma points, ψ represents the prior absolute phase, P ψ represents the error variance of the ψ, L represents the dimension of the state variable, and S represents the central difference step and S ≥ 1. The parameter S determines the spread of the sigma points around the prior mean ψ. The optimal setting for S is equal to the kurtosis of the prior random variable. Because the status variable is unidimensional, it produces three sigma points in this case. If the status variables are consistent with a Gaussian distribution, the optimum value of S is √ 3. However, the distribution of the unwrapped phase is not purely Gaussian. After processing the actual data, we realized that the unwrapping effect is relatively ideal when adjusting S = 3/2. All sigma points are transferred by the linear state equation ς i = f (σ i ), i = 0, . . . , 2L. The one-step predicted unwrapped phaseψ, the prediction error variance Pψ, and the optimized prediction error variance Pψ of the wrapped pixel are calculated as follows: where u is the parameter to regulate the prediction error variance. Its recommended adjustment range is 0.4 ≤ u ≤ 1, and in this article, u = 0.6. Q is the process noise variance resulting from the phase gradient estimation error. w (m) and w (c) represent the mean and variance weights, respectively, which are defined as follows: where w represents the weight coefficients in calculating the covariance approximation's first-and second-order statistical properties. Based on the results of (12), the unwrapped phase prediction information matrix Ω and the amount of prediction information ξ can be calculated as follows: Next, the predicted measurement value ϕ and the covariance P ψ ϕ between the predicted state value and the predicted measurement value are computed by the nonlinear measurement After that, update the information matrix and the amount of information whereΩ andξ are the updated information matrix and the updated information vector, respectively, R represents the measurement noise variance, which can be computed from the CC, and υ represents the residual measurement (i.e., ϕ − ϕ). Finally, the unwrapped phase ψ of the waiting pixel and its estimation error variance P ψ are as follows: In summary, the specific implementation steps of the phase cost indicator-guided CDIFPU approach can be constructed as follows: Step 1. Use (8) to estimate the phase quality for each wrapped pixel.
Step 2. Estimate the horizontal and vertical phase gradients for each interferometric pixel in the frequency domain while calculating the Q and R for each pixel. The detailed computation process is mentioned in [30].
Step 3. Choose the seed pixel with the lowest phase cost indicator value. Determine its unwrapped phase ψ and the estimated error variance P ψ . Create an adjacency list and add the neighboring pixels above and below the seed pixel.
Step 4. Unwrap the pixel with the lowest value in the adjacency list based on (11) to (17). Then remove the pixel from the adjacency list, add its neighboring wrapped pixels to it, and update the list order. Repeat this step until all the interferometric pixels are unwrapped.

IV. PERFORMANCE ANALYSIS
The long-baseline and large-area interferometric datasets are used to test the performance of the phase cost indicatorguided CDIFPU approach proposed in this article. The PU results are compared and analyzed against other methods, including L 1 -norm [34], minimum-cost flow algorithm (MCF) [27], statistical-cost network-flow algorithm (SNAPHU) [35], UKFPU, UIFPU, PU max-flow/min-cut (PUMA) [36], PGNet, refined D-LinkNet UKFPU [15], and adaptive square-root UKFPU (ASRUKFPU), which was proposed in our previous paper [37]. The accuracy of the PU is mainly assessed through the root-mean-square error (RMSE) between the unwrapped phase and the reference unwrapped phase.

A. Long-Baseline Interferometric Dataset
First, we chose a region of Mount Hua in Weinan City, Shaanxi Province, China, as the area-of-interest, and its Google image is shown in Fig. 3(a). The landform of the experimental area is mainly mountainous, with apparent changes in terrain gradients and topographical discontinuities, making it challenging to test the performance of the proposed approach in this area. The radar interferometry data were derived from TanDEM-X satellite images. The long-baseline interferogram shown in Fig. 3(c) was obtained through SAR image registration and flat-earth phase removal. The baseline length is 370.46 m with a data size of 1000 × 1000 pixels. It is clear that the fringes of the interferogram are extremely dense. Consequently, this region has become a representative of the performance assessment of single-and multibaseline PU algorithms, as seen in papers [15] and [38]. Fig. 3(b) shows this region's shuttle radar topography mission (SRTM) DEM data. Fig. 3(d) shows the simulated unwrapped phase generated from Fig. 3(b) to evaluate the differences in unwrapping results across different PU methods. The long-baseline interferometric data is chosen because the longer the vertical baseline, the smaller the ground elevation difference corresponding to an interference fringe when the wavelength and side-looking angle are fixed. It indicates that the interferometric phase is sensitive to the terrain height, which is beneficial for producing high-precision DEM.
We have unwrapped the interferogram using the L 1 -norm, UKFPU, UIFPU, and CDIFPU. The PU results are provided in Fig. 4(a)-(d). Fig. 4(e)-(h) refer to the unwrapping error distributions. Notably, the UKFPU used the PDV quality map to guide the unwrapping path, whereas the CDIFPU used the phase cost indicator. The UIFPU performed two experiments with these two indicators, respectively. Since the MCF is equivalent to the L 1norm to some degree, the MCF unwrapping result is almost identical to Fig. 4(a). Therefore, we have not made statistics about the MCF in this experiment. From the unwrapping results, we can see that the three QGPU algorithms (i.e., UKFPU, UIFPU, and CDIFPU) can retrieve the absolute phase accurately, while the unwrapping results based on the optimization are pretty different from the actual phase, especially on the right side of the interferogram. This is primarily because the optimization-based method easily causes a global unwrapping error. Although the QGPU algorithm produces an unwrapping error transmission, the phase quality heap-sort strategy will guide the unwrapping path to bypass the phase discontinuity region in most cases, minimizing the unwrapping error transmission. It is important to emphasize that the PU results of the UKFPU and UIFPU still have two conspicuous unwrapping error regions, as shown in Fig. 4(f) and (g). To quantitatively describe the unwrapping accuracy of the different methods, we calculated the RMSE for the unwrapped phase, as shown in Table I. In addition, we tracked the filter unwrapping times for different PU approaches, as seen in Table I. The time we tracked did not include the phase gradient and quality estimation time. The statistical time is the unwrapping or iteration time pixel-by-pixel along the path. Note  Table I, it can be seen that the unwrapping error of the L 1 -norm is significantly larger than the other three methods, and the unwrapping accuracy of the UKFPU and UIFPU are comparable, but the unwrapping efficiency of the UIFPU is higher than that of the UKFPU. Thanks to the phase cost indicator proposed in this article, the CDIFPU approach almost avoids the transmission of unwrapping errors on a large scale. Its unwrapping accuracy is approximately 25.2% higher than the UIFPU and 81.6% higher than the L 1 -norm.
To see the effect of different improvements on unwrapping accuracy more clearly, we have replaced the PDV quality map of the traditional UIFPU method with the phase cost indicator. The PU result is presented in Fig. 5(a), and the distribution of the unwrapping error is shown in Fig. 5(b). The unwrapping accuracy of this procedure is 2.6356 rad. Comparing Figs. 4(g) and 5(b) combining the PU accuracy of the two methods, it can be observed that the phase cost indicator plays a decisive role in reducing the transmission of the unwrapping error and significantly improves the accuracy of the conventional information filter PU method. Essentially, the phase cost indicator is more accurate and reliable for labeling poor-quality pixels than the PDV, improving the unwrapping path to minimize the transmission of unwrapping error. Fig. 6(a) shows the PDV quality map based on Fig. 3(c) computed from (2). Fig. 6(b) is the initial wrapped phase variance quality map calculated from (3), where some red textures are visible, corresponding to the edges of the interference fringes. Compared to the PDV, the noticeable feature of the phase variance quality map is that it labels all the edges of the interference fringes as low-quality pixels. Fig. 6(c) shows the Roberts gradient variance quality map calculated from (7), which is more consistent with the PDV quality map, with the intuitive difference that it labels more regions as low-quality pixels. The quality map of the phase cost indicator integrating the three phase quality estimation methodologies is presented in Fig. 6(d). We have statistically plotted the histogram of its and PDV quality value distribution, as shown in Fig. 7. It is clear that in comparison with the PDV quality map, the phase cost indicator map labels more high-quality pixels (0-0.2) as low-quality pixels (0.2-0.5). The mislabeling of low-quality pixels as high-quality pixels often leads to an error transfer prematurely on the unwrapping path. Ultimately,  Fig. 3(d). (f) Difference between (b) and Fig. 3(d). (g) Difference between (c) and Fig. 3(d). (h) Difference between (d) and Fig. 3(d). the two primary defects of the PDV are directly addressed by using the Roberts gradient operator to calculate the first-order phase derivative and performing variance information statistics on the wrapped phase, resulting in the most reliable performance of the phase cost indicator. Comparing Figs. 4(h) and 5(b) with the PU accuracy of the two corresponding methods in Table I, we can see that the CDIF slightly improves the unwrapping accuracy compared to the UIF. Although the difference is not particularly noteworthy, the CDIF is restricted by only one parameter, which simplifies the operation and increases the approach's stability. Furthermore, the unwrapping efficiency of the CDIF is approximately 9.4% greater than that of the UIF. Additionally, we have especially tested the PU performance of the state-of-the-art PUMA based on the interferogram, as shown in Table II. It is worth mentioning that in our previous research [15], the PU accuracy of the PGNet method based on this data is 9.5272 rad, and the PU accuracy of the UKFPU based on D-LinkNet is 3.4857 rad, both of which are significantly lower than that of the CDIFPU approach. The unwrapped phase and PU error map of these representative approaches are shown in Fig. 8. It is apparent from the comparison that the CDIFPU  approach can obtain better unwrapped results than these methods in regions with large gradient terrain.
In summary, the CDIFPU approach can considerably increase the unwrapping accuracy than the UIFPU method and eliminate the unwrapping error transfer phenomena of the QGPU algorithm in areas with obvious terrain gradient changes. The phase cost indicator can more consistently guide the unwrapping path setting than the PDV, and the employed CDIF has a slight improvement in unwrapping accuracy and efficiency compared with the UIF, both of which are effective and meaningful improvements.

B. Large-Area Interferometric Dataset
To comprehensively demonstrate the performance of the phase cost indicator-guided CDIFPU approach proposed in this article, we have chosen the terrain in and around Jingle County, Xinzhou City, Shanxi Province, China, as another areaof-interest. The topography in the region is dominated by the loess plateau with highly fractured mounds. In contrast to the first set of test areas, we purposely chose a large area, covering approximately 739 km 2 . For a more accurate assessment of the performance of the different PU methods, we use the NASA DEM as a reference, which is reprocessed data from the SRTM with higher elevation precision. Fig. 9(a) shows the elevation of the NASA DEM in the test area after reverse geocoding. Fig. 9(b) presents a slope map for the test area divided into six  levels. We notice that the topography of the test area is mainly a slope of 15°-25°, and more steep slopes are larger than 45°on the left side. The large slope of the land and the significant slope change may reflect the complexity of the topographic phase, so the performance of the PU approach for unwrapping this region should be superior. We collected the experimental data from a dual-reception image captured in this region by the TanDEM-X/TerraSAR-X satellite. Fig. 10(a) is the intensity of the SAR image, where the data size is 2400 × 1600 pixels, with a range resolution of 10.91 m and a cross-range resolution of 17.65 m. Fig. 10(b) shows the wrapped phase after removing the flat earth phase, where it can be seen that the interference fringes are extraordinarily dense. The simulated unwrapped phase generated based on Fig. 9(a) is represented in Fig. 10(c). Fig. 10(d) is the quality map derived from the phase cost indicator. Comparing it with the slope map, we see that the phase cost of the large slope pixels is generally high and is mostly concentrated on the left side of the interferogram. In addition, we unwrap the interferogram using the MCF in GAMMA software, SNAPHU in Doris software, UIFPU and ASRUKFPU, respectively. The error range and  Table III. Fig. 10(e)-(i) show the unwrapping results of the different PU methods, and Fig. 10(j)-(n) indicate the differences between the unwrapping results and Fig. 10(c). It should be emphasized that we did not undertake any filtering operations before the PU. This can be explained by the fact that we have implemented a multilook operation and overfiltering will introduce a loss of phase information, making the unwrapping result less good than without filtering. Furthermore, when a random region of the interferogram is unwrapped separately, the different PU approaches' accuracy is significantly less than the overall unwrapping result. Therefore, choosing the large-area interferometric data is more compelling to test the methods outlined in this research.
On the whole, the five PU approaches successfully recover the absolute phase, but there is a significant unwrapping error at the edge of the left test area. This part of the unwrapping error is enlarged significantly, as shown in the red box in Fig. 10(j)-(n). The actual area corresponding to the box is approximately 92 km 2 . It is easy to observe that the error distribution of the unwrapping results is not similar within the box. Specifically, the unwrapped phase is primarily underestimated by the MCF and SNAPHU, and the MCF unwrapping error is smooth and continuous. Meanwhile, the UIFPU and ASRUKFPU underestimate and overestimate the unwrapped phase, but with a smaller error margin. Remarkably, the CDIFPU approach has virtually no unwrapping error in this region. In addition, Table III also shows that the phase cost indicator-guided CDIFPU approach has the smallest distribution of unwrapping errors, and the unwrapping accuracy is approximately 40.4% greater than that of the UIFPU.
Analyzing how the unwrapping bias changes along the PU path can give us some idea of where the significant error occurs. Following paper [16], we count the mean-square error (MSE) of all unwrapped pixels once every time a pixel is unwrapped so that the advantage of the unwrapping path of the CDIFPU approach is readily visible. When a large phase jump pixel suddenly appears on the unwrapping path, the MSE of all unwrapped pixels will abruptly increase and is then transmitted to the following wrapped pixels. An effective unwrapping path should delay or even not appear the abrupt MSE value to prevent the accumulation of unwrapping errors early. In this experiment, we statistically plotted the MSE variation of the unwrapped pixels along the direction of the unwrapping path for the PDV-guided UIFPU and the phase cost indicator-guided CDIFPU, respectively, as shown in Fig. 11. Since the error distributions of the ASRUKFPU and UKFPU unwrapped results are similar, and the unwrapping paths are guided by the same phase quality map, the MSE trend plots for both are similar, so it is not necessary to count ASRUKFPU. As can be seen on the overall trend chart, the large unwrapping error happens mainly after unwrapping 2.5 million pixels. In the partially enlarged view of Fig. 11, we have marked the four MSE mutation locations of the UIFPU. It is easy to see that the abrupt error at the first location occurs after a delay of about 100 000 pixels in the path of the CDIFPU approach, while the abrupt errors at the other three locations do not appear in the path of the CDIFPU. This indicates that the phase cost indicator can guide the unwrapping path to circumvent the pixels with phase discontinuities and delay unwrapping these pixels, thus avoiding the generation or premature accumulation To sum up, the accuracy of the phase cost indicator-guided CDIFPU approach proposed in this article is superior to the mature network flow PU algorithms and UIFPU in unwrapping the large-area interferometric data. The phase cost indicator proposed in this article can guide the unwrapping path to delay unwrapping phase discontinuity pixels and is a reliable method for phase quality estimation.

V. DISCUSSION
The experimental results above demonstrate that the phase cost indicator can better assist the PU path bypassing the phase discontinuity region than the PDV. The proposed CDIF unwrapping process has improved accuracy and efficiency compared to the UIF. The CDIFPU approach can address the unwrapping bias of conventional UKFPU and UIFPU in some local large gradient terrain regions, which is also a continuation of our previous study [15]. Based on the long-baseline interferometric dataset, we have demonstrated the effectiveness of the phase cost indicator and the CDIF unwrapping process. However, we did not clarify the enhancement effects of the Roberts gradient variance and the original wrapped phase variance in the phase cost indicator compared to the PDV. The phase gradient calculated by the Roberts gradient operator is undoubtedly precise. Nevertheless, this does not mean that the unwrapping accuracy guided by the Roberts gradient variance quality map is much better than that of the PDV. Essentially, both are gradient variance maps, which primarily reflect the interferogram phase gradient information and cannot reflect other additional information, such as the texture of the interference fringes. Therefore, although the Roberts gradient is much more accurate, its performance in guiding the PU path is very similar to that of the PDV. What is remarkably different is the original wrapped phase variance quality map that can reflect the texture features of the interference fringes and mark the edges as low-quality pixels. Consequently, the phase cost of the edge pixel of the interference fringe is lower than the PDV value, which is the core that the CDIFPU approach is challenging to occur unwrapping bias in large gradient terrain regions. Because the PDV has shown reliable performance in practical applications, we have retained most of its information and reinforced these values with Roberts gradient variance. It should be noted that this does not lead to duplication of data on gradients but rather highlights the importance of the phase gradient as an indicator of phase quality assessment. This is also why we chose the PDV quality map as the standard to calculate the RMSE of the three quality maps.
Recent deep learning approaches have also proven valid for PU, which also considers improving path integration through the phase quality parameter, especially coherence. Since we cannot reproduce all these methods, we chose the deep learning PU methods, such as PGNet and refined D-LinkNet, that we can implement for comparison. The two advanced learning-based PU approaches formulate the phase gradient estimation as a multiclasses segmentation issue, which still requires traditional algorithmic model-based PU methods for the iterative processing of phase gradients. Nowadays, many learning-based PU approaches with excellent performance have emerged, some of which are capable of predicting branch-cuts from the residual map and then performing post-processing, such as the deep learning-based branch-cut method (i.e., BCNet) [39]. Some can even perform a regression directly from the wrapped phase to the absolute phase in one step, such as the PU method based on the conditional generative adversarial network (i.e., PU-GAN) [40]. Theoretically, the phase cost indicator map presented in this paper can be introduced as a conditional variable in the PU-GAN generator and discriminator. In other words, we can use external phase quality information to influence the unwrapped phase generation process to enhance further the model generalization performance. Our future research work will continue exploring such a program's feasibility.
For topographic mapping, a better way to evaluate the PU accuracy is to convert the unwrapped phase to an encoded DEM and verify it using the Ice, Cloud, and Land Elevation Satellite data or ground control points. However, due to the limited control data available to us, this validation work has not been performed in the article. Compared to the simulated unwrapped phase generated by the existing external DEM, the performance of the PU method can also be meaningfully explained, which is the frequently used verification method in most current studies. Nevertheless, it has to be recognized that the CDIFPU approach proposed in this article is practical, and its presence makes the results of the single-baseline PU methods relatively good, especially in large gradient topography regions. In the meantime, we hope that the method proposed in this article will provide reference significance for the data processing of new InSAR satellites, e.g., L-band spaceborne bistatic system LuTan-1.

VI. CONCLUSION
PU is an essential procedure in the InSAR processor. As a long-standing research topic, how to minimize the PU error in large gradient terrain areas has always been the focus of InSAR topographic mapping. To this end, we propose a phase cost indicator-guided CDIFPU approach from the perspective of path optimization. It is worth emphasizing that the phase cost indicator may attribute a more reasonable phase quality value to each pixel, thus guiding the unwrapping path to circumvent the phase discontinuity region. Furthermore, we also would like to highlight that the CDIF is regulated by only a single parameter, and its accuracy and efficiency are slightly better than the UIF, which is more appropriate for PU. According to the PU results based on Mount Hua data, the CDIFPU approach improves the accuracy by 81.6% and 25.2% compared to the L 1 -norm and UIFPU, respectively, and by 72.5% and 24.9% compared to learning-based PU approaches like PGNet and refined D-LinkNet, respectively. The CDIF also preserves the accuracy of the UIF while enhancing PU efficiency by 9.4%. The unwrapping results based on data from Loess Plateau show that the CDIFPU approach improves the accuracy by 53.5% and 40.4% compared with MCF and UIFPU, respectively, and is capable of preventing or delaying the occurrence of PU errors. The state-of-the-art PU approaches are challenging to replicate due to the lack of open-source code and parameter-setting instructions. Nevertheless, in future studies, it will be necessary to compare the performance of the CDIFPU with other advanced PU methods, especially learning-based PU approaches. In addition, we will also explore the possibility of combining the phase cost indicator map with PU-GAN to enhance the generalizability of the one-step PU models.