A Novel Adaptive InSAR Phase Filtering Method Based on Complexity Factors

I.   Introduction

Synthetic aperture radar (SAR) is one of the most vital technology in the field of microwave remote sensing. Interferometric synthetic aperture radar (InSAR) is an important branch of SAR. It provides the digital elevation model (DEM) and earth surface deformation [1], [2], which have been successfully and widely used in civil and scientific research fields.

Interferometric phase noise cannot be avoided due to the existence of thermal noise and various decorrelation factors, which increase the difficulty of phase unwrapping and ultimately affects the precision of DEM and deformation reconstruction [3]. Thus, phase filtering is crucial for improving the quality of SAR interferograms before phase unwrapping [4]. The first and foremost objective of InSAR phase filtering is suppressing the noise as much as possible while preserving fringe details adaptively.

With the rapid development of InSAR, a large number of phase filtering methods have been proposed by researchers, which can be generally divided into spatial domain filtering and transform domain filtering.

In the spatial domain filtering, the simplest filtering methods are the mean filter and the median filter. Early spatial domain filtering methods adopted fixed filter window [5], so they are poor in fringe detail preservation. Lee et al. proposed a filtering method that adaptively adjusts the size and orientation of the filter window to maintain a balance between denoising and fringe detail preservation [6]. However, its limited 16 window orientations may cause distortion to complicated fringes. In view of the problems in Lee’s filtering method, several papers propose improvements to the Lee’s filter in terms of directional window [7], fringe frequency [8], and phase gradient [9]. The aforementioned filtering methods are based on interferometric phase itself. In 2006, Vasile et al. proposed an intensity driven adaptive neighborhood (IDAN) method [10], which detects homogeneous pixels in the neighborhood by SAR intensity information, to assist the interferometric phase filtering. However, its performance will be undermined when there are not enough homogeneous pixels in the neighborhood or heterogeneous pixels are included. After that, a complex Markov random field (CMRF) filter, which estimates the phase by minimizing the local energy function in the window, was proposed [11], [12]. To adapt to the changes of fringe pattern, Li et al. proposed a variable window CMRF filter [13].

The spatial domain filtering methods introduced above are all based on locally adjacent pixels. In the early days, due to the limitation of computing resources, the efficiency of an algorithm was the most concerned index. In recent years, with the improvement of computer technology, non-local phase filtering methods have been developing rapidly [14], [15], so as to overcome the constraint of estimating the phase in a local window. With phase similarity calculated in a matching window, weighted averaging of similar pixels is performed, therefore non-local phase filtering can reduce noise while preserving structures [16], [17]. Buades et al. first proposed a non-local mean filtering method [18]. Then, Deledalle et al. applied the non-local idea to interferometric phase filtering based on the statistical characteristics of InSAR data and proposed the non-local InSAR (NL-InSAR) filtering method [19]. The patch size is adaptively selected based on the heterogeneity of local scenes. Li et al. proposed an improved non-local filter that uses a normalized probability density function to measure the similarity between the center pixel and the remaining pixels in the matching window [20]. However, the interferometric fringes are dense in steep terrain areas, it is difficult to select similar pixels, which in turn diminishes noise reduction. Therefore, introducing the fringe frequency compensation technology to NL-InSAR can reduce fringe density and increase the number of similar pixels in the search window, which is conducive to obtaining more reliable noise suppression results.

The transform domain filtering method rely on transforming the interferometric phase from the spatial domain to the frequency domain or the wavelet domain for filtering, and then transforming filtered results back to the spatial domain. In 1998, Goldstein et al. proposed a classic frequency domain filtering method [21], which achieved noise reduction by smoothing the frequency spectrum. Besides the frequency domain, also the wavelet domain has been considered for phase filtering. A complex wavelet interferometric phase filter (WInPF) is implemented in [22]. The useful signals are extracted and amplified by utilizing discrete wavelet transform. In the wavelet domain, the phase information and noise are easier to separate [23], however, which are highly depend on the wavelet decomposition and wavelet coefficients [24]. Bioucas-Dias et al. proposed the phase estimation using adaptive regularization based on local smoothing (PEARLS) algorithm [25], which adaptively determines the polynomial fitting window size by the intersection of confidence intervals (ICI) algorithm.

In Goldstein’s method, noise is suppressed based on the different spectrum characteristics between useful signal and the noise. However, the filter parameter $\alpha$ takes a fixed value from [0, 1], and $\alpha=0$ means no filtering applied. With the increase of $\alpha$ , the smoothing effect is enhanced but the fringe details may be damaged [26]. Therefore, to balance the fringe details preservation and noise suppression, efforts have been made by researchers to improve the selection of $\alpha$ . Baran et al. associate $\alpha$ with the coherence, and at the same time define the smoothing operator as the convolution with the mean kernel function [27]. Li et al. proposed that the noise standard deviation can be used for adaptive selection of $\alpha$ [28]. These improved Goldstein filtering methods using coherence or noise standard deviation are dependent on coherence, which is usually a biased estimation. In order to reduce the impact of biased coherence on the filtered phase, Zhao et al. proposed a pseudo coherence to determine $\alpha$ [29]. In addition, Sun et al. proposed to use the signal-to-noise ratio defined by the noise variance as a criterion to selection $\alpha$ [30]. To overcome the subjectivity of the selection of $\alpha$ , Song et al. proposed two improved Goldstein filter methods based on empirical mode decomposition and adaptive neighborhood respectively [31], [32]. Generally, there tend to be more noise in the dense fringe area, which increases the difficulty of filtering. Therefore, a Goldstein filter method based on local fringe frequency compensation was proposed in [26]. However, fringe frequency estimation in both noisy phase and residual phase makes it quite time-consuming.

