Enhanced Goldstein Filter for Interferometric Phase Denoising Using 2-D Variational Mode Decomposition

Denoising of interferograms is a vital step in the processing of synthetic aperture radar (InSAR) data. The primary goal is to filter the noise to the extent possible while retaining the fringes of the interferograms. Among the widely available classes of filters, the frequency-domain filters are still being used, owing to their robustness and generalizability to varying phase noise characteristics. This article deals with an enhancement to the well-known frequency-domain filter, i.e., the Goldstein filter, which is basically a phase filtering algorithm for interferometric products. The proposed extension to the Goldstein filter deals with deriving the tuning parameter based on the spatial frequency modes. This is achieved by using the mode-level characteristics rendered by the 2-D version of variational mode decomposition (2D-VMD) on the interferograms under test. The results of simulation and real interferogram data show that the proposed approach reduces the noise levels while minimizing the loss of signal.


I. INTRODUCTION
S YNTHETIC aperture radar interferometry (InSAR) is an effective means of obtaining digital elevation models (DEMs) and assessing ground deformation. Owing to its high spatial resolution and accuracy, it finds tremendous applications in monitoring of volcanoes [1], ground and glacier movements/thickness change [2], oil/water extraction, natural or human made subsidence due to excavations, mining [3], and so on [4]. All this is possible due to the changes seen in phase information in the synthetic aperture radar (SAR) data leading to interferograms depicting the changes in the target being viewed. Thus, the applications of SAR range from satellites at a macro-level to micro-level analysis as in [5]. The technique of interferometry uses the interference patterns created by two or more waves to assess the differences in their amplitude and phase. The interference pattern between these two images is analyzed for detecting the changes in the Earth's surface. However, the changes in the phase can be due to other unforeseen factors like propagation delay caused due to nonstationarity in the distribution of atmospheric temperature, pressure, water vapor, atmospheric turbulence, and decorrelation effects. This significantly contaminates the quality of deformation and topographic height in the interferograms. Also, the denoising of phase information is a primary prerequisite for the usage of InSAR in fields like deformation monitoring, topography mapping, and so on.
The noise correction in interferograms began with the advent of spatial filters. However, they suffered from the loss of spatial resolution and performed very poorly in regions of large slopes. The filtering schema shifted from the spatial domain to the transformed domain due to the need for robust filters. The major advantage of having transformed-domain filters is that the phase signatures are very well separated from the noise. Wavelet filters and frequency-domain filters were largely developed as an alternative to spatial filters. The former filters came with increased space and time complexity, and hence, the frequency-domain filters gained popularity due to the ease and simplicity of implementation. As a result, several frequency-domain-based approaches were developed, and the Goldstein filter was found to be the most effective in this regards [6].
Recently, deep-learning-based interferogram noise cleaning approaches are researched exclusively. Many popular convolutional neural networks (CNNs)-based approaches [7], [8], [9] and autoencoder-based approaches [10] are currently explored due to its powerful feature extraction capability. One similar approach using CNNs is DeepInSAR [11], which outperformed the popular Boxcar filter. However, such filters are still in the stage of infancy, due to the high amount of training data required to train the network and nongeneralizability to new unseen data of completely new distributions of noise. Frequency-domain filters have proven to be generalizable across different phase noise and they being unsupervised in nature, they are still being deployed in various interferogram products. Thereby, owing to their high efficiency and simplicity, the Goldstein filter is still a popular and widely used filtering algorithm for denoising of interferometric phase noise [12]. The main issue with the Goldstein filter lies in the fact that it uses a single tuning parameter across all the patches in the given interferogram, thereby applying similar stress in terms of filtering on all the regions of the interferogram. Thus, several variants of the Goldstein filters have been tried in order to improve the denoising capability. More details of these variants are discussed in Section II.
Owing to the wide applications of interferograms being hindered due to unforeseen noise in the data, the proposed work aims at building a robust filter by providing an extension to the Goldstein filter. This extension aims at overcoming the drawbacks of the existing filters by adaptively deriving the tuning parameter based on the noise profile of interferograms (as determined by the decomposition algorithms) under test. The main contributions of this article are as follows.
1) Proposing a novel and effective enhancement to the existing Goldstein filter for filtering the interferometric phase signals using 2-D variational mode decomposition (2D-VMD) algorithm. 2) Comparison of the proposed approach with other closely related advanced filters. 3) Testing the proposed approach on simulated and real world datasets. The remainder of this article is organized as follows. The enhancements proposed for the Goldstein filter, along with its mathematical background, are explained in Section II. The experimental results and comparison with other related works are discussed in Section III. This article concludes in Section IV along with the future roadmap of the study.

