Bias Removal for Goldstein Filtering Power Using a Second Kind Statistical Coherence Estimator

The adaptive Goldstein filter driven by InSAR coherence is one of the most famous frequency domain-based filters and has been widely used to improve the quality of InSAR measurement with different noise features. However, the filtering power is biased to varying degrees due to the biased coherence estimator and empirical modelling of the filtering power under a given coherence level. This leads to underor over-estimation of phase noise over the entire dataset. Here, the authors present a method to correct filtering power on the basis of the second kind statistical coherence estimator. In contrast with regular statistics, the new estimator has smaller bias and variance values, and therefore provides more accurate coherence observations. In addition, a piece-wise function model determined from the Monte Carlo simulation is used to compensate for the nonlinear relationship between the filtering parameter and coherence. This method was tested on both synthetic and real data sets and the results were compared against those derived from other state-of-the-art filters. The better performance of the new filter for edge preservation and residue reduction demonstrates the value of this method.


Introduction
Synthetic aperture radar interferometry (InSAR) has developed rapidly in the last three decades.Because of its accuracy, convenience, and efficiency, InSAR is an effective tool for surface deformation monitoring [1], forest management and classification [2,3], and ocean currents and earthquake inversion [4,5].Successful applications of this technique rely on the quality of the interferogram.However, SAR instruments suffer from both temporal and spatial decorrelation, leading to stochastic phase noises with different spatial features in InSAR measurements [6,7].The general consensus among previous studies is that phase noise likely causes error propagation in data processing (e.g., phase unwrapping) and increases the uncertainty of the final InSAR products [8].It is therefore imperative to provide effective filtering methods to enhance the quality of the interferogram and simultaneously preserve spatial information.
In recent years, many studies have tried to resolve the problem from spatial and frequency viewpoints.The Goldstein filter has proven to be a feasible method to improve the tradeoff between estimation accuracy and computational efficiency [7].It smoothens the Fourier spectrum of the interferogram by using the strong anti-noise capacity available in the frequency domain.Although this approach is promising, it was developed for geophysical applications that span large geographic regions.Therefore, the fringe pattern and spatial information cannot be preserved in detail in the interferogram, leading to blurred image features and loss of the phase signal.
In order to identify higher resolution surface changes, variants to overcome the shortcoming of classical Goldstein filtering have been developed based on the quality maps of interferometric signals [9][10][11][12].The key principle behind this is that strong noise areas are filtered more than low noise areas.The quality maps include coherence, theoretical phase standard deviation, the signal-to-noise ratio (SNR), and pseudo-correlation [9,10,13].A common limitation of these methods is that statistical estimators of these quality maps are usually highly biased and therefore under-estimate the phase noise in the real word.Although some of them mitigate the bias by modeling the bias with empirical functions, most of the models are determined using subjective criteria.To compensate for the impact of biased estimation, several methods also consider iterative filtering, especially over fast decorrelation areas [12][13][14].However, this behavior will likely cause artifactual fringe patterns because pure noise will contribute to the Fourier spectrum and filtering will confuse noise for fringes.
The primary objective of this paper was to develop, test, and validate a new method for InSAR interferogram filtering.This method resolves the issues mentioned above by accurately estimating the Goldstein filtering power.Specifically, (1) the weighted coherence estimator based on a second kind statistic is developed to calculate sample coherence.In contrast with regular statistics used in previous studies, this estimator has smaller bias and variance values, and can therefore infer better noise information for the interferogram; (2) a piece-wise function is established to quantify the filtering power using the Monte Carlo simulation.This model allows for heavy filtering over strong noise areas and preserves phase signals in areas with low noise.In the following sections, the authors outline the definition, describe the methods and results, and conclude with a discussion of lessons learned from this research.