It can be seen that neither non-local filtering nor Goldstein filtering can perform well in dense fringes of steep terrain. Facing the tradeoff between noise suppression and fringe detail preservation, researchers combine classic filtering methods with the fringe frequency compensation [26], [33]. The essence of local fringe frequency (LFF) compensation is to estimate the LFF of the noisy phase [34], execute the filtering in residual phase after removing LFF, and add back the removed phase to the filtered residual phase. Trouve et al. proposed a filtering algorithm based on LFF estimation to enhance the fringe preserving ability [33], [35]. The limitation of this algorithm lies in its fixed window size and linearity presumption of LFF. Aiming at this problem, Cai et al. proposed that the size and shape of LFF estimation window should be adaptively determined by the coherent accumulation principle, which is rather time-consuming [36]. Suo et al. proposed a high-order LFF estimation method based on weighted least squares phase unwrapping [37]. This algorithm breaks the linearity presumption and behaves better in preserving phase details. However, its pre-filtering and phase unwrapping process may introduce additional operation and errors. From the above analysis, the improved methods for LFF exhibit high time complexity and cannot perform adaptive LFF estimation for fringe frequency compensation.

For areas with sparse fringes in the interferogram, LFF compensation is unnecessary for noise suppression. However, the dense fringe area needs additional LFF estimation operation, and the filtering methods based on LFF compensation are time-expensive. In addition, the main drawback of the existing LFF estimation methods is that the linear or nonlinear LFF cannot be estimated adaptively based on the terrain slope, resulting in large residual LFF estimation error or time consumption. On the other hand, there are heavy noises in low-correlation area, where needs stronger filtering strength. However, the high-correlation area needs to reduce the filter strength to prevent the fringe details from being damaged. For an interferogram, it is difficult to achieve high-performance noise suppression efficiently with a single phase filtering methods. The complicated and changeable terrain in the interferogram will pose a huge challenge on these non-adaptive phase filtering methods. Therefore, to deal with these problems, it is necessary to balance noise suppression, fringe detail preservation and computational efficiency based on different image characteristics. Factors like noise level and terrain slope information should be taken into account to select the optimal filtering strategy and filtering parameters for different areas.

To address the aforementioned problem, a novel adaptive InSAR phase filtering method based on complexity factors is proposed. Initially, the pseudo coherence, normalized maximum phase gradient (MPG), and normalized phase derivative variance (PDV) are calculated as the complexity factors. Based on the three complexity factors, the complexity indicator is constructed to guide the adaptive selection of the suitable and effective filtering strategy for different areas in the interferogram. Then, the complexity scalar, calculated by the three complexity factors, is used to guide the adaptive LFF estimation and adaptive filter parameters in different filtering methods. Experiments on simulated and real data prove that for complicated and changeable terrain, the proposed method can not only effectively suppress noise and preserve phase fringe details, but also increase the calculation efficiency.

The rest of this article is organized as follows. The novel adaptive phase filtering method based on complexity factors is presented in detail in Section II. The proposed method is tested on both simulated and real SAR data, where the experimental results are compared with those of the slope adaptive filtering, improved Goldstein filtering, and improved NL-InSAR filtering methods in Section III. The conclusions are drawn in Section IV.

II.   Adaptive InSAR Phase Filtering Method Based on Complexity Factors

Each filtering method has certain limitations and unique applicability. Usually, a specific filtering method is selected according to the characteristics of the entire image. However, different areas in the interferometric phase image usually exhibit different interferometric fringe patterns due to different terrain scene. For example, there are sparse and dense fringes in complicated and changeable terrain. For sparse fringe, most filtering methods can obtain satisfactory filtering result. So, the filtering methods with higher calculation efficiency is preferred. For the dense fringe, the filtering becomes more difficult, taking into account the fringe detail preservation and noise suppression. It is difficult for a single filtering strategy to be the optimal filtering for all areas when both filtering performance and calculation efficiency are required at the same time.

Therefore, a novel adaptive InSAR phase filtering method based on complexity factors is proposed. The flowchart is shown in Fig.1. Firstly, the three complexity factors related to the noise distribution and terrain slope are employed to adaptively select the filter strategies. Based on the three complexity factors, the complexity indicator ${\rm{CF}}1(m,n)$ of interferogram is constructed to guide the adaptive selection of the most suitable and effective filtering strategy for different areas of the entire image. Then, the complexity scalar ${\rm{CF}}2(m,n)$ is calculated according to the three complexity factors, which can guide the adaptive LFF estimation and compensation. The adaptive filter parameters are determined by ${\rm{CF}}2(m,n)$ of filter window in different filter methods for both residual phase and noisy phase without removing LFF. In summary, the proposed method adaptively selects the filtering strategy and adjusts the filtering parameters based on the noise and slope characteristics of the interferogram.

Figure  1.  Flowchart of the proposed method.

DownLoad: Full-Size Img PowerPoint

1   Complexity factors

Many indexes, including coherence, pseudo coherence, PDV, MPG, and second-order phase gradient can be used to describe the interferometric phase quality [38]. In this article, the pseudo coherence coefficient, normalized PDV and normalized MPG are selected as complexity factors, which represent the noise level and slope information.

1   Pseudo coherence coefficient

The pseudo coherence coefficient $\gamma(m,n)$ represents the noise level of interferometric phase and can be calculated with interferometric phase alone. The calculation formula is expressed as [39]

where $\gamma(m,n)$ is calculated in a $k\times k$ window centered on $(m,n)$ . $\varphi_{i, j}$ represents the interferometric phase on pixel $(i,j)$ in the local window.

2   Normalized PDV

The PDV of an interferogram is defined as [38]

where $\varDelta _{i,j}^{x} = {\rm wrap}({{\varphi }_{i+1,j}} - {{\varphi }_{i,j}})$ and $\varDelta _{i,j}^{y} = {\rm wrap} ({{\varphi }_{i,j+1}} - {{\varphi }_{i,j}})$ are the phase derivatives in row and column direction on $(i,j)$ respectively, $\bar{\varDelta}_{m, n}^{x}$ and $\bar{\varDelta}_{m, n}^{y}$ are the local mean values of the phase derivatives for pixel $(m,n)$ in row and column respectively.

