BM3D Denoising for a Cluster-Analysis-Based Multibaseline InSAR Phase-Unwrapping Method

: Multibaseline (MB) phase unwrapping (PU) is a key processing technique in MB interferometric synthetic aperture radar (InSAR). As one of the most popular methods, the cluster analysis (CA)-based MBPU method often suffers from the problem of low noise robustness. Therefore, the block-matching and 3D ﬁltering (BM3D) algorithm, one of the most effective ﬁltering methods for image denoising, is applied to improve the performance of the method. Five different ﬁltering strategies for applying BM3D are proposed in the paper: interferogram ﬁltering (IFF), intercept ﬁltering (ICF), cluster number ﬁltering (CNF), unwrapped phase ﬁltering (UPF), and simultaneous ﬁltering (STF). In particular, while keeping the general structure of BM3D, four different similarity measures are deﬁned for interferograms, intercepts, clusters, and unwrapped phases to accommodate the special characteristics of different ﬁltering objects. Experiments on synthesized and real InSAR datasets prove their feasibility and effectiveness, and the experiment results show that (1) the PU accuracy and robustness of the CA-based MBPU method can be greatly improved by adding BM3D denoising; (2) simultaneous ﬁltering of interferograms, intercepts, cluster numbers, and unwrapped phases works best, but with the worst time complexity; (3) when ﬁltering is performed for only one object of the CA-based MBPU method, the ﬁltering effect of CNF and UPF is better than that of IFF and ICF; and (4), considering the three indicators of PUSR, NRSE, and time consumption, CNF and UPF should be the best choices.