Adaptive Goldstein Filtering
Goldstein and Werner first proposed an adaptive radar interferogram filter and defined it as the multiplication of the Fourier spectrum Z(u, v) of the interferogram patch by the smoothed modulus S{|Z(u, v)|} raised to the power of an exponent α: where H(u, v) denotes the filter response; S{•} denotes the smoothing operator; and u and v are spatial frequencies in each interferogram patch (e.g., 32 × 32 window).The filtering power α is used to control the smoothness of filtering, which increases with larger α values.According to the definition in [7], α ranges from 0 to 1.
Each interferogram patch has different noise features.The use of a uniform value of α will underor over-estimate the degree of noise, leading to significant changes to the structure of the interferogram.Therefore, adaptive behavior can be used to change α according to the quality of the interferogram [9].For example, α can be defined as a linear model: where parameter γ is used to evaluate the phase noise and usually measured by InSAR coherence [15,16].A higher value of γ implies a lower noise degree and vice versa.γ = 1 denotes a noise-free signal where the filter power α becomes zero and no filtering occurs.

Definition of the Coherence Estimator
One of the difficulties with Equation ( 2) is the fact that the estimator of coherence parameter γ is highly biased.This limitation further under-estimates the value of α, leading to under-filtering over low coherence regions.Therefore, it is necessary to obtain an accurate coherence estimation before filtering.The coherence estimator γ is defined as: where z 1 and z 2 are single complex images (SLC) used for interferogram formalization and n is the sample number in each estimate window.When weight w = 1, the equation converts to the conventional coherence estimator defined in [15].The effect of weight is to reduce the heterogeneity of neighboring pixels that need to be averaged.Therefore, according to the theory of WML estimation investigated in [17], γ is a weighted maximum likelihood (WML) estimator.

Weight Determination Weight
The study takes advantage of the similarity between the reference pixel and its neighbors to determine the weight for Equation (3).The similarity is designed from the two-sample Anderson-Darling test statistic by [18,19]: where F a (x), G b (x), and H N (x) with N = a + b are the empirical distribution functions (EDF) of a-samples, b-samples, and combined samples, respectively, and In this case, replacing a-samples and b-samples with m pixels in a patch centered at the reference pixel Θ re f and a neighbor Θ neig , respectively, to get a = b = m and N = 2m, and the test statistic, can be rewritten as: where F Θ re f (x) and G Θ neig (x) denote the EDF of reference and neighboring pixels, respectively, and are estimated from m samples in each patch, while H Θ re f ∪Θ neig (x) is the EDF of the combined samples in both patches.The Anderson-Darling (AD) test is a non-parametric test and a measurement of distance between two distributions [19].Looking at the numerator in Equation ( 5), the AD distance places more weight on observations towards the tails of the distribution.Thus, it is a better measure than its competitors, such as the Cramér-von Mises and Kolmogorov-Smirnov statistics.Given the fact that AD is proportional to the difference between two distributions, it can be used to measure the homogeneity between two pixels by using their samples in patches.
Let AD(i) be the value of the ith pixel in the coherence estimate window centered at a reference pixel P. Each AD(i) corresponding to P is estimated by sliding a window of size m (Figure 1a).This procedure is repeated until the estimation for the nth pixel is finished.Throughout this paper, a 5 × 5 window size is used to estimate AD values and m = 25.It is always less than that used to calculate coherence in Equation (3) and therefore n > m. Figure 1b presents footprints of similarity patches Θ neig (green) passed by the Anderson-Darling test under a significance level of α = 0.01, and all AD values corresponding to the reference pixel P have been shown in Figure 1c.Visually, pixels sharing similar backscattering features with the reference pixel (Figure 1b) have smaller AD values, while pixels with other textures have larger AD values., and all AD values corresponding to the reference pixel P have been shown in Figure 1c.Visually, pixels sharing similar backscattering features with the reference pixel (Figure 1b) have smaller AD values, while pixels with other textures have larger AD values.The weight by the reciprocal of the test statistic is defined as: ( ) It is worth pointing out that for patch ref neig  =  , ( ) 0 AD P = , and it makes no sense to apply weighting in coherence estimation.Instead, ( ) 0.1 AD P = is set according to 1000 simulated scenes with different textures.The method used to determine the optimal weight as 0.1 for the reference pixel is the same as in [20], where overall errors between true and estimated values from all possible weight values were explored for each scene and the weight with minimum error was selected.From a statistical viewpoint, there is no significance in estimating coherence by means of samples from different populations.Thus, the effect of weighting is to reduce spatial sample heterogeneity and to assign larger weights to those pixels that follow the same distribution as the reference pixel.Since SAR backscattering is usually used as a proxy for phase stability, the AD test is applied to intensity SAR [21].From the image processing viewpoint, same surface targets generally share similar backscattering properties, showing the same textures in SAR images.The selection of homogeneous pixels with respect to the reference pixel can preserve the image resolution and mitigate the edge-smearing phenomenon [22].However, because isolated targets may increase the variance of the estimated coherence image, the weighting strategy was developed to reduce the variance.