Since $\gamma (m,n)\in [0,1]$ , in order to ensure the consistency for all complexity factors, the ${{\rm{PDV}}}'(m,n)$ is normalized to ${\rm{PDV}}(m,n)$ by

where $i \in M,\; j \in N$ , $M$ and $N$ represent the number of rows and columns of the interferogram respectively. And $\min \left({\rm{P D V}}^{\prime}(i, j)\right)$ and $\max \left({\rm{P D V}}^{\prime}(i, j)\right)$ represent the minimum and maximum PDV of the interferogram respectively.

3   Normalized MPG

The MPG of an interferogram is defined as [39]

where ${{\rm{MPG}}}'(m,n)$ is the maximum phase gradient of row and column direction in a $k\times k$ window centered on $(m,n)$ .

Similarly, ${{\rm{MPG}}}'(m,n)$ is normalized to ${\rm{MPG}}(m,n)$ by

where $\min \left({\rm{M P G}}^{\prime}(i, j)\right)$ and $\max \left({\rm{M P G}}^{\prime}(i, j)\right)$ represent the minimum and maximum MPG of the interferogram respectively.

It is known that the pseudo coherence coefficient $\gamma (m,n)$ is inversely proportional to the noise level. The ${\rm{PDV}}(m,n)$ and ${\rm{MPG}}(m,n)$ are related to the slope of the terrain, and their values are proportional to the steepness of the terrain. So pseudo coherence coefficient, normalized PDV, and normalized MPG represent the noise level and terrain slope information and are selected as complexity factors to guide the adaptive phase filtering of the interferogram.

2   Adaptive selection of filtering strategy

For interferograms with complicated and changeable terrain background, the noise and terrain slope changes greatly. It is necessary to select different filtering strategies for different areas with different terrain characteristics. Hence, an adaptive selection of filtering strategy is proposed based on the complexity factors.

$\gamma (m,n)$ characterizes the phase noise. The pseudo coherence coefficient indicator ${{\gamma }_{1}}(m,n)$ is calculated by

where ${{\gamma }_{{\rm{mean}}}}$ is the mean value of pseudo coherence coefficient of the interferogram. ${{\gamma }_{1}}(m,n)=0$ means that the phase noise level of the pixel $(m,n)$ is relatively low.

Similarly, PDV indicator ${{{\rm{PDV}}}_{1}}(m,n)$ and MPG indicator ${{{\rm{MPG}}}_{1}}(m,n)$ of pixel $(m,n)$ are calculated by

where ${{{\rm{PDV}}}_{{\rm{mean}}}}$ and ${{{\rm{MPG}}}_{{\rm{mean}}}}$ are the mean value of the normalized PDV and normalized MPG of an interferogram. ${{{\rm{PDV}}}_{1}}(m,n)=1$ and ${{{\rm{MPG}}}_{1}}(m,n)=1$ mean that the slope of the terrain near the pixel $(m,n)$ is relatively steep.

According to the relationship between the three complexity factors $\gamma (m,n)$ , ${\rm{PDV}}(m,n)$ , and ${\rm{MPG}}(m,n)$ of each pixel $(m,n)$ and the corresponding mean value ${{\gamma }_{{\rm{mean}}}}$ , ${{{\rm{PDV}}}_{{\rm{mean}}}}$ , and ${{{\rm{MPG}}}_{{\rm{mean}}}}$ , ${{\gamma }_{1}}(m,n)$ , ${{{\rm{PDV}}}_{1}}(m,n)$ , and ${{{\rm{MPG}}}_{1}}(m,n)$ are calculated, which are binary parameters of 0 or 1 with no dimension.

The complexity indicator ${\rm{CF}}1(m,n)$ can be calculated by

Only when ${{{\rm{PDV}}}_{1}}(m,n)$ and ${{{\rm{MPG}}}_{1}}(m,n)$ both equal to 1, the value of ${{{\rm{PDV}}}_{1}}(m,n) \times {{{\rm{MPG}}}_{1}}(m,n)$ is 1, which makes the steepness indicator more robust and reliable. ${{\gamma }_{1}}(m,n)=1$ means that the noise of the pixel $(m,n)$ is relatively large. Therefore, the three discrete values 0, 1, and 2 of ${\rm{CF}}1(m,n)$ can be obtained by the formula (9), which guides the adaptive selection of the filter strategy.

Moreover, the window size is very important for the filtering. Here, the adaptive window size is determined by the three complexity factors. Firstly, it is important to set a basic filter window size ${w_{\rm{indow}}}$ , then the filter window size of different filter strategies is adjusted based on ${w_{\rm{indow}}}$ and ${\rm{CF}}1(m,n)$ . ${w_{\rm{indow}}}$ is calculated by

where

${{{\rm{{\rm{PDV}}}}}_{{\rm{std}}}}$ and ${{{\rm{MPG}}}_{{\rm{std}}}}$ are the standard deviation of the normalized PDV and normalized MPG of the interferogram respectively. $\left\lceil \cdot \right\rceil$ rounds up its arguments to the nearest integer.

The ${\rm{ PDV}}\_r$ and ${\rm{MPG}}\_r$ are the maximum value of Z-score normalized ${\rm{PDV}}$ and Z-score normalized ${\rm{MPG}}$ respectively. Z-score normalization [40] could convert ${\rm{PDV}}$ and ${\rm{MPG}}$ to the same magnitude. The ${\rm{PDV}}\_r$ and ${\rm{MPG}}\_r$ are related with the largest terrain slope pixel of the interferogram. Sun et al. state that larger terrain slopes lead to more severe baseline decorrelation and lower signal noise ratio (SNR) in the interferogram [41]. The larger window size should be set for the lower SNR areas. Therefore, the ${\rm{PDV}}\_r$ and ${\rm{MPG}}\_r$ , determined by the largest terrain slope, are used to jointly set the basic radius of filter window size.