Introduction
Interferometric synthetic aperture radar (InSAR) can acquire multiple synthetic radar images of an observed region and measure their interferometric phases to obtain the corresponding terrain height information (such as a digital elevation model, DEM) or deformation information [1][2][3][4][5]. However, the measured interferometric phase is wrapped into the principal value interval, which is modulo 2π. In order to obtain the correct absolute interferometric phase, phase unwrapping (PU) must be performed [6]. The accuracy of the target elevation is directly related to the PU accuracy, so PU has always been a key processing step in InSAR applications.
Researchers first proposed single-baseline phase-unwrapping (SBPU) methods to reconstruct the unwrapped phase, and several representative works are introduced in [7][8][9]. However, the traditional SBPU method must assume that the actual phase jump between adjacent pixels is less than π to ensure the uniqueness of the solution, which is called the phase continuity assumption. When the phase jump between adjacent pixels of the observed terrain is greater than π, the accurate phase-unwrapping results cannot be obtained using the SBPU technique.
Therefore, multibaseline phase-unwrapping (MBPU) algorithms are further proposed to overcome the limitation of SBPU. According to [6], there are two main categories of MBPU methods: parametric-based methods and nonparametric-based methods. The parametricbased methods construct a statistical framework according to the probability density function of the interferometric phase and regard the elevation or elevation difference as the parameters to be estimated. In [10], a PU method for reconstructing highly discontinuous ground elevation profiles using ML technology was proposed. In [11], an improved technique based on ML estimation was proposed to reconstruct the elevation information. In this method, a multifrequency independent phase dataset obtained by nonoverlapping band-pass filtering the interferometric SAR raw data pair was used, and the unknown surface was approximated by the local plane assumption. In [12], a maximum a posteriori (MAP) estimation method based on multifrequency interferometry and a Markov random fields (MRFs) model was proposed. In [13], the MAP estimation method was improved based on the total variation model and graph-cut-based optimization algorithms. In [14], an enhanced MBPU method by combining an unscented Kalman filter with an enhanced joint phase gradient estimator was proposed. In [15], an extended Kalman filter (EKF)based MBPU approach, which has the ability of handling phase discontinuities, was proposed. In [16], a nonlocal denoising method based on patch similarities and totalvariation regularization was proposed to realize resolution-preserving denoising and robust MBPU for DEM reconstruction.
The nonparametric-based methods do not need the probability density function of the interferometric phase and directly use multiple interferograms with different baseline lengths to reconstruct the absolute phase. In [17], H. Yu et al. proposed a two-stage programming approach (TSPA), which transplants the SBPU framework to MBPU. In [18,19], two refined TSPA methods were presented. The first method is abbreviated as LPM-TSPA [18], and it improves the performance of stage one of the TSPA by approximating the terrain surface to a local plane, which allows more interferometric phase information to be used to estimate the ambiguity number gradient. The second refined TSPA-based MBPU method proposed in [19] uses an unscented Kalman filter (UKF) to improve the performance of stage two of the TSPA. As one of the most popular methods, the cluster analysis (CA)-based MBPU method has also been widely studied. In [20], the ambiguity vectors corresponding to each pixel are solved according to the closed-form robust Chinese remainder theorem (CRT), and then clustering is performed according to the ambiguity vectors. In order to improve the noise robustness of the conventional CRT-based MBPU method, a fast CA algorithm was proposed in [21], which clusters pixels with the same ambiguity vector into a group according to the recognizable mathematical pattern, i.e., intercept information, and the ambiguity vectors of all pixels in the same cluster can be estimated through the group-by-group clustering center information. In [22], a refined CA-based method was proposed by linearly combining multiple InSAR interferograms with different baseline lengths, which makes the new combined interferogram have a larger ambiguity interval. In [23], the recognizable mathematical patterns were expanded into row, line, and intercept, which further improves the noise robustness of the CA-based method. In [24], a new closed-form robust CA-based MBPU and filtering algorithm (MBPUF) was proposed to improve the efficiency, height reconstruction accuracy, and noise robustness of the CA-based MBPU method.
In addition, the application of deep learning in PU is becoming more and more extensive. In [25], deep learning was applied to the CA algorithm, and an unsupervised deep convolutional neural network named CANet was proposed. In [26], a novel deep convolutional neural network (DCNN), abbreviated as PGNet, was proposed for estimating the phase gradient information instead of phase continuity assumptions. In [27,28], by reformulating the definition of the problem of directly obtaining the continuous original phase as the wrap count obtained on each pixel through semantic segmentation, a novel framework called PhaseNet that uses a deep, fully convolutional neural network to unwrap the phase network was proposed. In [29], a comprehensive overview of InSAR phase unwrapping based on artificial intelligence was provided, which reviews the SBPU and MBPU methods based on artificial intelligence. In addition to PU, deep learning also has many applications in SAR and InSAR, and readers can refer to [30][31][32][33][34].
As a representative MBPU method, the CA-based MBPU method has the following problems: (1) it has poor noise robustness and is sensitive to phase noise, i.e., small phase noises may cause large unwrapping errors [23]; (2) when the intercepts of the cluster with fewer pixels and the cluster with more pixels are very close, the two clusters will be merged into one cluster due to the presence of phase noise, resulting in wrong clustering results [23]; and (3) it does not have phase-filtering capability [24]. In order to solve the above problems, this paper introduces filtering techniques into the CA-based MBPU method. In SAR and InSAR applications, traditional methods mainly use filtering techniques for interferograms [35][36][37][38] or SAR images [39,40].
However, in this paper, we consider the application of the current state-of-the-art BM3D technology to the four objects involved in the CA-based MBPU method: interferogram, intercept, cluster number, and unwrapped phase. The main advantages of the BM3D technology are as follows [41]: (1) it has achieved the most advanced denoising performance in terms of peak signal-to-noise ratio and subjective visual quality; (2) it can be adapted to various noise models by modifying the coefficient variance calculation in the basic and Wiener parts of the algorithm; and (3) it preserves relatively complete information about the structure of the original image. Therefore, five different filtering strategies for applying BM3D are proposed: interferogram filtering (IFF), intercept filtering (ICF), cluster number filtering (CNF), unwrapped phase filtering (UPF), and simultaneous filtering (STF). In particular, while keeping the general structure of BM3D, four different similarity measures are defined for interferograms, intercepts, clusters, and unwrapped phases to accommodate the special characteristics of different filtering objects.
The rest of this paper is organized as follows. First, Section 2 reviews the principle of the CA-based MBPU method and points out the shortcomings of this method. Section 3 briefly reviews the structure of the BM3D algorithm. Section 4 introduces the application of the BM3D algorithm to the CA-based MBPU method. In Section 5, experiments on simulated and real data are introduced to prove the effectiveness and feasibility of the proposed method, and the results are discussed. Finally, the conclusions are given in Section 6.