Bias Removal Using Second Kind Statistic Estimator
Nicolas and Anfinsen (2002) introduced second kind statistics based on the characteristic functions that substitute the Fourier transformation with the Mellin transformation [23].It has been demonstrated that the second kind variances are considerably lower than regular variances for probability density functions defined on R + and slightly higher than the Cramer-Rao bound [23].By deriving the second kind first characteristic function, the Mellin transformation allows the user to define the following second kind moments [23,24], The weight by the reciprocal of the test statistic is defined as: It is worth pointing out that for patch Θ re f = Θ neig , AD(P) = 0, and it makes no sense to apply weighting in coherence estimation.Instead, AD(P) = 0.1 is set according to 1000 simulated scenes with different textures.The method used to determine the optimal weight as 0.1 for the reference pixel is the same as in [20], where overall errors between true and estimated values from all possible weight values were explored for each scene and the weight with minimum error was selected.
From a statistical viewpoint, there is no significance in estimating coherence by means of samples from different populations.Thus, the effect of weighting is to reduce spatial sample heterogeneity and to assign larger weights to those pixels that follow the same distribution as the reference pixel.Since SAR backscattering is usually used as a proxy for phase stability, the AD test is applied to intensity SAR [21].From the image processing viewpoint, same surface targets generally share similar backscattering properties, showing the same textures in SAR images.The selection of homogeneous pixels with respect to the reference pixel can preserve the image resolution and mitigate the edge-smearing phenomenon [22].However, because isolated targets may increase the variance of the estimated coherence image, the weighting strategy was developed to reduce the variance.

Bias Removal Using Second Kind Statistic Estimator
Nicolas and Anfinsen (2002) introduced second kind statistics based on the characteristic functions that substitute the Fourier transformation with the Mellin transformation [23].It has been demonstrated that the second kind variances are considerably lower than regular variances for probability density functions defined on R + and slightly higher than the Cramer-Rao bound [23].By deriving the second kind first characteristic function, the Mellin transformation allows the user to define the following second kind moments [23,24], where m v denotes the second kind moment of the order v and p x (u) is the probability density function (PDF) of a specified random variable.Given PDF of sample coherence γ [24], where γ and n are the true coherence and the number of samples used to estimate sample coherence γ in Equation (3), respectively, and 2 F 1 (•) denotes the ordinary hypergeometric function.Let v = 1 and one can deduce the second kind mean by: and the equivalent expectation of second kind mean is then defined as: Figure 2 represents the expectation of sample coherence as a function of true coherence, in which the dashed line indicates the bias of the second kind statistic coherence estimator and is calculated by combining Equations ( 8)- (10) under each true coherence γ, while the dash-dotted line denotes the bias of the regular coherence estimator and is calculated by Equation ( 6) in [16].It can be seen that the sample coherence from the regular statistic is more biased than from the second kind statistic in Equation (10), especially over low coherence regions.However, the bias of the second kind estimator cannot be completely removed by inverting (10) because the analytical expressions for Equation ( 9) are not existent.If a sufficiently large number k of independent sample coherence γ are averaged: the second kind mean estimate γ tends to be distributed normally around m 1 , with the variance var( γ)/k.Therefore, a bias corrected estimate can be obtained at E( γ) ≈ e γ, γ : E( γ) can be substituted with its ML estimate e γ in Equation (10), which is inverted according to Equation (12) to deduce an unbiased coherence estimate γ.Tables of inversion are first calculated using Equations ( 8)- (10) to reduce the computational burden.Similar tables were obtained in [16] using simulations in the open source SHPS-InSAR toolbox that is available at http://mijiang.org.cn/index.php/software/.In reality, the sample mean in Equation ( 11) is calculated from the area in the coherence map covered by an effective patch in Goldstein filtering (patch minus the overlapping area: 32 × 32 window size with 28 overlapping pixels).Thus, k = 128 so that γ approximates to E( γ).
The performance of the unbiased estimate at k = 128 with its standard deviation is shown in purple in Figure 2.