Correlation is inversely proportional to the noise level, therefore $\gamma \_r$ , determined by correlation coefficient, is used to adjust the filter window size. Here, 0.8 and 0.4 are chosen experimentally [1]. The quality of the interferimetric phase is good if $0.8 \lt {{\gamma }_{{\rm{mean}}}} \leq 1$ . On the contrary, $0 \leq {{\gamma }_{{\rm{mean}}}} \leq 0.4$ means that the noise is heavy and the filter window size needs to be increased to suppress noise. With the decrease of ${{\gamma }_{{\rm{mean}}}}$ , the filter window size will increase accordingly, to enhance noise suppression.

According to formula (9), ${\rm{CF}}1(m,n)$ has three discrete values, 0, 1, and 2, which represents the filtering difficulty of filter window centered on pixel $(m,n)$ . Therefore, different filtering strategies and different filtering window size centered on $(m,n)$ are chosen for each pixel according to ${\rm{ CF}}1(m,n)$ .

1) For ${{\rm{CF}}1(m,n)=1}$ , Goldstein filter based on LFF compensation

${\rm{CF}}1(m,n)=1$ means that the noise level and terrain slope in the pixel $(m,n)$ are moderate. By tradeoff between noise suppression and computation burden, the LFF compensation is needed, so the Goldstein filter method based on LFF compensation is adopted, and the LFF estimation window size is equal to filter window size according to formula (10). Considering that the noise and slope information are both moderate in this case, the filter window size of the Goldstein filter method based on LFF compensation equals to the basic filter window size ${w_{{\rm{indow}}}}$ , which means

2) For ${{\rm{CF}}1(m,n)=2}$ , NL-InSAR filter based on LFF compensation

${\rm{CF}}1(m,n)=2$ means that the noise level and terrain slope in the pixel $(m,n)$ are relatively higher. Usually, this means dense fringes and heavy noise, which increases the filtering difficulty compared with the first case ${\rm{CF}}1(m,n)=1$ . In this case, an NL-InSAR filtering method based on LFF compensation is adopted, which has great fringe detail preservation and noise suppression capabilities while heavy calculation burden.

For the NL-InSAR filtering method based on LFF compensation, the LFF estimation window size is also equal to ${w_{\rm{indow}}}$ . To ensure the fringe continuity of adjacent pixel, the filter window size span of different filter methods should not be too large, otherwise it will bring additional LFF estimation error. In the case ${\rm{CF}}1(m, n) =2$ , due to the fringe is denser than that in the case ${\rm{CF}}1(m,n)=1$ , to ensure the accuracy of LFF estimation, the filter window size needs to be slightly reduced. Therefore, the filter window size is calculated by

3) For ${{\rm{CF}}1(m,n)=0}$ , Goldstein filter

${\rm{CF}}1(m,n)=0$ means that the noise and terrain slope are lower. Compared with ${\rm{CF}}1(m,n)=2$ and ${\rm{CF}}1(m,n)=1$ , the filtering difficulty is greatly reduced. In this case, most filtering methods can achieve satisfactory noise suppression effects while preserving the fringe details. Therefore, the Goldstein filtering method is chosen for its high calculation efficiency, and the filter window size can be slightly increased. So, the filter window size is

3   Adaptive LFF estimation and compensation

In the case of ${\rm{CF}}1(m,n)=2$ and ${\rm{CF}}1(m,n)=1$ , the LFF is estimated and removed before filtering, so the noise can be effectively filtered out with little loss of fringes detail. In the traditional LFF estimation method [33], the size of the estimation window is fixed and LFF is assumed to be linear. However, for complicated and changeable terrain, LFF is generally nonlinear, in which case linear LFF compensation may cause terrain information loss during filtering. An adaptive LFF estimation and compensation method is proposed according to the complexity factors. To deal with the adverse impact of phase noise on LFF, prefiltering is employed before LFF estimation, then nonlinear or linear LFF estimation is performed adaptively to the terrain.

The complexity scalar ${\rm{CF}}2(m,n)$ , which reflects the noise level and terrain slope information, is thus defined as

where

where ${{\gamma }_{2}}(m,n)$ , ${{{\rm{PDV}}}_{2}}(m,n)$ and ${{{\rm{MPG}}}_{2}}(m,n)$ represent the mean value of $\gamma (i,j)$ , ${\rm{PDV}}(i,j)$ , and ${\rm{MPG}}(i,j)$ in the $k\times k$ filter window centered on pixel $(m,n)$ respectively. And the values of ${{\gamma }_{2}}(m,n)$ , ${{{\rm{PDV}}}_{2}}(m,n)$ , and ${{{\rm{MPG}}}_{2}}(m,n)$ are within [0, 1]. Therefore, ${\rm{CF}}2(m,n)$ is a continuous value within [0, 1], which is used to instruct the adaptive LFF estimation and calculate the filter parameters.

1   Prefilter with variable window size before LFF estimation

When ${\rm{CF}}1(m,n)=2$ and ${\rm{CF}}1(m,n)=1$ , it is known that the fringes are relatively dense and LFF compensation needs to be performed in the filter window. To improve the accuracy of LFF estimation, adaptive mean filter is implemented as prefilter, whose window size is determined by ${\rm{CF}}2(m,n)$ as

Since prefilter is performed within the filter window ${\rm{mid}}\_{\rm{w}}$ or ${\rm{min}}\_{\rm{w}}$ , therefore, it is necessary to ensure that the prefilter window size ${\rm{prefilter}}\_{\rm{win}}$ is smaller than ${\rm{mid}}\_{\rm{w}}$ or ${\rm{min}}\_{\rm{w}}$ . Also ${\rm{prefilter}}\_{\rm{win}}$ cannot be too large, otherwise it will reduce the estimation sensitivity. In addition, ${\rm{prefilter}}\_{\rm{win}}$ should be small to ensure the continuity of the LFF estimation value [26]. So ${\rm{prefilter}}\_{\rm{win }}$ varies from $3\times 3$ to $7\times 7$ , which is small enough for most interferogram.

