A Refined Phase Unwrapping Method for High Noisy Dense Fringe Interferogram Based on Adaptive Cubature Kalman Filter

Cubature Kalman ﬁlter phase unwrapping (CKFPU) is an eﬀective algorithm in unwrapping the interferograms. The local phase slope estimation is a key factor that aﬀects the unwrapped accuracy. However, the estimation accuracy of local phase slop is relatively low in high noisy and dense stripes areas, which usually leads to the unsatisfactory unwrapped results. In order to eﬀectively solve this issue, the rewrapped map of the unwrapped phase (obtained by CKFPU algorithm), which is a ﬁltered interferogram with clearer fringes and more detailed information, is proposed in this paper to improve the phase slope estimation. In order to solve the problem of imprecise error variance for the new phase slope estimation, an adaptive factor is introduced into the CKFPU algorithm to increase the stability and reliability of the phase unwrapping algorithm. The proposed method is compared with the standard CKFPU algorithm using both simulated and real data. The experimental results validate the feasibility and superiority of the proposed method for processing those high noise dense fringe interferograms.


Introduction
Phase unwrapping is one of the most critical steps in InSAR data processing and its precision directly affects the precision of elevation or deformation measurements [1][2][3][4]. In recent years, numerous phase unwrapping methods have been proposed that can be generally divided into two types. e first type is to calculate the unwrapped phase based on path-following or the minimum norms' theory, such as branch-cut, region-growing, minimum discontinuity, minimum cost flow (MCF) networks, weighted leastsquares, and unweighted least-squares [5][6][7][8][9][10][11][12][13]. For these methods, noise filtering is necessary before unwrapping the noisy interferograms (usually called prefiltering). e second type is filter-based unwrapping methods, such as extended Kalman filter phase unwrapping (EKFPU), unscented Kalman filter phase unwrapping (UKFPU), and cubature Kalman filter phase unwrapping (CKFPU) [14][15][16][17][18][19][20][21][22][23][24][25][26][27]. For high noisy and dense stripes areas in interferograms, the prefiltering step can usually introduce additional problems. e information of pixels in high noisy and dense stripes areas is usually lost or distorted and cannot be recovered. e filter-based phase unwrapping algorithms that result in simultaneous noise filtering and phase unwrapping have been considered to be appropriate methods [16-18, 20, 27]. Since this type of algorithm is actually a state estimation method and the noise can be removed at the same time of phase unwrapping, it can avoid the phase loss and the distortion caused by prefiltering before phase unwrapping.
e Kalman filter phase unwrapping method was proposed by Krämer and Loffeld in 1996 [14]. Since then, several scholars have made numerous improvements, especially in the mathematical model for observation equation (sine and cosine function) and the stochastic model [16][17][18][19][20][21][22][23][24][25][26][27]. For the former, the algorithms have been improved in accuracy from the nonlinear function perspective, such as UKFPU, CKFPU, and relative algorithms [20,[22][23][24][25][26][27]. For the latter, the improvement is mainly reflected in the statistical model of noise, which is not assumed to be Gaussian distribution, such as particle filter phase unwrapping (PFPU) and unscented particle filter phase unwrapping (UPFPU) algorithms [19,21]. Additionally, filter-based phase unwrapping algorithms have also been applied in multibaseline InSAR and three-dimensional phase unwrapping (for DEM generation) [15,18,25,26]. In fact, the determination of the statistical model is a difficult problem due to the complexity of the measured noise data. So far, the theoretical and applied research of filter-based phase unwrapping algorithms has mainly focused on EKFPU, CKFPU, and UKFPU, with a Gaussian model hypothesis for noise, and has achieved good results [16,20,[22][23][24].
e CKFPU is considered to be an effective and simple algorithm for interferograms [23,27]. However, in high noisy and dense stripes areas, the accuracy of local phase slope estimation is relatively low, which usually leads to unsatisfactory unwrapped phase [23,27]. eoretically, the local phase slope estimation is a key factor that affects the accuracy of unwrapped phase. erefore, it is extremely important to improve its accuracy. It is generally known that filtering is a good way to improve the accuracy of phase slope estimation for the high noisy fringes. However, for dense fringe areas, filtering the interferograms before phase unwrapping usually leads to distorted fringe information. erefore, the main objective of this study is to propose a proper strategy that can retain as much detailed information as possible when the interferogram is filtered. e filterbased phase unwrapping algorithms can simultaneously perform noise filtering and phase unwrapping, thus usually retaining more details compared to the prefiltered phase unwrapping methods [27]. In this paper, the CKFPU filtering algorithm is employed first. en, the rewrapped map of the CKFPU result is used to estimate the local phase slope. However, the error variance σ 2 w of the new phase slope estimation method is usually not explicit (empirical model) [16,[19][20][21]. In order to avoid this problem, an adaptive factor is added to the CKFPU model to increase the stability and performance of the algorithm. e effectiveness of the proposed algorithm is validated by experimental results using both simulated and real data. is remainder of this paper is organized as follows. Section 2 introduces several primary preliminaries for phase unwrapping algorithm. Section 3 elaborates the improved phase slope estimation method and the adaptive CKFPU phase unwrapping algorithm. Simulated and real data tests are carried out in Section 4, while discussion and conclusion are presented in Sections 5 and 6, respectively.

Preliminary
2.1. Phase Unwrapping eory and Model. For phase unwrapping, the objective is to obtain the unwrapped phase φ(k) from the modulo 2π mapped phase φ(k).
If there are no phase errors, the phase difference between two adjacent pixels is is the true discrete phase slope and |2π is the modulo 2π operation. Under the assumption that |δ φ (k)| ≤ π, equation (1) is changed to Accordingly, the true unwrapped phase can be easily calculated by For most practical interferograms, especially obtained for high noise areas, phase errors always exist. e phase difference is always as follows: where ε φ (k) refer to the wrapped phase errors. erefore, the unwrapped phase calculated by (3) is not reliable because the phase slope is biased based on equation (2). Translating the phase unwrapping to state estimation is a suitable method [16,[19][20][21]. e phase unwrapping model is as follows [16,23]: where x(k) is the estimated unwrapped phase, y(k) is the measurement vector containing real and imaginary parts of the complex interferogram pixel k, δ → φ (k) and w(k) are the phase slope estimation and estimation error, respectively, and v 1 (k) and v 2 (k) are the measurement errors.

CKF Phase Unwrapping
Algorithm. For the unwrapping model shown in equation (5), the CKFPU is considered to be an effective phase unwrapping method [27]. In this algorithm, w(k) and v(k) are assumed to be Gaussian white noises.
e steps of one-dimensional CKFPU algorithm are as follows: (1) One-step prediction: the one-step prediction value x →− (k + 1) and its prediction error variance P − xx (k + 1) are 2 Mathematical Problems in Engineering where x → (k) is the estimated unwrapped phase of pixel k, P + xx (k) is its estimation error variance, and Q(k) is the error variance of the phase slope.
(2) Measurement update: firstly, the cubature points and the weights are calculated as Secondly, the error covariance is factorized by the Cholesky decomposition method as follows: irdly, the cubature points are obtained by Next, the propagated cubature points are calculated as en, the predicted measurement y →− (k + 1), error covariance matrix P − yy (k + 1), and cross-covariance matrix P − xy (k + 1) are as follows: where R(k) is the covariance matrix of the measurement noise. us, the state estimation x → (k + 1) and its covariance matrix P + xx (k + 1) are expressed as where K(k) is the cubature Kalman gain matrix. In this paper, the two-dimensional CKFPU is combined with a path-following strategy whose quality index is Fisher Distance and omnidirectional local phase slope estimation considering the quality of all the pixels in a 3 × 3 window. Details can be referred to [20,23,27].
From the whole process, it is clear that the local phase slope estimation is a key step that affects the accuracy of the unwrapped phase [16,[19][20][21].

Local Phase Slope Estimation.
Under the assumption that InSAR signals are stable in a small local area (B m × B n ), a complex sinusoidal signal can be expressed as [20,[28][29][30] where m � 1, 2, . . . , B m and n � 1, 2, . . . , B n ; m and n refer to the row and the column indices, respectively; B m and B n are the row and the column lengths of the local window, respectively; a(m, n) is the random amplitude of the complex signal; f x and f y refer to the true local frequency in the row and the column directions, respectively; w(m, n) is the corresponding noise. It is almost impossible to obtain the true phase slope directly from the noisy interferograms, and the higher the noisy, the lower the estimation accuracy. e local frequency is estimated by a commonly used maximum likelihood estimator in this paper. e phase slope from pixel (a, s) to pixel (m, n) in a small area is calculated by where f → x(a,s) and f → y(a,s) are the estimated frequency for pixel (a, s). Accordingly, the estimation error variance σ 2 w is Mathematical Problems in Engineering 3 where , r is the coherence coefficient. Noise is the main factor affecting the accuracy of the estimated frequency [27]. If there is little noise before performing maximum likelihood estimator, the estimated frequency is usually satisfying. However, it is an intractable problem that reduces the noisy before frequency estimation for those high noise dense fringe areas. at is because the original fringe information can usually be distorted.

Improved Phase Slope Estimation
Method. Equations (6)- (11) show that the phase slope estimation is one of the most important factor that affect the accuracy of unwrapped phase, while noise is the main factor [27]. For high noise dense fringe interferograms, if the original fringe information can be retained as much as possible after noise is filtered, the accuracy of fringe estimation can be improved. e filter-based phase unwrapping algorithms can simultaneously perform noise filtering and phase unwrapping, thus usually retaining more details compared to the prefiltered phase unwrapping methods [27]. erefore, the rewrapped maps of the filter-based unwrapped maps (here, the CKFPU unwrapped phase is selected) are the refined interferograms that can be used to improve the phase slope estimation. e rewrapped phase can be obtained as where φ(m, n) is the CKFPU unwrapped phase of pixel (m, n) and φ(m, n) is the corresponding rewrapped phase. en, the local phase slope estimation is implemented based on the rewrapped maps.
Figures 1(a) and 1(b) show the original noisy interferogram and the rewrapped map of the CKF unwrapped phase, respectively. It is obvious that the rewrapped map is much clearer, and the stripe shape is almost the same as the original even in the dense areas. us, the phase slope estimation δ → ϕ (k) can be more precise. Based on equations (13) and (14), the estimation error variance σ 2 w should be updated correspondingly because of the coherence coefficient r. For a single interferogram, pseudocoherence coefficient r is used and calculated as where φ i,j is the rewrapped phase of pixel (i, j) and k is the diameter centered on pixel (m, n). It should be pointed that the calculation of estimation error variance σ 2 w is approximate and there is no exact model yet. In order to address this problem, an adaptive CKF phase unwrapping algorithm is presented.

Adaptive CKF Phase Unwrapping Algorithm.
For the filter-based phase unwrapping method, the measurement model is clear (sine and cosine function) and the observation vector is assumed to be reliable. According to the literature, the precision of the unwrapped phase (state estimation) mainly depends on the phase slope estimation in the state model [27]. For the improved phase slope estimation method in Section 3.1, the error variance σ 2 w is not explicit. erefore, an adaptive factor is added to the CKFPU model to increase the stability and performance of the algorithm. According to the flow of CKF algorithm, the Kalman gain matrix K(k + 1) is the main influence on the precision of the estimation, and it depends on the covariance matrix P − yy (k + 1) and cross-covariance matrix P − xy (k + 1). us, an adaptive factor is applied to modify P − yy (k + 1) and P − xy (k + 1). e detailed steps are as follows [31].
Its theoretical covariance matrix is P − yy (k + 1) and the estimated covariance matrix is V T . e adaptive factor α k is expressed as where ΔV k is calculated as and c is an empirical constant with the optimal value of 1.0 [31]. It is obvious that α k ≤ 1. en, the adaptive process is as follows: 4 Mathematical Problems in Engineering x In case of the adaptive CKFPU process, if the phase slope estimation in the state model is not reliable and the noise statistics are not accurate, the adaptive factor can automatically adjust the weight of the prediction value x →− (k + 1). In general, the estimated covariance matrix T will be greater than the theoretical covariance matrix, and the adaptive factor α k is less than 1, which decreases the weight of the prediction value x →− (k + 1) and increases the weight of the measurement value y(k + 1).

Experiments
In this section, both simulated and real data are used to validate the effectiveness of the proposed methods. Compared with the prefiltered phase unwrapping methods, the CKFPU is an appropriate method for high noise dense fringe interferograms [23,27]. erefore, it is used in the comparative experiment in this paper. In order to ensure a fair comparison, Fisher Distance (FD) was selected to guide the unwrapping path and the maximum likelihood method was used to estimate the phase gradient [27,33]. e statistical characteristics of the process and the observation noises are given in Section 2. e proposed method is called improved adaptive cubature Kalman filtering phase unwrapping (IACKFPU) in this paper. e estimation results are evaluated with two indexes. One is the difference d between the unwrapped phase φ u and the priori phase φ p [32], which is defined as e other is the misfit and is defined as where N is the number of total samples and σ is the standard deviation, which is used to normalize the misfit [33]. Obviously, the smaller the value of the both indexes, the better the algorithm performance. For a simulated interferogram, the priori value is the true phase before wrapped. For the real data of topography interferogram, the simulated unwrapped phase based on a preexistent DEM is used as the priori value. Additionally, for the deformation interferogram, the collected and investigation materials and GPS measurement data are used for comparison.   Table 1. e scene shown in Figure 2(a) was generated by the peaks function in MATLAB and had large terrain fluctuation and high steepness. e true phase, the interferogram with noise, and the coherence map for the scene are shown in Figures 2(b)-2(d), respectively. e noise was also generated using the software, reflecting the degree of geometrical incoherence.
at is, the more serious the geometric incoherence (the denser the strips), the greater the noise. It is clear that the data are characterized by complex and dense stripes. e unwrapped phase, the rewrapped phase, and the error phase maps of the two methods (CKFPU and IACKFPU) are shown in Figure 3. e error maps are calculated by the first index d. For the IACKFPU method, the local phase slope estimation is based on Figure 3(c). It can be seen that the fringe details of Figures 3(c) and 3(d) are preserved well and a lot of the noise has been removed compared with the original interferograms (Figure 2(c)). It can be seen from Figures 3(a), 3(b), 3(e), and 3(f ) that, in areas where the stripes are sparse and the noise is low, the unwrapping results are continuous and the errors are small for both methods. In other words, the pixels with large errors are distributed in areas with dense fringe and high noise. Additionally, the distinctions of the unwrapped results are obvious for the two methods, mostly distributed in areas with dense fringe and high noise (the red ellipses in Figures 3(a) and 3(b)). Moreover, it can be seen from Figures 3(e) and 3(f ) that the errors for the IACKFPU method in most areas with dense fringe and high noise are decreased compared with the CKFPU method (marked by black ellipses). It should be pointed out that there are a few areas where the stripes are quite dense and even incoherent appearing opposite (marked by purple rectangles). It can be seen from the rewrapped map that the stripes are aliasing and their shapes change rapidly. us, the gradient estimation from the rewrapped map is not reliable. Actually, neither the CKFPU nor the IACKFPU method can obtain satisfying unwrapped phase in those incoherent areas and even the proposed method can obtain the worse result.
at is to say, the proposed method cannot bring advantage into play in those incoherent areas but can improve unwrapped precision in those coherent areas with dense fringe and high noise. Table 2 and Figure 4 show the quantitative results to further compare the unwrapping performance of the CKFPU and the proposed method. Table 2 shows the expectation and the standard deviation of the two indexes (d and χ 2 ). Mean(·) represents the expectation and Std(·) is the represents standard deviation. From equations (25) and (26), it is known that d is in radians and χ 2 is a nondimensional quantity. e error and the misfit statistical results in Table 2 clearly demonstrate that the IACKFPU method is superior to the CKFPU method. e results in Table 2 are calculated based on all pixels of the unwrapped and priori phase maps, which can reflect the performance of the algorithms in general.
Furthermore, misfit values along the unwrapping path can reflect where the errors occur and how the errors propagate. An optimal unwrapping algorithm should not create any large misfits in the beginning of the path to avoid or reduce early error propagation. at is to say, the latter the large misfit value appears, the better the algorithm is. Figure 4 shows the misfit curves along the path guided by FD for both CKFPU and IACKFPU methods. It can be seen from the figure that the misfit values of the two methods are extremely small and basically similar before about the 4 × 10 4 th pixel. After about the 4.525 × 10 4 th pixel, the misfit values of the IACKFPU method are much smaller than that of the CKFPU method. Moreover, the bigger misfit values of the IACKFPU method appear later than the CKFPU method (the red line is increasing later than the blue one marked in black rectangle in Figure 4). Figure 5 shows the phase unwrapping order in twodimensional map, which can more clearly reflect the unwrapping order of each pixel. e redder the pixel is, the later it is unwrapped. Figures 2(c) and 4 indicate that the pixels in areas with sparser stripes are almost unwrapped earlier and their misfit values are smaller. is is because the phase quality is determined by the degree of geometrical incoherence for the simulated data (the sparser the stripes are, the higher the phase quality is). e above results and analysis demonstrate that the areas with sparser stripes have better phase quality and they are unwrapped first. Both methods show about the same precision in these areas. is is because, for the stripes in high quality areas, the gradient estimation based on the original interferogram is as good as the result based on the rewrapped map. erefore, both methods show about the same precision in regions with good phase quality (areas with sparser stripes) and the error or misfit values are low. However, the IACKFPU method appears superior in areas with dense stripes. at is because the phase quality is extremely low, which can lead to unreliable gradient estimation based on original interferogram. In other words, the proposed method has obvious advantage in areas with high noisy and dense stripes.

Real Data.
In this section, two sets of experimental TerraSAR-X SAR data are used for testing the algorithms, which are covering Gujiao coal mining area, Taiyuan, acquired on January 4, 2013, and January 15, 2013. One is for topography and the other is for deformation. e parameters     Table 3. Most of the area is covered with extensive vegetation. So, in this paper, 2 × 2 multilooks of real data are performed. It has been pointed out that CKFPU method is superior to those prefiltered phase unwrapping methods for those high noisy and dense fringe interferograms [23]. en, the datasets are processed by the CKFPU and IACKFPU methods.

Dataset for Topography Interferogram.
e interferogram and its coherence map with 300 × 300 pixels processed by GAMMA software are shown in Figures 6(a) and 6(b), respectively. It can be seen that the interferogram is with high noise and the coherence is low in most areas. is is because the study area is covered with extensive vegetation and being featured by large mountains with steep slopes and ravines. An unwrapped phase map (Figure 6(i)) is simulated using a relief DEM with 4 m resolution, which is used to    e unwrapped phase maps, the error maps, and their histograms for CKFPU and IACKFPU methods are shown in Figures 6(c)-6(h), respectively. It can be seen that there are some differences in the unwrapped phase and the error maps, especially in the areas marked by the ellipses. e unwrapped phase in area circled by black ellipse based on the IACKFPU method is clearly more continuous than the result from the CKFPU method. In red ellipse marked area, it can be seen that there are fewer discontinuous small areas (small blue areas) for the IACKFPU method. In addition, the differences in the purple rectangles are also obvious. Combining the coherence maps, it can be seen that they located in the north-south low-quality band in the center of the maps. Because of the extremely low coherence, there is no important comparative significance. Moreover, from the histograms of the error maps (Figures 6(g) and 6(h)), the number of low error values for IACKFPU is larger than that for the CKFPU (circled by red rectangles). Furthermore, the misfit values along the FD unwrapping path for the two methods shown in Figure 6(j) and the statistical results for the evaluation indexes shown in Table 4 also verify the advantages and the reliability of the proposed IACKFPU method.

Dataset for Deformation Interferogram.
e data are covering the 18a203 working face which was active during the data acquisition period.
ere was about 20-30 cm (80-100 radians in the line of sight direction) of subsidence during this repeated cycle. From the interferogram and its coherence map (Figures 7(a) and 7(b)), it is clear that there is fast subsidence with severe decorrelation in the center of the working face and it has high noise. For DInSAR in this paper, a two-pass method was used and the relief DEM was selected as the external DEM. e preprocessing was implemented by GAMMA software to obtain the SAR interferogram and coherence map. e phase unwrapping step was completed by our own software. e unwrapped phase maps and their rewrapped maps for CKFPU and IACKFPU methods are shown in Figures 7(c)-7(f ), respectively. It is clear that the unwrapped phase for the IACKFPU method is superior to the result for the CKFPU method in terms of continuity and the unwrapping ability. e former obtains a more continuous unwrapped phase and a greater maximum subsidence (more than 75 radians) which is closer to the actual subsidence of working face (there is about 80-100 radians in the areas with dense fringes). In addition, the fringes are more continuous and clear in the edge of the subsidence area and more noise is removed for the IACKFPU method.

Discussion
Both the simulated and the real data experiments demonstrate that the proposed method is superior to the standard CKFPU method. Meanwhile, it can be seen that the superiority of the IAKFPU method for the topography data experiment is relatively less significant than the simulated data and the deformation data experiments.
is is because there is no significant change for the fringe density and the fringes are relatively sparse. From the coherence map ( Figure 6(b)), it can be seen that the coherence is not very low in most of the fringe area. us, this will not lead to extremely unreliable gradient estimation results based on the CKFPU method.
Specifically, both methods show about the same precision in regions with good phase quality. e reason for this is that the gradient estimation is almost the same based on the original interferogram and rewrapped map. Additionally, the IACKFPU method shows a little superiority in the area where the phase quality is slightly poor. is is because the relative reliable gradient estimation can be obtained based on original interferogram in the CKFPU method. When the phase quality is low, the performance superiority of the IACKFPU method becomes obvious due to the precise gradient estimation based on the noisy filtered rewrapped maps. In fact, if the interference stripes are with high noise but not density, other prefiltered methods can be used and it will not lead to distorted and aliased stripe information. However, for those high noisy and density stripe interferograms, the traditional prefiltered step usually causes distortion of some detailed information of the stripes, leading to worse gradient estimation results. While, the filterbased phase unwrapping method is a fine filtering method that can retain more details, and thus, the rewrapped maps can be used to improve the phase slope estimation. erefore, the proposed method is suitable for processing such high noise dense fringe interferograms and shows a significant superiority and reliability.

Conclusions
An improved adaptive cubature Kalman filtering phase unwrapping (IACKFPU) method for high noise and dense stripe interferograms is proposed. e rewrapped map of the filter-based unwrapped maps (the CKFPU unwrapped phase is selected in this paper) is used to improve the phase slope estimation, and an adaptive factor is added to the CKFPU model to increase the stability and the performance of the algorithm. e Fisher distance (FD) is used to calculate the quality maps to guide path tracking. e experimental results using both simulated and real datasets demonstrate that the proposed method has significant superiority and reliability and is suitable for processing high noise and dense fringe interferograms.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare no conflicts of interest.