II. PROPOSED ENHANCEMENT TO GOLDSTEIN FILTER
This section briefs the steps involved in the proposed enhanced Goldstein filter. We initially give a quick background of the Goldstein filter. Since we use mode decomposition exclusively in our study, we provide a detailed mathematical review of the mode decomposition algorithm used in the study. This is followed by the pipeline designed in the study.

A. GOLDSTEIN FILTER ANALYSIS
Goldstein and Werner [6] proposed the well-known adaptive radar interferogram filter. The idea was based on multiplying the Fourier spectrum (r, s) of a given small patch of interferogram with its smoothed value S{| (r, s)|} being raised to an exponent α as given by where S{·} is the smoothing operator, r and s are the spatial frequencies, and H(r, s) is the spectrum of the filtered interferogram (filter response). Small portions of the interferogram under test are taken as patches and are overlapped to avoid the discontinuities at the boundaries. The tuning parameter α plays a major role in determining the performance characteristics of the filtering algorithm [13]. Its value is chosen empirically in the range of [0, 1]. For α = 0, the multiplication factor being unity has no filtering done on the original spectrum (r, s). Larger values of α lead to significant filtering. The selection of optimal tuning parameter value of α is the major issue while dealing with Goldstein filter as high values of α leads to the loss of resolution in the filtered phase [6]. On the contrary, smaller values of α lead to the phases remaining still noisy [13]. Hence, determining the value of α is a major issue as the noise characteristics vary spatially across the interferogram under test.
There have been several attempts in exploring the ideal value of α for Goldstein filters. Baran et al. [14] obtained the tuning parameter by harnessing the relation between the phase standard deviation (PSD) and the coherence using the probability density function (PDF). This extension to Goldstein filter ensured that the interferograms with high coherence were benefited with the filtering process as it minimized the loss of resolution. However, in real-life scenarios, the interferograms come from nonstationary and nonlinear target scenes. This distorts the PDF due to the inhomogeneity in the interferograms which is a result of inaccurate estimation of coherence [15]. A recent approach in obtaining the tuning parameter α is by Song et al. [15] in which the authors use the bidimensional empirical mode decomposition (BEMD). The major shortcoming of this approach is the inability of BEMD in distinguishing the noise components from the actual information inherent in the data under test. This is shown through simulation studies in the next section. The approach in [13] deals with obtaining the tuning parameter α using the coherence information and the residual noisy phase content. The authors mainly use the dominant frequencies in the residual noisy phase. When the noise spectrum spans a larger bin, the dominant frequency cannot solely represent the entire spectrum. Hence, there is a need to aggregate all the information of the noise spectrum effectively in determining the tuning parameter. The proposed work aims at achieving this goal by using the entire noise profile (obtained through the advent of mode decomposition) and analyzing its signal-to-noise ratio (SNR) as discussed in the subsequent sections.

B. 2-D VARIATIONAL MODE DECOMPOSITION
2D-VMD [16] is an extension of VMD [17]. The latter basically decomposes a 1-D signal into a set of K number of subcomponents u K (called modes). Each mode u k is compact around a center frequency ω k and the sum of all the modes (k = 1, 2, . . . , K) equals the original signal. 2D-VMD is the 2-D extension of VMD with a goal of decomposing 2-D signal x (like images) into K modes, u K [16]. The objective function is similar to 1D-VMD, i.e., to have all the modes with a compact center frequency with the bandwidth as minimum as possible; and the sum of all the modes equals the original 2-D signal under test. The 2D-VMD is variational, fully adaptive, and nonrecursive in nature that decomposes the signals in a well-defined mathematical manner. It searches for the optimal center frequency and bandwidth across modes in an iterative fashion to have effective separation in the frequency domain of the signal under test. The working principle of 2D-VMD is provided in the following steps. For more details on the algorithm, the reader is requested to refer the work in [16]. Every mode u k is supposed to be compact around the center frequency ω. The center frequency for a given mode u k is denoted by ω k and the bandwidth of every mode centered around ω k is computed as follows [16].
1) 2-D Hilbert transform is used to get the 2-D analytic signal of the mode u k (x). The 2-D analytic signal of the kth mode u k is represented by u AS,k as where the operator * signifies the convolution operation, δ(·) is the unit impulse function, j = √ −1 is the imaginary unit, π = 22/7, and ·, · is the inner product operation.
2) The frequency spectrum corresponding to the mode u k (x) is shifted to its baseband portions. This is achieved by mixing it with an exponential tuned around the computed center frequency as is the gradient. 4) The K number of modes is computed by solving the following constrained variation problem: where f (x) is the original signal (image), and β k is the penalty factor.
The constrained problem in (3) is solved by introducing the penalty factor β and Lagrangian multiplier λ. This makes the problem unconstrained in nature and also ensures the correctness of signal reconstruction even in the presence of noise. This ensures that each individual mode is around a compact frequency and thereby it avoids mode-mixing problem [16]. The Lagrangian is thus given as The alternating directive method of multipliers, an iterative suboptimization technique, is used to solve the above minimization problem and hence, 2D-VMD is nonrecursive in nature. Iterative updation of u n+1 k , ω n+1 k , and λ n+1 k is done to get the optimal solutions. The expression of u n+1 k is given by On similar grounds, the expression for ω n+1 k is given as where ω ∈ k , k = {ω|ω · ω k ≥ 0}. Thus, the overall working of 2D-VMD can be summarized as in [16]. 1) Step 1: Initialization of u 1 k , ω 1 k , and n, where n is the number of iterations.