Basic Principle of the CA-Based MBPU Method
For an InSAR system, the terrain height can be obtained from the unwrapped phase and other system parameters [18]: where h(p) is the terrain height of the p-th pixel, λ is the working wavelength of the InSAR system, r(p) represents the slant range of the p-th pixel for the main channel, θ o is the look angle, B i is the length of the orthogonal baseline for the i-th interferogram, and ψ i (p) is the flattened absolute interferometric phase of the p-th pixel for the i-th interferogram. The measured phase of the target obtained by the system, namely the wrapped phase, can be expressed by [7]: where ϕ i (p) is the wrapped phase of the p-th pixel of the i-th interferogram (from 0 to 2π in this paper) and k i (p) is the unknown integer ambiguity of the p-th pixel of the i-th interferogram. It can be seen from (2) that the purpose of PU is to find the correct value of k i (p) to obtain the correct unwrapped phase ψ i (p). For simplicity, only the case of two different baselines is considered for now. Then, the following equation can be obtained according to (1) and (2): i.e., It can be seen from (4) that the information of each pixel is reflected by a linear equation. Because B 1 and B 2 are constants, the slope of the straight line is also a constant, while the intercept varies with ϕ 1 (p) and ϕ 2 (p). According to [21], the term ambiguity vector [k 1 (p), k 2 (p)] is defined as the ambiguity numbers corresponding to a pixel in two different interferograms. Additionally, a cluster is defined as a set of pixels with the same ambiguity vector. Obviously, the intercepts of the pixels in a cluster are the same since their corresponding straight lines pass through the same integer point [k 1 , k 2 ] and with the same slope: B 2 /B 1 . Therefore, the cluster to which a pixel belongs can be determined by its corresponding intercept information.
Through the above analysis, the main steps of the CA-based MBPU method can be summarized as follows: (1) the intercepts of all pixels are calculated according to the input interferograms; (2) the cluster to which each pixel belongs is obtained by using the CA algorithm; (3) the cluster ambiguity vectors of each cluster are calculated according to CRT and the intercept of each cluster's centerline; and (4) the final unwrapped phases of each pixel are acquired by the cluster ambiguity vectors and wrapped phases. Although the CA-based MBPU methods can meet the requirements of PU to a certain extent, the poor noise robustness makes the algorithm prone to produce erroneous clustering results. It should be noted that when the clustering is finished, a natural number is assigned to each cluster, which is called the cluster number. A cluster number is a natural number we artificially set for each cluster, and the purpose is to number the clustering results. Different cluster numbers represent different clusters. Therefore, we can use Figure 1 to show the influence of noise on the clustering results.   Figure 1b is the noisy clustering result generated by the CA algorithm, and there are some wrong clustering results due to the existence of noise. In addition, although the CA algorithm improves the noise robustness and reduces the time complexity to a certain extent when the cluster with fewer pixels is very close to the cluster with more pixels, these two clusters will be merged into one cluster due to the presence of phase noise [23], resulting in erroneous clustering results. Therefore, the BM3D algo-  and there are some wrong clustering results due to the existence of noise. In addition, although the CA algorithm improves the noise robustness and reduces the time complexity to a certain extent when the cluster with fewer pixels is very close to the cluster with more pixels, these two clusters will be merged into one cluster due to the presence of phase noise [23], resulting in erroneous clustering results. Therefore, the BM3D algorithm can be applied to the CA-based MBPU method to improve PU accuracy. According to the main steps of the CA-based MBPU method, BM3D can be applied to the following four objects: (1) the interferograms; (2) the intercepts; (3) the cluster numbers; and (4) the unwrapped phases.

BM3D Algorithm and Similarity Measures Improvement
In this section, we first review the basic processing procedure and then define four similarity measures for the four different objects to be filtered in this paper: interferograms, intercepts, clusters, and unwrapped phases.

BM3D Algorithm
The basic flow chart of the BM3D algorithm is shown in Figure 2. According to [30], the BM3D algorithm was originally proposed for Gaussian white noise, which combines nonlocal methods with wavelet domain shrinkage and Wiener filtering. Then, the main steps of the BM3D algorithm can be described as follows [30].  (1) Block-matching grouping: Find blocks that are similar to the reference block and named as similar blocks, and the similarity is determined by calculating the predefined distance between the reference block and the block to be matched. Then, the (1) Block-matching grouping: Find blocks that are similar to the reference block and named as similar blocks, and the similarity is determined by calculating the predefined distance between the reference block and the block to be matched. Then, the similar blocks are stacked into a 3D array according to the similarity, which is called a 3D group. (2) Collaborative hard thresholding: Perform 3D transformation processing on the 3D groups, attenuate noise through hard thresholding, and then use 3D inverse transformation to obtain the estimated values of the 2D similar blocks. (3) Aggregation: Since there are often overlapping parts between the blocks, the result is the same pixel often being contained in several different blocks. Therefore, the pixel value of each pixel is repeatedly estimated, and a weighted average of these multiple estimates is needed to obtain the basic estimate for each block. (1) Block-matching grouping: Find the position of similar blocks by means of blockmatching in basic estimation. Block-matching grouping consists of two parts, which can obtain two 3D groups, one from the noise image and the other from the basic estimated image. (2) Collaborative Wiener filtering: Apply 3D transformation to the above two 3D groups, use the 3D group in the basic estimation as the energy spectrum of the real signal, use the energy spectrum to perform collaborative Wiener filtering on the noise image, and finally, the processed data are inversely transformed and returned to the original position of the similar block to obtain the final estimated value. (3) Aggregation: Because the blocks obtained after grouping and filtering may overlap each other, weighted average processing on pixels with multiple estimated values is performed to obtain the final estimation.