${\rm{CF}}2(m,n)$ is proportional to the noise level and terrain slope in the filter window. Conversely, correlation coefficient is inversely proportional to the noise level. The boundary values 0.2 and 0.6 of ${\rm{CF}}2(m,n)$ are set according to the boundary values 0.8 and 0.4 of correlation coefficient in formula (13) respectively. Therefore, the areas of heavy noise will be prefiltered with larger windows.

It should be noted that prefilter with variable window size will improve the LFF estimation accuracy without losing details of the interferogram as it is not used for phase filter but only for LFF estimation.

2   Adaptive prominent fringe components estimation

The fringes are often nonlinear in complicated terrain. If only the linear fringe is removed, the residual fringe will hamper the phase filtering. Several linear or nonlinear LFF estimation methods have been proposed over the past years. However, the linear LFF estimation methods cannot accurately compensate the fringe frequency in areas with complicated terrain, which limits the ability of noise suppression in dense fringe areas. The nonlinear LFF estimation method [37] breaks through the first-order limitation of fringe frequency by performing weighted least squares phase unwrapping, which may affect the filtering result. Here, prominent fringe components estimation method is presented and the linear or nonlinear fringe is compensated adaptively.

The prominent fringe component is estimated by extracting the prominent frequency components of phase in the frequency domain. Implement 2-D FFT in local window to obtain the interferometric phase spectrum as follows:

where ${{I}_{n}}=\exp \left( {\rm{j}}{{\varphi }_{n}} \right)$ is the complex form of the noisy phase. The components with the amplitude of the spectrum greater than the threshold are retained, and lower than the threshold is set to zero.

where $b \in \left( 0,\max\left( |S\left( u,v \right)| \right) \right]$ . $b = \max \left( |S\left( u,v \right)| \right)$ means extracting the linear phase. The smaller the value of $b$ , the richer the nonlinear phase details. But if the value of $b$ is too small, the extracted phase will contain part of the phase noise. Therefore, it is necessary to select a reasonable threshold to effectively distinguish the prominent fringe spectrum from the noise spectrum.

It is proposed that $b=\max (|S(u, v)|) \times97\%$ is beneficial to estimate the nonlinear fringes more precisely in [42]. When ${\rm{CF}}1(m,n)=2$ and ${\rm{CF}}1(m,n)=1$ , the interferometric fringes are complicated and changeable, therefore, in our methods, the threshold is set to

Sorting the spectrum amplitude in descending order, the value of $X$ is determined by the complexity of the terrain, which is more adaptable to different scenarios. When $0 \leq {\rm{CF}}2(m,n) \leq 0.2$ , the prominent fringe is approximately equal to the linear fringe, so we set $X=1$ , that is $b=\max (|S(u, v)|) \times99\%$ , so as to avoid the noise spectrum. In the case of more complicated terrain, $0.6 \lt {\rm{CF}}2(m,n) \leq 1$ , the nonlinear phase details are richer, so the prominent fringe can be better extracted by taking $X=3$ , which means $b=\max (|S(u, v)|) \times 97\%$ .

The complex form of the prominent phase component in the local window is shown as

The residual phase ${{\varphi }_{r}}$ is obtained by removing the prominent phase ${{\varphi }_{m}}$ from the noisy phase ${{\varphi }_{n}}$ , which can be expressed as

When ${\rm{ CF}}1(m,n)\;=\;2$ or ${\rm{CF}}1(m,n)\;=\;1$ , through the above steps, the adaptive LFF compensation are performed in the filter window centered on pixel $(m,n)$ . Finally, the filtered interferometric phase ${{\varphi }_{f}}$ consists of the removed phase ${{\varphi }_{m}}$ and the filtered residual phase ${{\varphi }_{rf}}$ .

4   Adaptive filter parameters

Complexity scalar ${\rm{CF}}2(m,n)$ is also used to calculate the adaptive filter parameters in the different filtering methods.

1   For ${{\rm{CF}}1(m,n)=2}$ , residual phase NL-InSAR filtering based on LFF compensation

For ${\rm{CF}}1(m,n)=2$ , NL-InSAR filtering is performed on the residual phase. In the NL-InSAR filter method [19], the weights are determined by the similarity between the neighboring pixel blocks, and the calculation is related to the Euclidean distance of neighboring pixel blocks and smoothing parameter $h(m,n)$ . The filtering effect of NL-InSAR filter method is greatly affected by the smoothing parameter $h(m,n)$ , which is closely related to the decay speed of the weight function.

Usually a smaller $h(m,n)$ is beneficial to the preservation of fringes details, and a larger $h(m,n)$ can improve the noise suppression ability. In [42], the smoothing parameter $h(m,n)=10 \sigma_{n} \cdot \gamma \cdot(1+{f}'{x}^{2}+ {f}'{ y}^{2})^{-1 / 2}$ , which is jointly defined by the noise standard deviation $\sigma_{n}$ , the correlation coefficient $\gamma$ and the residual phase fringe frequency $({f}'{x},{f}'{ y})$ . The value of $\gamma \cdot\left(1+{f}'{x}^{2}+ {f}'{ y}^{2}\right)^{-1 / 2}$ is within (0,1). Since the residual phase fringe frequency $({f}'{x},{f}'{ y})$ need to be estimated, this method is time-consuming, and may introduce fringe frequency estimation errors.

Therefore, to ensure the fringe detail preservation and noise suppression effect at the same time, the smoothing parameter $h(m,n)$ is adaptively calculated by