2)
Step 2: Updation of u k and ω k as per (5) and (6), respectively, in the frequency domain.

3)
Step 3: Updation of the Lagrangian multiplier λ(x) is done using a standard gradient ascent with fixed time step τ as

4)
Step 4: Stopping of the iteration when the stopping criteria is met, else return to step 2. Since, 1D-VMD lacks the ability to exploit the spatial information of interferograms, we use the 2D-VMD algorithm in this study. Because of its robustness in decomposing a 2-D signal into its constituents, the 2D-VMD algorithm is used exclusively in fields like multi/hyperspectral classification [18], denoising [19], fusion of multispectral and panchromatic images [20], radar-based human motion recognition [21], and so on. The most closely related state-of-the-art technique related to this work is the one proposed in [15] which is termed as Goldstein-EMD hereafter. Since, this method's core component is the usage of BEMD, we throw some light on why this article is focused on the usage of 2D-VMD. Fig. 1 shows the decomposition of a sample image rendered by BEMD and 2D-VMD algorithms. The implementation of these algorithms is done in MATLAB. The input image has certain dominant entities like oval, rectangle, and jagged ellipses. Consider mode 1 from BEMD and VMD. The former has a mixture of multiple entities like the oval, jagged ellipses, and some portion of rectangle being distorted. However, in case of 2D-VMD, we see the two entities, i.e., the oval and the rectangle appear distinct and clear. In mode 2, we see that one jagged ellipse is distinct and prominent in 2D-VMD, whereas in case of BEMD, distorted versions of multiple entities appear. The same is the case with modes 3, 4, and 5 wherein the 2D-VMD truthfully recovers the entities either individually or prominently in the given mode, when compared to that of its BEMD counterparts. Hence, it can be observed that 2D-VMD is better in terms of retrieving the modes which retain the appearance of the entities in a more clear and precise manner. However, BEMD generates modes which distort the entities significantly. Also, the issues due to mode mixing, i.e., 2 or more modes having same components is clearly visible in the figure related to BEMD. This is usually evident in case of intermittency signals and noises and refers to the conditions wherein the modes comprise of different frequencies instead of compact central frequencies in any given mode. Also, the lack of mathematical background in BEMD and other shortcomings like the end effects, high dependence on interpolation, finding of extrema, and stopping criterion makes BEMD inefficient in terms of decomposition of modes [22]. Moreover, the strong theoretical foundations of 2D-VMD make it an enhanced version of BEMD. It is to be noted that 1D-VMD is not designed to explore the spatial information of 2-D data like that of interferograms. Hence, 2D-VMD is used exclusively in this study, as discussed in the subsequent sections. Fig. 2 illustrates the steps in the proposed method. Note that the input image used (with triangle, oval, rectangle, and noise as embedded entities) and its decomposed mode shown in the figure are for illustration purposes only. For a given patch I p of the original image I, we decompose it into K modes using the 2D-VMD algorithm, which yields u K number of modes, which signifies