Similarity Measures Improvement
In the block-matching grouping of the original BM3D algorithm, the similarity between two blocks is usually expressed by the inverse of some distance measured between them [30]. The smaller the distance is, the more similar the two blocks are. Since there are four different processing objects involved in the CA-based MBPU method, interferograms, intercepts, cluster numbers, and unwrapping phases, four different similarity measures are defined to accommodate the special characteristics of different filtering objects.
For the interferogram, the similarity of two pixels can be measured by the cosine value of their wrapped phase difference: where ϕ 1 and ϕ 2 represent the wrapped phases corresponding to two different pixels, respectively. For the intercept, the similarity of two pixels is measured by the following formula: where b 1 and b 2 represent the intercepts of two different pixels, max 1<p<n (b p ) represents the largest intercept value among all the pixels, and min 1<p<n (b p ) represents the smallest intercept value among all the pixels. For the cluster, the greater the difference between two cluster numbers, the greater the difference between the two pixels, so the similarity of the two pixels can be measured by: where m 1 and m 2 represent the cluster numbers of two different pixels, max 1<p<n (m p ) represents the largest cluster number among all the pixels, and min 1<p<n (m p ) represents the smallest cluster number among all the pixels. For the unwrapped phase, the similarity can be measured by: where ψ 1 and ψ 2 represent the unwrapped phases of two different pixels, max 1<p<n (ψ p ) represents the largest unwrapped phase among all the pixels, and min 1<p<n (ψ p ) represents the smallest unwrapped phase among all the pixels. According to the definitions of Equations (5)-(8), the closer the value of S is to one, the stronger the similarity of the two pixels; on the contrary, the closer the value of S is to zero, the weaker the similarity of the two pixels.

BM3D Denoising for the CA-Based MBPU Method
Considering that it can perform denoising well and retain some original details of an image, the BM3D algorithm is applied to the CA-based MBPU method.
According to the main steps of the CA-based MBPU method described in Section 2, the BM3D algorithm is applied to the following four objects separately or simultaneously: (1) the wrapped phases, i.e., interferogram; (2) the intercepts; (3) the clustering results, i.e., cluster numbers; and (4) the unwrapped phases. There can be five different filtering strategies in total, and the corresponding basic flow charts are shown in Figure 3. Among them, the first four filtering strategies filter only one object, while the fifth filtering strategy filters all four objects at the same time. Moreover, all the five filtering strategies contain four main steps: data preprocessing, BM3D denoising, denormalization, and post-processing. In the following, we will describe these four main steps for applying the BM3D algorithm to the CA-based MBPU method in detail.

Data Preprocessing
The interferogram, intercept, cluster number, and unwrapped phase datasets to be filtered are first obtained according to the CA-based MBPU methods.
In order to uniformly use BM3D in the four objects of the CA-based MBPU methods, they need to be normalized, and the value range of these inputs should be changed to zero-one. Therefore, the normalization of these four objects can be implemented using the following four equations: where ϕ norm (p) represents the normalized wrapped phase of pixel p, b norm (p) represents the normalized intercept of pixel p, max 1<p<n (b(p)) represents the maximum intercept, min 1<p<n (b(p)) represents the minimum intercept, m norm (p) represents the normalized cluster number of pixel p, M represents the total number of clusters, ψ norm (p) represents the normalized unwrapped phase of pixel p, max 1<p<n (ψ(p)) represents the maximum unwrapped phase, and min 1<p<n (ψ(p)) represents the minimum unwrapped phase.

BM3D Denoising
The BM3D algorithm absorbs the processing idea of nonlocal filtering in the NL-means algorithm. In the filtering process, it further uses the information of block groups in the whole image. While keeping the general structure of BM3D, four different similarity measures are defined in Equations (5)-(8) for the four different processing objects involved in the CA-based MBPU method to accommodate the special characteristics of different filtering objects: interferograms, intercepts, cluster numbers and unwrapping phases. Moreover, five different filtering strategies for applying BM3D are proposed: interferogram filtering (IFF), intercept filtering (ICF), cluster number filtering (CNF), unwrapped phase filtering (UPF), and simultaneous filtering (STF). According to the definitions of Equations (5)- (8), the closer the value of S is to one, the stronger the similarity of the two pixels; on the contrary, the closer the value of S is to zero, the weaker the similarity of the two pixels.