Modeling Filtering Power Using Unbiased Coherence
Unbiased coherence can quantify the quality of the interferogram more accurately and the mitigation of heterogeneity in the coherence map further increases the spatial resolution of filtering.This ensures that the coherence map can guide Goldstein filtering correctly.However, the filtering

Modeling Filtering Power Using Unbiased Coherence
Unbiased coherence can quantify the quality of the interferogram more accurately and the mitigation of heterogeneity in the coherence map further increases the spatial resolution of filtering.This ensures that the coherence map can guide Goldstein filtering correctly.However, the filtering power over noise areas is still biased due to the unknown relationship between the filtering parameter used and the coherence magnitude.Therefore, the impact of filtering on different fringe patterns is investigated under each coherence level and the optimal filtering power is determined by minimizing the difference between the noise-free and filtered interferograms based on the Monte Carlo simulation.Specifically, the noise-free interferogram was first simulated using the universal multifractal technique with random multifractal parameter settings [25].The noise interferogram was then generated by adding phase noise with specified phase standard deviation, which is a function of coherence [26].Given 1000 noise interferograms with different patterns under a coherence level γ ∈ {0, 0.1, . . . ,1}, each of them was filtered by Goldstein filtering with discrete parameter (α) values from 0 to 1 and evaluated against the noise-free image using root mean square error (RMSE).The optimal filtering power α was finally determined by maximizing the number of times that α corresponding to minimum RMSE in each of 1000 trials was considered optimal.
Figure 3 presents the coherence level versus optimal filtering power.Compared with the empirical model (dashed line) of the Baran filter presented in previous studies [9], non-linear characteristics can be clearly observed in this case.This explains one of the reasons why under-filtering generally occurs in the Baran filter.To correct this bias, the piece-wise function model (solid line) was deployed to fit the trajectory by iterative reweighted least squares and can be expressed as: Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 15

Filter Development
The modification of bias correction for filtering power leads to a newly developed Goldstein filter, which can be presented schematically as follows: (1) For each image pixel P, a coherence estimate window of size n is defined.The similarity between P and its neighbors is compared individually on an average intensity SAR image with patch size m (Figure 1a), and the weight of the coherence estimate window is confirmed using Equation (6).The sample coherence for pixel P is then estimated using Equation ( 3).(2) The procedure is repeated until the sample coherence of the last pixel in the whole image is estimated and coherence map is generated.(3) Goldstein filtering is performed on the interferogram in which the coherence sample with size k in each filtering patch is first averaged to obtain  using Equation (11), and is further corrected to unbiased coherence  according to Equation (12).