C. PROPOSED METHOD
where the first term corresponds to clean signal modes (u S ), and the second term corresponds to the noisy modes (u N ).
After the decomposition of the patch I p into u K modes, the modes corresponding to noise (u N ) contain valuable information of the noise characteristics while the phase information is present in the u S modes. This is illustrated as four modes in the example shown in Fig. 2, wherein each mode has either of the entities (triangle in mode 1, oval in mode 2, and rectangle in mode 3) and the noise entity is obtained as the last mode. For a given image patch with m × n dimension, letĪ p be the mean pixel intensity for the image patch I p , andū k be the mean pixel intensity for the kth mode. The noisy mode (u N ) is identified from the K number of available modes using the following minimization problem as: (10) which signifies that the mode k ∈ K with least correlation coefficient with that of the corresponding patch I p is the one to be considered as noisy mode (u N ). There can be scenarios wherein more than one mode contributes as noise. In such cases, the modes are summed up (called summed noisy mode u N ) and considered for computation of the variance. An important measure to quantify the noise in a given patch is computing the local SNR. Let σ 2 φ be the variance of phase in a given image patch. In order to derive an SNR which signifies how much the given patch I p is nonstationary, we compute the ratio of variance of I p along with the variance of the noisy mode (u N ). Thus, the SNR ξ of the pth patch is given by Larger the nonstationarity in a given patch I p , higher the values of ξ p . As these values span arbitrarily, we normalize them in the range of [0, 1] as This ensures that patches with higher SNR are filtered less when compared to the ones with lower SNR values. Alternatively, the normalization process can be carried out as to account for any possible underfitting in the filtering process. We then compute the tuning parameter α as This modification to the tuning parameter ensures the avoidance of under filtering when the values are very less (say < 0.2). The smoothing parameter S{·} is now raised with this tuning parameter in the Goldstein filter in (1) per patch as where H p (r, s) is the spectrum of interferogram patch I p .

III. RESULTS AND DISCUSSION
The efficiency of the proposed modification to Goldstein filter is tested through simulations and is analyzed quantitatively using the metrics of evaluation, discussed in the next section. Qualitative analysis is shown using visualizations of the filtering outputs. The real-world data is also subjected to filtering, and appropriate evaluation metrics are used to assess the performance of filtering algorithms.

A. EVALUATION METRICS
We use metrics of evaluation like the PSD, edge preservation index (EPI), and the correlation of the filtered interferogram with that of the ground-truth interferogram. PSD is indicative of the extent of smoothness of phase and smaller values indicate less noise and thereby smoother phase values. It is defined as where the linear phase ramp is represented byφ(m, n) in a given moving window. The EPI signifies the capability of edge and fringe preservation. For a given filtering algorithm, the values of EPI closer to 1 indicate better preservation of edges and fringes in the filtered interferograms. It is calculated as per [15] as where φ r (m, n) is the reference phase from the ground truth, and φ f (m, n) is the phase of the filtered interferogram. Another metric named residue signifies the inconsistencies or the jumps in phase values and smaller residue values are indicative of less noise. This number is a vital component during the unwrapping of phase and hence, is a primary measure of evaluating the quality of an interferogram. All these metrics can be used to evaluate the simulated data due to the presence of the ground truth. However, in case of real world interferograms, only the residue and PSD can be used for evaluation of the filtering capabilities of the algorithms under test.