BM3D Denoising for the CA-Based MBPU Method
Considering that it can perform denoising well and retain some original details of an image, the BM3D algorithm is applied to the CA-based MBPU method.
According to the main steps of the CA-based MBPU method described in Section 2, the BM3D algorithm is applied to the following four objects separately or simultaneously: (1) the wrapped phases, i.e., interferogram; (2) the intercepts; (3) the clustering results, i.e., cluster numbers; and (4) the unwrapped phases. There can be five different filtering strategies in total, and the corresponding basic flow charts are shown in Figure 3. Among them, the first four filtering strategies filter only one object, while the fifth filtering strategy filters all four objects at the same time. Moreover, all the five filtering strategies contain four main steps: data preprocessing, BM3D denoising, denormalization, and post-processing. In the following, we will describe these four main steps for applying the BM3D algorithm to the CA-based MBPU method in detail.

Interferogram Filtering (IFF)
Please refer to Figure 3a. In this method, the noisy interferogram is filtered directly by BM3D and the other steps are essentially unchanged. The wrapped phases are damaged due to the presence of noise, so filtering is performed on them. The noise interferogram is used as the input noise image of the BM3D algorithm to obtain the filtered wrapped phase ϕ f iltered (p), and then the filtered interferogram is used for the following other steps. This method is abbreviated as the IFF method.

Intercept Filtering (ICF)
Please refer to Figure 3b. Since the CA-based MBPU methods use the intercept as a clustering index, when they are contaminated by noise, the clustering results will be affected. Therefore, the intercepts contaminated by noise can be filtered using BM3D to obtain the filtered intercept b f iltered (p), and then the filtered intercept is used for subsequent processing. This method is referred to as the ICF method in this paper.

Cluster Number Filtering (CNF)
Please refer to Figure 3c. Due to the presence of noise, errors will occur on the cluster numbers. Therefore, BM3D can be used after the clustering is completed. The cluster number distribution map is used as the input noise image of the BM3D algorithm, and the filtered cluster numbers m f iltered (p) are used for the following other steps. This method is abbreviated as the CNF method in this paper.

Unwrapped Phase Filtering (UPF)
Please refer to Figure 3d. After the CA algorithm is completed, the unwrapped phase is restored. Due to the presence of noise, the obtained unwrapped phase contains errors. Therefore, BM3D can be used on the unwrapped phase, and the unwrapped phase is used as the input of the BM3D algorithm. Additionally, the filtered unwrapped phase ψ f iltered (p) is obtained at last. This method is abbreviated as the UPF method in this paper.

Simultaneous Filtering (STF)
Please refer to Figure 3e. In fact, BM3D can be used to filter the interferogram, the intercept, the cluster number, and the unwrapped phase at the same time, so the best filtering results should be obtained. This method is abbreviated as the STF method in this paper.

Denormalization
After filtering, the obtained interferogram, intercept, cluster number, and unwrapped phase are, respectively, denormalized so that the subsequent clustering can be processed normally and the correct unwrapped phase can be obtained, that is, processed according to the following four equations: where ϕ f inal (p) represents the final wrapped phase, b f inal (p) represents the final intercept, m f inal (p) represents the final cluster number, and ψ f inal (p) represents the final unwrapped phase.

Post-Processing
After being processed by the IFF method, the data obtained by filtering the interferogram need to be used in the subsequent steps of the CA-based MBPU methods to obtain the intercept, obtain the cluster number from the intercept, obtain the ambiguity number from the cluster number, and finally obtain the unwrapped phase. After being processed by the ICF method, the filtered intercept is used in the subsequent steps of the CA-based MBPU methods to obtain the cluster number, after which the ambiguity number is obtained from the cluster number and the unwrapped phase is finally obtained. After being processed by the CNF method, the obtained cluster number is used in the subsequent steps to obtain the ambiguity number and, finally the unwrapped phase. After being processed by the UPF method, the filtered unwrapped phase is obtained.