Filter Development
The modification of bias correction for filtering power leads to a newly developed Goldstein filter, which can be presented schematically as follows: (1) For each image pixel P, a coherence estimate window of size n is defined.The similarity between P and its neighbors is compared individually on an average intensity SAR image with patch size m (Figure 1a), and the weight of the coherence estimate window is confirmed using Equation ( 6).The sample coherence for pixel P is then estimated using Equation ( 3).
(2) The procedure is repeated until the sample coherence of the last pixel in the whole image is estimated and coherence map is generated.(3) Goldstein filtering is performed on the interferogram in which the coherence sample with size k in each filtering patch is first averaged to obtain γ using Equation (11), and is further corrected to unbiased coherence γ according to Equation ( 12).(4) The filtering power is adjusted according to Equation ( 13) and Goldstein filtering is performed on such patches.( 5) Steps ( 3) and ( 4) are repeated until the whole interferogram has been filtered.
For this paper, the coherence window size in step ( 1) is set to 15 × 15 (n = 225) and the similarity patch is 5 × 5 (m = 25).In steps ( 3) and ( 4), the filtering window is set to be 32 × 32 with the overlaps between windows set at 28 pixels, leading to the averaged coherence sample number k = 128.The flowchart of the proposed algorithm is shown in Figure 4.

Synthetic Data
To investigate the effectiveness of the newly developed filter, the authors compared the results from both synthetic and real data against those of state-of-the-art filters, including the Baran filter [9], the adaptive Goldstein filter derived from unbiased coherence under regular statistics (AGFC) [11], and iterative pseudo-correlation (AGFP) [13].Without loss of generality, the classical window size 77  was used for AGFC and AGFP to estimate map quality.
The simulation starts from bi-temporal SLC imagery 1 z and 2 z , and proceeds as follows: (1)   generate image patterns with different spatial distributions.Given the independence between intensity and phase pattern, a panchromatic Gaofen-2 image with 400 400  pixels over the Shandong area was used as a noise-free SAR intensity images ( I ) (Figure 5a), and the coherence map (  ) was simulated by normalizing the Gaofen-2 image in the interval [0,1] (Figure 5b).The noise-free interferometric phase  was simulated using the universal multifractal technique (Figure 5c); (2)   generate noise image according to a covariance matrix: where 1 x and 2 x are independent complex Gaussian random variables with zero mean and unit *