B. SIMULATION RESULTS
We first evaluate the performance of the proposed pipeline on simulated dataset. The interferograms are generated as per [23] and [24] which is based on generating wrapped phase interferograms considering user-defined probabilities of built up areas, slopes, atmospheric turbulence (i.e., fractal  Perlin noises), distorted 2-D Gaussian surfaces, deformation by earthquakes, and completely decorrelated areas like water. We compare our proposed approach with that of Goldstein filter (with α = 0.5); Goldstein filter based on EMD (Goldstein-EMD) [15] and using the improved Goldstein filter from [13] (Improved-Goldstein). Song et al. [15] showed the robustness of Goldstein-EMD filter over the Baran filter [14] and hence, the latter is excluded in the study. Since, this study focuses largely on unsupervised methods, we refrain from exploring any deep learning methods in interferogram denoising. Fig. 3(a) shows the unfiltered interferogram phase which is subjected to filtering and the resulting cleaned images are shown in Fig. 3(b)-(d) for Goldstein, Goldstein-EMD, Improved-Goldstein, and proposed approach, respectively.
The residues, PSD, EPI, and the correlation coefficient values are given in Table 1. It is to be noted that the PSD values of the proposed filter are lesser than the existing methods used for comparison. Also, the PSD value of the proposed filter is very close to that of the ground-truth interferogram. The EPI value of the proposed method is closer to 1 when compared to the Goldstein and the Goldstein-EMD filters, which indicates that it is better in retaining the edge and fringes information. Although the PSD and EPI characteristics of the Improved-Goldstein filter are better, however, the overall correlation coefficient is lesser and this results in the distorted  filtered output as seen in Fig. 3(d). Thus, the proposed filter not only has better smoothing performance but is robust in terms of retaining the edge and fringe information in a better way.
In order to visually inspect the filtering abilities of the interferograms by different algorithms, we provide the cross section view of the interferograms. A sample cross section view from row 2 is shown in Fig. 4(a) which corresponds to the simulated interferogram's wrapped phase values (in radians). It can be observed that the filtered phase by the proposed method matches closely to that of the groundtruth phase values in comparison to the rest of the methods. Fig. 4(b) shows the offset between the ground truth and the filtered phase by different approaches over the region marked in circle in Fig. 4(a). This portion mainly corresponds to the regions where phase distortion commonly occurs. Thus, the proposed technique is capable of preserving the phase information in the interferograms.

C. EVALUATION OF REAL DATASET
We evaluated the performance of the proposed filter and the existing ones with a real interferogram dataset obtained using the LICSAR tool [25]-based data from Sentinel-I. The interferograms under study consist of data collected from the ascending passes of the satellite. Four earthquake eventsrelated interferograms from Las Brisas (Mexico), Namie   (Japan), Marotta (Italy), and Tonopah (Nevada) are considered which contain the phase change information pertaining to the deformations. The details of the study locations are provided in Table 2.
For the sake of brevity, we only present the visualization of one particular study location. Fig. 5 shows the interferogram from Location -2, i.e., Namie (Japan) and the results of filtering by different approaches. A particular region of interest (marked as patch "A") in the figure has been zoomed upon and shown in Fig. 6 to provide the details of filtering capabilities of different approaches. The quantitative results for each of the study locations used are presented through Tables 3 and 4. It is to be noted that the proposed method is efficient in reducing the residues from the unfiltered images when compared to other robust filters. Also, the reduction in PSD is higher for the proposed method. The Improved-Goldstein filter [13] performed better in the reduction of residues and PSD, albeit, it is observed that certain unwanted edges are created in the filtered image as is evident in the zoomed version of the images (Fig. 6). Hence, in case of filtering of real interferograms which comes with lack of ground truth, only the residue and PSD information are incapable of determining the overall efficiency of the filtering algorithm under test. A deeper visualization of each of the algorithms under test has to be carried out in selecting the best possible approach for the target application.
Though the classical Goldstein filter and its variants are effective in terms of denoising the phase noise, their performance is often curtailed due to the inability to retain the edge information near the fringe-dense regions. This largely depends on the inaccuracy in obtaining proper values of the tuning parameter α. In interferograms, the values of backscatter are uniform in a region of local homogeneity and can be treated as stationary. The phase values extracted from such stationary locations or patches have similar morphological characteristics like edges, textures, and noise. The proposed technique harnesses this information by segregating the common elements as separate modes from the noisy ones. The α values thus derived from such locations are common across similar patches and thus filtering is done in accordance to the intensity or the severity of the noisy component only.
The simulation and real-world data results show that the proposed technique avoids the inaccurate estimation of α owing to contaminated PDFs in the interferograms under test, which is basically due to decorrelation, atmospheric screening, and inhomogeneity issues. The proposed modification is thus proved to perform efficiently when compared to the classical Goldstein filter and the state-of-the-art modified Goldstein filter.

IV. CONCLUSION AND FUTURE SCOPE
This study is centered around designing an enhanced version of the widely used Goldstein filter. The major contribution of the study lies in effectively computing the tuning parameter α which quantitatively assesses the roughness or the texture of the patch of the interferogram. Through the advent of 2D-VMD, the useful information and the noisy components are segregated into disjoint orthogonal modes. This ensures that the noisy modes are separate and thus the roughness is computed based on the noise content and the value of α is derived relatively with this information. Simulation studies and the testing done on real-world interferograms confirm the robustness of the proposed technique over closely related robust filters. As a future roadmap, we intend to extend the inclusion of the proposed pipeline in the phase unwrapping case of interferogram analysis, in order to do away with errors in the unwrapping phase if any.