Results and Discussion
Experiments on synthesized and real InSAR datasets are shown in this section to verify the feasibility and effectiveness of the improved CA-based MBPU method by adding BM3D denoising.
The first experiment establishes a simple simulation scene, which shows smooth and discontinuous regions. The height of the scene was settled to 35 m and 80 m, respectively (the DEM is shown in Figure 4a), and the corresponding coherence coefficients were set to 0.5 and 0.6, respectively. The noiseless interferograms generated by simulation are shown in Figure 4b,c. There are only two clusters, and the corresponding cluster intercepts are 1 and 1/3, respectively. The real intercept distribution is shown in Figure 4d. The real cluster number distribution is shown in Figure 4e. The simulated noiseless unwrapped phase is shown in Figure 4f,g. The simulated noisy interferograms are shown in Figure 4h,i. As can be seen from the figure, the phase distribution becomes very ambiguous due to the addition of noise. The noisy intercept distribution is shown in Figure 4j. The noisy cluster number distribution obtained by the CA algorithm is shown in Figure 4k. The unwrapped phases obtained by the CA-based MBPU method are shown in Figure 4l,m.
Next, the BM3D algorithm was applied to the CA-based MBPU method to verify the effectiveness of the proposed method. First, the wrapped phases (interferograms) of long and short baselines with noise are taken as the inputs of BM3D and then filtered. The two filtered interferograms are shown in Figure 5a,b. They are used in the CA algorithm to obtain the cluster number distribution, as shown in Figure 5c. Second, BM3D is used to filter the noisy intercepts. The filtered intercept distribution map is shown in Figure 5d. Additionally, it is then applied to the CA algorithm to generate a cluster number distribution map, which is shown in Figure 5e. Third, the cluster numbers generated by clustering were filtered, as shown in Figure 5f. Fourth, the recovered unwrapped phases were filtered. The filtered unwrapped phases corresponding to the long and short baselines are shown in Figure 5g,h. Finally, BM3D is applied to the above four objects simultaneously, and the final unwrapped phases corresponding to the long and short baselines are shown in Figure 5i,j. It can be seen from the experimental results that the performance greatly improved after adding BM3D separately or simultaneously. are 1 and 1/3, respectively. The real intercept distribution is shown in Figure 4d. The real cluster number distribution is shown in Figure 4e. The simulated noiseless unwrapped phase is shown in Figure 4f,g. The simulated noisy interferograms are shown in Figure  4h,i. As can be seen from the figure, the phase distribution becomes very ambiguous due to the addition of noise. The noisy intercept distribution is shown in Figure 4j. The noisy cluster number distribution obtained by the CA algorithm is shown in Figure 4k. The unwrapped phases obtained by the CA-based MBPU method are shown in Figure 4l,m.  Next, the BM3D algorithm was applied to the CA-based MBPU method to verify the effectiveness of the proposed method. First, the wrapped phases (interferograms) of long and short baselines with noise are taken as the inputs of BM3D and then filtered. The two filtered interferograms are shown in Figure 5a,b. They are used in the CA algorithm to obtain the cluster number distribution, as shown in Figure 5c. Second, BM3D is used to (d) true intercept distribution; (e) true cluster number distribution; (f) simulated noiseless unwrapped phase for the long baseline (unit: rad); (g) simulated noiseless unwrapped phase for the short baseline (unit: rad); (h) simulated noisy interferogram for the long baseline (unit: rad); (i) simulated noisy interferogram for the short baseline (unit: rad); (j) noisy intercept distribution; (k) noisy cluster number distribution; (l) unwrapped phase without filtering for the long baseline (unit: rad); (m) unwrapped phase without filtering for the short baseline (unit: rad). tion map, which is shown in Figure 5e. Third, the cluster numbers generated by clustering were filtered, as shown in Figure 5f. Fourth, the recovered unwrapped phases were filtered. The filtered unwrapped phases corresponding to the long and short baselines are shown in Figure 5g,h. Finally, BM3D is applied to the above four objects simultaneously, and the final unwrapped phases corresponding to the long and short baselines are shown in Figure 5i,j. It can be seen from the experimental results that the performance greatly improved after adding BM3D separately or simultaneously.  Then, the performance was quantitatively analyzed as follows. We used the phase unwrapping success rate (PUSR) and normalized reconstruction square error (NRSE) to check the accuracy of the results. First, the correctness of the unwrapped result was judged according to the correctness of the final obtained ambiguity vectors corresponding to all pixels. PUSR is defined as the percentage of pixels whose integer ambiguity is correctly restored in an interferogram to the total pixels. The higher the success rate of phase unwrapping, the higher the accuracy of representing the ambiguity number. When the accuracy of the ambiguity number is more correct, it shows that the accuracy of the clustering results is higher. Therefore, the success rate of phase unwrapping can be used to reflect the accuracy of the clustering results. Table 1 shows the PUSR of various methods in the experiment. It can be seen from Table 1 that the PUSR has been greatly improved after filtering. Obviously, it can be seen from Table 1 that the filtering effect of STF is the best because the STF filters all of the four objects. However, if only one object is filtered, the filtering effect of UPF and CNF is better than that of IFF and ICF. Second, the accuracy of Then, the performance was quantitatively analyzed as follows. We used the phase unwrapping success rate (PUSR) and normalized reconstruction square error (NRSE) to check the accuracy of the results. First, the correctness of the unwrapped result was judged according to the correctness of the final obtained ambiguity vectors corresponding to all pixels. PUSR is defined as the percentage of pixels whose integer ambiguity is correctly restored in an interferogram to the total pixels. The higher the success rate of phase unwrapping, the higher the accuracy of representing the ambiguity number. When the accuracy of the ambiguity number is more correct, it shows that the accuracy of the clustering results is higher. Therefore, the success rate of phase unwrapping can be used to reflect the accuracy of the clustering results. Table 1 shows the PUSR of various methods in the experiment. It can be seen from Table 1 that the PUSR has been greatly improved after filtering. Obviously, it can be seen from Table 1 that the filtering effect of STF is the best because the STF filters all of the four objects. However, if only one object is filtered, the filtering effect of UPF and CNF is better than that of IFF and ICF. Second, the accuracy of the unwrapped result is judged through NRSE = ĥ − h 2 / h 2 , used in [11][12][13]21,25], whereĥ is the height array estimated from the unwrapped phase array of an interferogram and h is the reference height array. As shown in Table 2, the NRSE of before filtering, IFF, ICF, CNF, UPF, and STF is 0.1812, 0.0342, 0.0221, 0.089, 0.0080 and 0.0045, respectively. Therefore, IFF, ICF, CNF, UPF, and STF can effectively improve the PU accuracy. Table 3 shows the time consumption of the five filtering strategies. As can be seen from Tables 2 and 3, if the time consumption is not considered, STF can better reduce the normalized square error and improve the PU accuracy than the other filtering strategies. However, the STF method has the highest time complexity. For the other four filtering strategies, they all have similar time complexity, but the filtering effect of UPF and CNF is better than that of IFF and ICF.
The second experiment tested the performance of different filtering strategies on a more realistic height profile. Figure 6a shows the reference DEM of Isolation Peak in Colorado, USA. The corresponding ambiguity heights were set to 32.1 m and 53.5 m, respectively. The simulated two noiseless interferograms are shown in Figure 6b,c. The real intercept distribution is shown in Figure 6d. The real cluster number distribution is shown in Figure 6e. The noiseless unwrapped phases for the long and short baselines are shown in Figure 6f,g. The noisy interferograms are shown in Figure 6h,i. As can be seen from the figure, the phase distribution became very ambiguous due to the noise. The noisy intercept distribution is shown in Figure 6j. The noisy cluster number distribution obtained by the CA algorithm is shown in Figure 6k. The noisy unwrapped phases for the long and short baselines obtained by the CA-based MBPU method are shown in Figure 6l,m.  Figure 7 shows the experimental results about Isolation Peak (CO, USA) with BM3D denoising. The two filtered interferograms of IFF are shown in Figure 7a,b, and the corresponding cluster number distribution is shown in Figure 7c. The filtered intercept distribution of ICF is shown in Figure 7d, and the corresponding cluster number distribution is shown in Figure 7e. The filtered cluster number distribution of CNF is shown in Figure 7f. The filtered unwrapped phases of UPF corresponding to the long and short baselines are shown in Figure 7g,h. The final unwrapped phases of STF corresponding to the long and short baselines are shown in Figure 7i,j. The relevant quantitative indexes for Experiment 2 are shown in Tables 1-3, respectively.
In the last experiment, a small-scale double-baseline TanDEM-X InSAR dataset with 736 × 191 pixels was tested by the proposed method. The scenario of this dataset is a mountainous area in Tongchuan City, Shaanxi Province, China. Figure 8a is the reference DEM obtained by the Space Shuttle Radar Topographic Mapping Mission (SRTM). The main parameters of the TanDEM-X InSAR dataset are shown in Table 4. The two noisy TanDEM-X InSAR interferograms are shown in Figure 8b,c. The noisy intercept distribution is shown in Figure 8d. The noisy cluster number distribution obtained by the CA algorithm directly is shown in Figure 8e. The filtered interferograms of IFF are shown in Figure 8f,g. They are used to obtain the cluster number distribution. Figure 8h shows the cluster number distribution of IFF. The filtered intercept distribution of ICF is shown in Figure 8i, and the corresponding cluster number distribution of ICF is shown in Figure 8j. The filtered cluster number distribution of CNF is shown in Figure 8k. The filtered unwrapped phase of UPF corresponding to the long and short baselines are shown in Figure 8l   The second experiment tested the performance of different filtering strategies on a more realistic height profile. Figure 6a shows the reference DEM of Isolation Peak in Colorado, USA. The corresponding ambiguity heights were set to 32.1 m and 53.5 m, respectively. The simulated two noiseless interferograms are shown in Figure 6b,c. The real intercept distribution is shown in Figure 6d. The real cluster number distribution is shown in Figure 6e. The noiseless unwrapped phases for the long and short baselines are shown in Figure 6f,g. The noisy interferograms are shown in Figure 6h,i. As can be seen from the figure, the phase distribution became very ambiguous due to the noise. The noisy intercept distribution is shown in Figure 6j. The noisy cluster number distribution obtained by the CA algorithm is shown in Figure 6k. The noisy unwrapped phases for the long and short baselines obtained by the CA-based MBPU method are shown in Figure 6l,m.   responding cluster number distribution is shown in Figure 7c. The filtered intercep tribution of ICF is shown in Figure 7d, and the corresponding cluster number distrib is shown in Figure 7e. The filtered cluster number distribution of CNF is shown in F 7f. The filtered unwrapped phases of UPF corresponding to the long and short bas are shown in Figure 7g,h. The final unwrapped phases of STF corresponding to the and short baselines are shown in Figure 7i,j. The relevant quantitative indexes for E ment 2 are shown in Tables 1-3,   In the last experiment, a small-scale double-baseline TanDEM-X InSAR datase 736 × 191 pixels was tested by the proposed method. The scenario of this dataset is a m tainous area in Tongchuan City, Shaanxi Province, China. Figure 8a is the reference obtained by the Space Shuttle Radar Topographic Mapping Mission (SRTM). The parameters of the TanDEM-X InSAR dataset are shown in Table 4. The two noisy DEM-X InSAR interferograms are shown in Figure 8b,c. The noisy intercept distrib is shown in Figure 8d. The noisy cluster number distribution obtained by the CA Figure 7. Experiment about Isolation Peak (CO, USA) with BM3D denoising. (a) Filtered wrapped phase of IFF for the long baseline (unit: rad); (b) filtered wrapped phase of IFF for the short baseline (unit: rad); (c) cluster number distribution of IFF; (d) filtered intercept distribution of ICF; (e) cluster number distribution of ICF; (f) filtered cluster number distribution of CNF; (g) filtered unwrapped phase of UPF for the long baseline (unit: rad); (h) filtered unwrapped phase of UPF for the short baseline (unit: rad); (i) unwrapped phase of STF for the long baseline (unit: rad); (j) unwrapped phase of STF for the short baseline (unit: rad). filtered cluster number distribution of CNF is shown in Figure 8k. The filtered unwrapped phase of UPF corresponding to the long and short baselines are shown in Figure 8l,m. The unwrapped phases of STF corresponding to the long and short baselines are shown in Figure 8n,o. The relevant quantitative indexes for Experiment 3 are also shown in Tables 1-3, respectively. Obviously, the results of this experiment are in perfect agreement with the results of Experiment 1 and Experiment 2, thus strongly demonstrating the effectiveness and reliability of the filtering strategy proposed in this paper.  (d) intercepts of all the pixels obtained by the linear combination of (b,c); (e) cluster number distribution obtained by the CA algorithm directly; (f) filtered wrapped phase of IFF for the long baseline (unit: rad); (g) filtered wrapped phase of IFF for the short baseline (unit: rad); (h) cluster number distribution of IFF; (i) filtered intercept distribution of ICF; (j) cluster number distribution of ICF; (k) filtered cluster number distribution of CNF; (l) filtered unwrapped phase of UPF corresponding to the long baseline (unit: rad); (m) filtered unwrapped phase of UPF corresponding to the short baseline (unit: rad); (n) unwrapped phase of STF corresponding to the long baseline (unit: rad); and (o) unwrapped phase of STF corresponding to the short baseline (unit: rad). The three experimental results in Tables 1-3 demonstrate the effectiveness of the BM3D filtering strategies for the CA-based MBPU method on the simulated and real InSAR datasets.

Conclusions
MBPU is one of the key steps of InSAR processing. In order to improve the PU accuracy and the robustness of the CA-based MBPU method, five different filtering strategies for applying BM3D are proposed: interferogram filtering (IFF), intercept filtering (ICF), cluster number filtering (CNF), unwrapped phase filtering (UPF), and simultaneous filtering (STF). In particular, while keeping the general structure of BM3D, four different similarity measures are defined for interferograms, intercepts, clusters, and unwrapped phases to accommodate the special characteristics of different filtering objects. Experiments on simulated and real InSAR data prove the effectiveness and feasibility, and the experimental results show that (1) the PU accuracy and robustness of the CA-based MBPU method can be greatly improved by adding BM3D denoising; (2) simultaneous filtering of interferograms, intercepts, cluster numbers, and unwrapped phases works best, but with the worst time complexity; (3) when filtering is performed for only one object of the CA-based MBPU method, the filtering effect of CNF and UPF is better than that of IFF and ICF; and (4) considering the three indicators of PUSR, NRSE, and time consumption, CNF and UPF should be the best choice. The above conclusions can serve as a guideline for the future application of the BM3D algorithm in MBPU.

Conflicts of Interest:
The authors declare no conflict of interest.