Synthetic Data
To investigate the effectiveness of the newly developed filter, the authors compared the results from both synthetic and real data against those of state-of-the-art filters, including the Baran filter [9], the adaptive Goldstein filter derived from unbiased coherence under regular statistics (AGFC) [11], and iterative pseudo-correlation (AGFP) [13].Without loss of generality, the classical window size 7 × 7 was used for AGFC and AGFP to estimate map quality.
The simulation starts from bi-temporal SLC imagery z 1 and z 2 , and proceeds as follows: (1) generate image patterns with different spatial distributions.Given the independence between intensity and phase pattern, a panchromatic Gaofen-2 image with 400 × 400 pixels over the Shandong area was used as a noise-free SAR intensity images (I) (Figure 5a), and the coherence map (γ) was simulated by normalizing the Gaofen-2 image in the interval [0,1] (Figure 5b).The noise-free interferometric phase φ was simulated using the universal multifractal technique (Figure 5c); (2) generate noise image according to a covariance matrix: where x 1 and x 2 are independent complex Gaussian random variables with zero mean and unit variance.The corresponding noise interferogram after the conjugation operator z 1 z * 2 is shown in Figure 5d.A result of the filtering process is shown in Figure 6.Clearly, the filtering results strongly depend on estimated parameters.For example, there is no significant difference between the Baran filter (i) and AGFC (j), because they use the same linear model and similar coherence magnitudes (>0.3) to estimate the filtering power (e and f).However, larger variances of coherence can be observed in (a) and (b).On the contrary, the coherence estimated by the new method (d) better for both edge preservation and magnitude because of the reduction of heterogeneity and also due to bias removal.More accurate coherence estimation in conjunction with the piece-wise model also leads to enhanced contrast in the filtering power (h), and therefore significant noise mitigation in (l).Although the filtering power of AGFP in (g) is still under-estimated and does not follow the noise features in Figure 5d, the filtered interferogram in (k) is smoother than those estimated from other methods.This is expected since AGFP is an iterative filtering process that removes noise from the previous result.
Quantitatively, the RMSEs between each filtered interferogram and true values are counted to evaluate the performance of different filters.As expected, Figure 6i,j have the largest RMSEs, at 1.10 rad and 1.09 rad, respectively.The RMSE of Figure 6l is 0.49 rad, slightly larger than that filtered by AGFP in Figure 6k (0.46 rad) because the iterative filtering enhanced smoothness at the cost of lost resolution.To further confirm this, the authors performed 100 iterations for each method using Equation ( 14), in which the universal multifractal technique can generate phase patterns with varying details.The results are shown in Figure 7.It can be seen that the Baran filter and AGFC underestimate strong noise areas due to their biased filtering power and therefore have larger mean values for RMSEs.The iterative AGFP has a smaller mean but larger variance.The reason for this is that reduplicative filtering results in the loss of signal and consequent image blurring, and the degree of loss is related to the complexity of phase patterns in each simulation.In testing, the new filter performs better than all other filtering methods.A result of the filtering process is shown in Figure 6.Clearly, the filtering results strongly depend on estimated parameters.For example, there is no significant difference between the Baran filter (i) and AGFC (j), because they use the same linear model and similar coherence magnitudes (>0.3) to estimate the filtering power (e and f).However, larger variances of coherence can be observed in (a) and (b).On the contrary, the coherence estimated by the new method (d) looks better for both edge preservation and magnitude because of the reduction of heterogeneity and also due to bias removal.More accurate coherence estimation in conjunction with the piece-wise model also leads to enhanced contrast in the filtering power (h), and therefore significant noise mitigation in (l).Although the filtering power of AGFP in (g) is still under-estimated and does not follow the noise features in Figure 5d, the filtered interferogram in (k) is smoother than those estimated from other methods.This is expected since AGFP is an iterative filtering process that removes noise from the previous result.
Quantitatively, the RMSEs between each filtered interferogram and true values are counted to evaluate the performance of different filters.As expected, Figure 6i,j have the largest RMSEs, at 1.10 rad and 1.09 rad, respectively.The RMSE of Figure 6l is 0.49 rad, slightly larger than that filtered by AGFP in Figure 6k (0.46 rad) because the iterative filtering enhanced smoothness at the cost of lost resolution.To further confirm this, the authors performed 100 iterations for each method using Equation ( 14), in which the universal multifractal technique can generate phase patterns with varying details.The results are shown in Figure 7.It can be seen that the Baran filter and AGFC under-estimate strong noise areas due to their biased filtering power and therefore have larger mean values for RMSEs.The iterative AGFP has a smaller mean but larger variance.The reason for this is that reduplicative filtering results in the loss of signal and consequent image blurring, and the degree of loss is related to the complexity of phase patterns in each simulation.In testing, the new filter performs better than all other filtering methods.

Real Data
Two SAR subsets acquired from Sentinel-1A with TOPS mode and TanDEM-X bistatic strip map mode were used to evaluate the performance of the new filter over incoherent and coherent scenarios, respectively.The results are shown in Table 1.During pre-processing, the TOPS SAR images were co-registered, using precise satellite orbits and an external digital elevation model, followed by a

Real Data
Two SAR subsets acquired from Sentinel-1A with TOPS mode and TanDEM-X bistatic strip map mode were used to evaluate the performance of the new filter over incoherent and coherent scenarios, respectively.The results are shown in Table 1.During pre-processing, the TOPS SAR images were co-registered, using precise satellite orbits and an external digital elevation model, followed by a