The larger ${\rm{CF}}2(m,n)$ means denser fringes and heavier noise. ${\rm{CF}}2(m,n)$ varies from 0 to 1, so the value range of $h^{\prime}(m,n)$ is within (0.7,1). Here the minimum value of $h^{\prime}(m,n)$ is set as 0.7 because the residual phase after LFF compensation is so sparse that the smoothing parameter $h(m,n)$ can be larger. For the areas with larger ${\rm{CF}}2(m,n)$ , $1-{\rm{CF}}2(m,n)$ is smaller, which is more conducive to preserving fringe details. Therefore, ${\rm{CF}}2(m,n)$ is used to adjust the smoothing parameter $h(m,n)$ to gain a tradeoff between fringe preservation and noise suppression.

2   For ${{\rm{CF}}1(m,n)=1}$ , residual phase Goldstein filtering based on LFF compensation

Goldstein filtering is performed on the residual phase after removing the LFF. Goldstein filter method [21] converts the interferometric phase from the spatial domain to the frequency domain, and then smooths the frequency spectrum. The filter parameter ${{\alpha }_{r}(m,n)}$ indicates the extent to which the spectrum is smoothed.

For Goldstein filter, the smoothing effect becomes more intense with the increase of the filter parameter ${{\alpha }_{r}(m,n)}$ , whereas the ability of fringe detail preservation reduces. In order to balance the details loss and noise suppression ability, the traditional Goldstein filter parameter is usually set to 0.5. Therefore, in our method, ${{\alpha }_{r}(m,n)}$ is adaptive to ${\rm{CF}}2(m,n)$ by

Given that the area with ${\rm{CF}}1(m,n)=1$ has moderate noise and slope, The LFF compensation is performed before residual phase filter and the density of fringes has decreased, so the filter parameter ${{\alpha }_{r}(m,n)}$ can be a little larger. $\alpha _{r}^{\min }=0.5$ is set to ensure the filter intensity, and $1-{\rm{CF}}2(m,n)$ is also used to adjust the fringe detail preservation ability. The filter parameter ${{\alpha }_{r}(m,n)}\in \left( {\alpha _{r}^{\min} },1 \right)$ , therefore the filter intensity can be adaptively adjusted according to ${\rm{CF}}2(m,n)$ .

3   For ${{\rm{CF}}1(m,n)=0}$ , noisy phase Goldstein filtering

For ${\rm{CF}}1(m,n)=0$ , the noise level and terrain slope of the area are low. To improve the calculation efficiency, Goldstein filter is adopted in noisy phase without removing LFF, its filter parameter ${{\alpha }_{n}(m,n)}$ is adaptively determined by

where $\gamma _{{\rm{mean}}}^{{\rm{win}}}(m,n)$ is the mean value of pseudo coherence coefficient in the filter window centered on pixel $(m,n)$ . The Goldstein filter parameter in [27] is only changed with the noise level. To make the adaptive filter parameter more robust, in our method, the filter parameter ${{\alpha }_{n}(m,n)}$ changes with the noise level and terrain slope simultaneously because of $\gamma _{mean}^{win}(m,n)\times (1 - {\rm{CF}}2(m,n))$ . In the case ${\rm{CF}}1(m,n)=0$ , $\gamma _{mean}^{win}(m,n)$ is larger and ${\rm{CF}}2(m,n)$ is smaller. The minimum value of ${{\alpha }_{n}(m,n)}$ is $1-\gamma _{mean}^{win}(m,n)$ , which can prevent over-filtering.

III.   Results and Analysis

In this section, to validate the proposed method, experiments are performed on both simulated and real interferograms, and results are compared with those of several recognized and representative methods. The effectiveness of the three proposed adaptive strategies based on the complexity factors is demonstrated by the first experiment with simulated data. In the second experiment, the superiority of the proposed method in terms of noise suppression and fringe detail preservation is verified by comparison with the three adaptive filter methods—The slope adaptive filter [33] and the improved Goldstein filter [27] can implement adaptive phase filter based on fringe frequency and correlation coefficient, respectively. The improved NL-InSAR filter [42] is one of the latest adaptive filter method proposed in 2021. Finally, the real data are processed to further verify the robustness and superiority of the proposed method in comparison with three existing methods.

1   The first experiment to validate the adaptive strategies

Two SAR single-look complex images are simulated according to certain SAR geometry and DEM data [43]. The noise flattened phase and corresponding real phase are shown in Fig.2. In order to verify the effectiveness of the three adaptive strategies in the proposed method, the three adaptive strategies are sequentially replaced with fixed strategies, results are shown in the Fig.3.

Figure  2.  Simulated data 1.

DownLoad: Full-Size Img PowerPoint

Figure  3.  Filtered results using different adaptive strategies.

DownLoad: Full-Size Img PowerPoint

Clearly, the filter result of proposed method shown in Fig.3(d) contains more fringe detail than those in Fig.3(a), (b), and (c). The first fixed strategy is to adopt NL-InSAR filtering method based on LFF compensation for the interferogram, but the LFF estimation and filter parameters are still adaptive. It can be seen that without adaptive selection of filtering method, result in Fig.3(a) presents more phase residues. The second fixed strategy is that the LFF estimation window and the prefilter window are fixed to $8\times 8$ and $5\times 5$ respectively, $X=1$ is adopted to extract prominent fringe, the filter results is shown in Fig.3(b). The third fixed strategy is that the filter parameters $h(m,n)$ , ${{\alpha }_{r}(m,n)}$ , ${{\alpha }_{n}(m,n)}$ are fixed to 0.5, and Fig.3(c) shows the filter results. Without adaptive LFF estimation and adaptive filter parameters, the fringe detail in Fig.3(b) and Fig.3(c) are preserved not as good as in Fig.3(a) and Fig.3(d).

In order to evaluate the filtered results, the performance of each filter is assessed by the number of phase residues, the edge preservation index (EPI) [31] and the root mean square errors (RMSE) [26]. The EPI and RMSE are calculated by

where ${{\varphi }_{f}}(m,n)$ and ${{\varphi }_{{\rm{real}}}}(m,n)$ represent the filtered phase and the real phase respectively. $M\times N$ represents the total number of pixels in the interferogram.

An EPI closer to 1 means better fringe and edge preservation. ${\rm{EPI}} \gt 1$ means that fringe details are false or more noise on the filtered phase, whereas ${\rm{EPI}} \lt 1$ means the fringe details are damaged. Smaller RMSE represents a better noise suppression effect.

The evaluation results are shown in Table 1. It is obvious that the three adaptive strategies in the proposed method have improved fringe details preservation and noise suppression.

Table  1.  Evaluation results of simulated data1

Interferogram	Residues	EPI	RMSE (rad)
Real phase	1	1	0
Noisy flattened phase	3270	7.8684	1.1425
The first fixed strategy	4	1.0620	0.1078
The second fixed strategy	4	1.0995	0.1096
The third fixed strategy	3	1.1029	0.1068
Proposed method	1	1.0165	0.1011

 | Show Table

DownLoad: CSV

2   The second experiment to validate superiority

In this part, the superiority of the proposed method is evaluated on simulated data. The slope adaptive filter, the Improved Goldstein filter and the Improved NL-InSAR filter are implemented as comparison. The noisy flattened phase and corresponding real phase and are shown in Fig.4.

Figure  4.  Simulated data 2.

DownLoad: Full-Size Img PowerPoint

As can be seen in Fig.4, the fringe density is variable, which is sparse on the right side and dense on the left side. The terrain is relatively flat in sparse fringe areas and the terrain is more complicated in dense fringe area, which increases the difficulty of filtering.

In order to facilitate comparison, the window size of the other three filtering methods is same with that of ${\rm{CF}}1(m,n)=1$ in the proposed method, i.e., ${\rm{win}}\_{\rm{size}}= {\rm{mid}}\_{w}$ . In Fig.5, we show the corresponding filtered results of different methods. In each group, the left image is the filtered phase, and the right one is the corresponding phase error between the filtered and real phase.

Figure  5.  Filtered result and corresponding phase error.

DownLoad: Full-Size Img PowerPoint

As shown in Fig.5(a) and Fig.5(b), in the flat area on the right, the slope adaptive filter has a better noise suppression effect, but in the dense fringe area on the left, the fringe details are damaged, so the phase error is larger. In Fig.5(c) and Fig.5(d), although the overall phase error of the improved Goldstein filter is smaller, it can be clearly seen that the fringe edge preservation ability is poor, and the noise in some areas is still large. For the improved NL-InSAR filter, as it can be observed in Fig.5(e) and Fig.5(f), there are some error points at the edges of the image. Comparing the result of the proposed method in Fig.5(g) and Fig.5(h) with the above three filter methods, the proposed methods show a good performance in noise suppression and fringe details preservation. The phase error diagrams clearly show that the proposed method outperforms other filter methods.

The quantitative evaluation results are shown in Table 2. The slope adaptive filter has obvious over-filter phenomenon, and the EPI is far less than 1. Although its residue number is close to that of the proposed method, it is at the cost of fringe detail loss, resulting in a larger RMSE. The RMSE of the improved Goldstein filter has been reduced, but there are more residues and poor fringe detail preservation in areas with steep terrain and low coherence. For the improved NL-InSAR filter, the fringe details are damaged in the dense fringe area, resulting in residual points. Compared with the other three methods, the EPI of the proposed method is closer to 1, hence, the ability of fringe details preservation is much better. Moreover, the residues and RMSE of the proposed method are the smallest because of an excellent noise suppression performance, proving that the proposed method achieves the best balance between noise suppression and fringe preservation compared with the other three methods.

Table  2.  Evaluation results of simulated data 2

Interferogram	Residues	EPI	RMSE (rad)	Time (s)
Real phase	1	1	0	–
Noisy flattened phase	7258	7.9945	1.0807	–
Slope adaptive filter	1	0.9328	0.3020	36
Improved Goldstein filter	6	1.1004	0.2069	9
Improved NL-InSAR filter	2	1.0684	0.1946	45
Proposed method	0	1.0410	0.1750	30

 | Show Table

DownLoad: CSV

As shown in Fig.4(b), “A” represents the area where phase distortion often occurs. The cross-section of phase error in “A” is extracted to validate the robustness of the proposed method in filtering the steep terrain region. As clearly shown in Fig.6, the filtered phase error of the proposed method is much closer to zero than the other filter methods, which proves that the proposed method has a better performance on the edge preservation than the other three methods.

Figure  6.  Phase error of different filter methods in cross-section A.

DownLoad: Full-Size Img PowerPoint

3   Experiments with real data

In this part, two sets of real data are employed to investigate the performance of the proposed method.

1   ERS SAR data

ERS SAR images over the ETNA Volcano in September and October 2000 is used as the test data. The interferometric noisy phase and the enlarged area in the white rectangle are shown in Fig.7. It can be seen that the fringe in Fig.7(b), contaminated with heavy noise, represents the complicated terrain of ETNA Volcano. And the mean pseudo coherence coefficient is only 0.5152.

Figure  7.  ETNA Volcano.

DownLoad: Full-Size Img PowerPoint

The filtering results of Fig.7(d) with the slope adaptive filter, the improved Goldstein filter, the improved NL-InSAR filter, and the proposed method are shown in Fig.8. In Fig.8(a), due to the LFF estimated by the slope adaptive method is not accurate enough, resulting in damage of the fringes edges, which causes more residues. As can be seen in Fig.8(b) for the improved Goldstein filter, the fringes in the dense fringe area are ambiguous, especially in areas with a heavy noise. Comparing Fig.8(c) with Fig.8(a) and Fig.8(b), it is seen that the fringe preservation of the improved NL-InSAR filter is much better than the slope adaptive filter and the improved Goldstein filter in dense fringe areas. In Fig.8(d), the proposed method shows a better performance in fringe detail preservation and the fringe in steep terrain is the most continuous.

Figure  8.  Filtered results of Fig.7(b) with four methods.