Real Data
Two SAR subsets acquired from Sentinel-1A with TOPS mode and TanDEM-X bistatic strip map mode were used to evaluate the performance of the new filter over incoherent and coherent scenarios, respectively.The results are shown in Table 1.During pre-processing, the TOPS SAR images were co-registered, using precise satellite orbits and an external digital elevation model, followed by a refinement of shifts to estimate the azimuth offset between the two images.The TanDEM-X data was co-registered using a cross-correlation algorithm at a 0.1 pixel accuracy.Both datasets were compensated for topographic phase contribution before coherence estimation was conducted.co-registered using a cross-correlation algorithm at a 0.1 pixel accuracy.Both datasets were compensated for topographic phase contribution before coherence estimation was conducted.Compared with the Baran filter and AGFC in Figure 9a,b, the noise is significantly mitigated in interferograms from AGFP and for the new method in Figure 9c,d.This can also be confirmed from the sum of phase difference (SPD) [11] and the number of unwrapping residues in Table 2, in which a smaller SPD and residue value demonstrate the better performance of the new method for filtering for noise suppression [9,10].Despite having minimum values for AGFP, it seems to yield spurious fringe patterns on the right in Figure9c, where pure noise dominates.This issue is understandable as Compared with the Baran filter and AGFC in Figure 9a,b, the noise is significantly mitigated in interferograms from AGFP and for the new method in Figure 9c,d.This can also be confirmed from the sum of phase difference (SPD) [11] and the number of unwrapping residues in Table 2, in which a smaller SPD and residue value demonstrate the better performance of the new method for filtering for noise suppression [9,10].Despite having minimum values for AGFP, it seems to yield spurious fringe patterns on the right in Figure 9c, where pure noise dominates.This issue is understandable as pure noise will contribute to the Fourier spectrum due to the missing estimation of spatial frequencies.This becomes an important drawback for iterative methods that cannot accurately threshold the boundary between completely decorrelated and low coherence areas.Therefore, the authors conclude that although the noise may be further reduced by iteration, the result shows that the bias removal for Goldstein filtering power can avoid artifacts and simultaneously mitigate noise.Finally, the authors present the results over a coherent area to demonstrate the effectiveness of the new filter for edge preservation.To this end, a mountainous area with very high coherence was selected and is shown in Figure 10. Figure 11 shows the filtered interferograms using different methods, where the differences among filters in Figure 11a,b,d are not distinguishable, while the result from AGFP seems to over-smooth the interferogram.This is identical to the previous case of synthetic data.The authors then deployed the receiver operating characteristics (ROC) to quantify the impact of the filtering effect on edge preservation.To do this, the Canny operator was first used to binarize the image.The edge of the original interferogram (Figure 10a) was assumed to be the positive instance and its complement the negative instance.Then, all true positive and false positive instances were counted for binary images estimated from each interferogram.By thresholding the Canny detector with various magnitudes for each filtered interferogram, the ROC curves were constructed as shown in Figure 12.The more the ROC curves toward the upper left, the better the  Finally, the authors present the results over a coherent area to demonstrate the effectiveness of the new filter for edge preservation.To this end, a mountainous area with very high coherence was selected and is shown in Figure 10. Figure 11 shows the filtered interferograms using different methods, where the differences among filters in Figure 11a,b,d are not distinguishable, while the result from AGFP seems to over-smooth the interferogram.This is identical to the previous case of synthetic data.The authors then deployed the receiver operating characteristics (ROC) to quantify the impact of the filtering effect on edge preservation.To do this, the Canny operator was first used to binarize the image.The edge of the original interferogram (Figure 10a) was assumed to be the positive instance and its complement the negative instance.Then, all true positive and false positive instances were counted for binary images estimated from each interferogram.By thresholding the Canny detector with various magnitudes for each filtered interferogram, the ROC curves were constructed as shown in Figure 12.The more the ROC curves toward the upper left, the better the edge preservation.The curve estimated from ADFP is located in the bottom right corner, showing the worst performance for the quality of edge preservation.It is also important to note similar ROC features for other methods in the case.This is due to the fact that there is no bias on high coherence levels in Figure 2 and weak differences between the linear model and piece-wise model in Figure 3, which results in similar filtering powers.

Conclusions
Bias correction for Goldstein filtering power is crucial for interferometric data post-processing and the subsequent multitemporal analysis.Having an accurate filtering power can mitigate both underand over-estimation of phase noise in interferograms showing different features.This paper presents a method for bias removal from accurate coherence estimation and model development that links coherence and filtering power.Based on the weighted coherence estimator that combines a second kind statistic with an Anderson-Darling test statistic, the unbiased sample coherence with a better spatial distribution is deduced.Driven by sample coherence, the optimal filtering power for each filtered patch is then determined from a piece-wise model, which detects the best filtering performance under each coherence level using a Monte Carlo simulation.The authors finally evaluate the effectiveness of the bias corrected Goldstein filter using both synthetic and real data and compare it against other adaptive filters that fix filter power by biased coherence, unbiased coherence with regular statistics, and iteration.A series of experimental results validated the improvements resulting from the newly developed method for noise reduction and edge preservation.
passed by the Anderson-Darling test under a significance level of 0.01  =

Figure 1 .
Figure 1.Similarity measurement using the patch-based Anderson-Darling test: (a) The procedure of AD estimation; (b) some similar patches (green) with respect to the reference patch (red); (c) the corresponding AD values of the whole window.

Figure 1 .
Figure 1.Similarity measurement using the patch-based Anderson-Darling test: (a) The procedure of AD estimation; (b) some similar patches (green) with respect to the reference patch (red); (c) the corresponding AD values of the whole window.

Figure 3 .
Figure 3.The best filtering power  under each coherence level.The dashed line denotes the empirical model presented in [9].

Figure 3 .
Figure 3.The best filtering power α under each coherence level.The dashed line denotes the empirical model presented in [9].

15 Figure 6 .
Figure 6.The results estimated from different methods.Each column corresponds to a filtering procedure: first (Baran filter), second (AGFC), third (AGFP), and last (new filter).Each row corresponds to a parameter: first (coherence), second (filtering power), and third (filtered phase).

Figure 7 .
Figure 7.The boxplot of RMSEs between the truth and interferogram filtered by different methods.

Figure 7 .
Figure 7.The boxplot of RMSEs between the truth and interferogram filtered by different methods.

Figure 7 .
Figure 7.The boxplot of RMSEs between the truth and interferogram filtered by different methods.

Figure
Figure8apresents the original differential interferogram over the Jiuzhaigou area, where an earthquake with magnitude 7.0 occurred on 8 August 2017.After burst and sub-swath merging of the TOPS image, the non-linear deformation signals can be observed in the whole image.However, the low coherence in Figure8bhinders a clear interpretation of the fringes.The purpose of this trial was to test the capability of the new filter to correct under-estimation of the Goldstein filtering power over an incoherent area.

Figure
Figure8apresents the original differential interferogram over the Jiuzhaigou area, where an earthquake with magnitude 7.0 occurred on 8 August 2017.After burst and sub-swath merging of the TOPS image, the non-linear deformation signals can be observed in the whole image.However, the low coherence in Figure8bhinders a clear interpretation of the fringes.The purpose of this trial was to test the capability of the new filter to correct under-estimation of the Goldstein filtering power over an incoherent area.

Figure 12 .
Figure 12.The ROC curves for evaluating the filtering performances on edge preservation over the coherent area.

Figure 12 .
Figure 12.The ROC curves for evaluating the filtering performances on edge preservation over the area.

Table 1 .
SAR dataset and study area.

Table 1 .
SAR dataset and study area.

Table 2 .
Statistics of SPD and residue reduction in interferograms estimated from different filters using Sentinel-1A data over the Jiuzhaigou area.

Table 2 .
Statistics of SPD and residue reduction in interferograms estimated from different filters using Sentinel-1A data over the Jiuzhaigou area.