DownLoad: Full-Size Img PowerPoint

A quantitative evaluation is also performed to compare the filtered results. The number of residues, the sum of phase difference (SPD) [44] and the phase standard deviation (PSD) [45] are employed as metrics. Compared with residues, the SPD and PSD can more accurately reflect the smoothness of the filtered phase. It is generally believed that smaller number of residues, SPD and PSD indicate a smoother phase with less noise.

The SPD for the interferogram is the sum of ${\rm{APD}}(m,n)$ , which is expressed by

where ${\rm{APD}}(m,n)$ is calculated as

where ${{\varphi }_{f}}(m,n)$ represents the filtered phase, and ${{\varphi }_{f}}(m+i,n+j)$ is the filtered phase of eight adjacent pixels.

The PSD of the interferogram is calculated by

where ${{\bar{\varphi }}_{f}}(m,n)$ is the linear phase in a window of size $3\times 3$ , and $M\times N$ represents the size of interferogram.

The evaluation results are shown in Table 3. As can be seen from Table 3, the number of residues are reduce by all methods. However, because of the fixed window used in slope adaptive filter, the loss of detail is severe in the dense fringe area, so its residues are more than the proposed method. The improved Goldstein filter and the improved NL-InSAR filter have the problem of under-filtering in the low-coherence area, resulting in a large SPD and PSD, and the residues of improved Goldstein filter is far more than that of the proposed method. Again, the proposed method shows a good improvement in terms of noise suppression, the number of residues, SPD and PSD are greatly reduced.

Table  3.  Evaluation results of real data1

Interferogram	Residues	SPD $\times \;10^4$ (rad)	PSD (rad)	Time (s)
Noisy phase	7455	6.3222	1.4953	–
Slope adaptive filter	19	1.8761	0.3273	40
Improved Goldstein filter	41	1.9855	0.4592	12
Improved NL-InSAR filter	13	1.9320	0.4289	51
Proposed method	4	1.8317	0.3109	25

 | Show Table

DownLoad: CSV

2   Millimeter-wave airborne InSAR data

The millimeter-wave airborne InSAR data are employed to conduct another experiment. The interferograms are provided by Beijing Institude of Radio Measurement, and the test site is situated in Zhaotong, Yunnan Province, Midwest China. The lower left corner of the interferogram is a residential area containing a lot of architectural details, and the upper right corner is a mountainous area, with partial shadows and layovers. The interferometric noisy phase and the area bounded by rectangle are shown in Fig.9.

Figure  9.  Millimeter-wave SAR data.

DownLoad: Full-Size Img PowerPoint

The slope adaptive filter, the improved Goldstein filter, the improved NL-InSAR filter and the proposed method are performed on Fig.9(b), and the filtered phase of all four methods are shown in Fig.10. The quantitative evaluation is given in Table 4.

Figure  10.  Filtered results of Fig.9(b) with four methods.

DownLoad: Full-Size Img PowerPoint

Table  4.  Evaluation results of real data 2

Interferogram	Residues	SPD $\times\;10^4$ (rad)	PSD (rad)	Time (s)
Noisy phase	6619	4.4042	1.4268	–
Slope adaptive filter	19	0.8478	0.3624	32
Improved Goldstein filter	20	0.7319	0.2900	12
Improved NL-InSAR filter	21	0.8113	0.3243	38
Proposed method	5	0.7286	0.2858	29

 | Show Table

DownLoad: CSV

Since the airborne data have relatively high SNR and sparse fringes, all four methods have achieved good noise suppression effect. From Fig.10 and Table 4, the fixed window in the slope adaptive filter causes the fringe details of the residential area to be completely filtered out. And the fringe edges are destroyed, resulting in residues. Compared with the slope adaptive filter, the SPD and PSD of the improved Goldstein filter and the improved NL-InSAR filter are reduced. However, in Figs.10(b) and (c), the dense fringes in the upper right corner where there is a mountainous shaded area, is severely damaged because the improved Goldstein filter and the improved NL-InSAR filter do not perform LFF compensation before filtering. The filtered result in Fig.10(d) and its quantitative evaluation in Table 4 show that the proposed method not only effectively suppresses noise, but also has the best performance in preserving fringe details.

It can be seen from the calculation efficiency in the tables, the filter speed of the improved Goldstein filter is always the fastest due to the frequency domain filter. The improved NL-InSAR filter has the slowest filter speed because twice LFF estimation is performed. The slope adaptive filter also requires LFF estimation for the entire interferogram. The proposed method adopts adaptive filter strategy, LFF estimation is not needed when $CF1(m,n)=0$ , which represents the noise level and fringe density are small. Therefore, the calculation efficiency is improved compared with the slope adaptive filter and the improved NL-InSAR filter.

IV.   Conclusions

With the increasing resolution of SAR imaging, there will be much more rich terrain types and detailed terrain information in the interferogram than ever. What is more, heavy noise caused by SAR imaging geometry and complicated terrain background brings a huge challenge to phase filtering. However, most of the current phase filtering methods cannot simultaneously take into account the three aspects of suppressing noise suppression effectively, preserving terrain details adaptively and improving calculation efficiency. Thus, a novel adaptive InSAR phase filtering method based on complexity factors is proposed in this paper. The complexity factors can characterize the noise level and terrain slope information of the interferogram effectively, which are used to guide the adaptive selection of suitable and effective filtering strategies for different areas in the interferogram.

The proposed method is presented in detail firstly, and after that, the proposed method is tested on simulated data and real data sets from ETNA Volcano and Yunnan Province, mountainous area in western China. By comparing its performance with the other three recognized and representative phase filtering method, it has been demonstrated that the proposed method offers the best filtering results. The adaptive selection of filtering strategy could improve calculation efficiency. Moreover, the adaptive LFF estimation and adaptive filter parameters based on complexity factors not only can effectively suppress noise, but also have excellent performance in preserving fringe details, the effectiveness and superiority of proposed method are validated.