Study of Systematic Bias in Measuring Surface Deformation With SAR Interferometry

This article investigates the presence of a new interferometric signal in multilooked synthetic aperture radar (SAR) interferograms that cannot be attributed to the atmospheric or Earth-surface topography changes. The observed signal is short-lived and decays with the temporal baseline; however, it is distinct from the stochastic noise attributed to temporal decorrelation. The presence of such a fading signal introduces a systematic phase component, particularly in short temporal baseline interferograms. If unattended, it biases the estimation of Earth surface deformation from SAR time series. Here, the contribution of the mentioned phase component is quantitatively assessed. The biasing impact on the deformation-signal retrieval is further evaluated. A quality measure is introduced to allow the prediction of the associated error with the fading signals. Moreover, a practical solution for the mitigation of this physical signal is discussed; special attention is paid to the efficient processing of Big Data from modern SAR missions such as Sentinel-1 and NISAR. Adopting the proposed solution, the deformation bias is shown to decrease significantly. Based on these analyses, we put forward our recommendations for efficient and accurate deformation-signal retrieval from large stacks of multilooked interferograms.


I. INTRODUCTION
A S AN established geodetic technique for Earth surface deformation monitoring, the accuracy of the interferometric synthetic aperture radar (InSAR) time-series analysis should be well quantified and the potential error sources must be known. Any uncertainty in the accuracy of the InSAR products compromises their reliability in sensitive applications.
Persistent scatterer interferometry (PSI) is among the pioneering techniques for improving the accuracy of the InSAR [1] in deformation retrieval. Exploiting the phase stable persistent scatterers (PS) within the time series, PSI avoids a major limitation of InSAR, namely, the signal decorrelation [2]. Using the high signal to noise ratio (SNR) PS measurements, another error source of InSAR is mitigated through the separation of the atmospheric signals from the deformation. The PSI technique has been perfected since its invention, and its accuracy has been studied thoroughly [3].
The low density of PS in nonurban areas motivated the invention of complementary techniques to PSI. These methods exploit partially decorrelating areas in time-series analysis. Referred to as distributed scatterers (DS), such areas pertain to an ensemble of natural scatterers that share similar scattering characteristics. A variety of methods have been put forward to allow the use of DS in deformation estimation, with small BAseline subset algorithm (SBAS) [4] and SqueeSAR [5] as the overarching approaches. The shortcoming of natural scatterers of DS is their inherent phase noise caused by signal decorrelation. Common to all DS techniques, spatial averaging, or multilooking, is employed to reduce this stochastic noise in the interferograms.
The purpose of this article is to investigate the accuracy of the multilooked interferograms with regard to the DS techniques. Different studies have been dedicated to the validation of DS with independent geodetic techniques, such as Global Navigation Satellite Systems (GNSS) measurements (see [6]). Being spatially sparse, such independent measurements restrict the comprehensive study of the DS behavior. Here, we consider a different validation approach to reveal a peculiar systematic signal in the multilooked interferograms. If unattended, the mentioned signal can severely bias the deformation estimates of DS. We investigate the accuracy of the deformation estimates in the presence of this signal to highlight the role of different DS techniques in either exacerbation or mitigation of the deformation bias.
Following the theoretical background of Section II, we design different comparison approaches in Section III to shed light on the following propositions.
1) The multilooked interferograms reveal a systematic signal that cannot be explained by the topographic or atmospheric variations and interfere in the accurate estimation of the deformation. Such signals are short-lived and decay with the temporal baseline, rendering the short temporal baseline interferograms to be more error-prone. Hereafter, we refer to this effect as the fading signal to reflect that it is inherently a short-lived but physical phase contribution. phenomena. The understanding of the source of this signal and its modeling is a current research topic [7]- [9]. Lacking a generic model to explain the behavior of the fading signal makes the calibration of these errors intricate. 5) Using the temporal data redundancy within the time series, i.e., by exploiting all possible interferograms, the fading signal is significantly mitigated in the multilooked interferograms. This mitigation improves the accuracy of both displacement time series and displacement velocity estimates. After the careful examination of these propositions in Sections IV-VI, we put forward our suggestions for achieving accurate deformation monitoring with DS in Section VII.
The focus of this article is on observing the fading signal and its mitigation using the data-adaptive approaches. Occasionally, different hypotheses are put forward to explain the physical phenomena behind the fading signals; however, the in-depth validation of the theses is beyond the scope of this article.

A. Single-Look Versus Multilooked Observations
In a time series of n SAR acquisitions, each pair of a so-called master and slave acquisitions allows the formation of an interferogram. The latter is the phase difference between the two complex-valued images, relative to the master acquisition. For each choice of master acquisition, n − 1 commonmaster interferograms exist within the time series. In total, n(n − 1)/2 multimaster interferograms may be formed for each time series. As the target of geodetic applications, Earth surface deformation is estimated from all or a subset of these interferograms. Therefore, the accuracy of the deformation estimates is governed by the quality and number of the exploited interferograms.
In this work, we distinguish between two types of observations in the interferograms: single-look versus the multilooked interferometric phases. An example of the former is PS, where the single complex-valued pixels are exploited within the time series. The latter observations are the result of spatial averaging common for the DS regions.
For the case of PS, the multimaster interferograms of the single-look observations are highly redundant. A mere subset of n −1 common-master interferograms allows the reconstruction of all possible interferogram combinations within the time series, for example Here, W (x) = mod{x + π, 2π} − π, m indexes the master acquisition; here and after, the amplitude of the interferograms is disregarded for brevity of notation. φ imk is an indicator for the consistency of phase components within the time series [8].
As the second category of observations, DS are characterized by the homogeneous areas that undergo signal decorrelation. Therefore, the single-look observations within the DS region have low SNR. The remedy in improving the SNR is to perform spatial averaging within the homogeneous DS region and form multilooked observations. All advanced DS-InSAR techniques employ multilooking, although some techniques include single-look observations [10], [11] or use data-adaptive and, thereby, feature-preserving methods [5], [12]. Although effective in noise reduction, multilooking changes the statistical properties of the interferograms. It reduces the redundancy among the n(n − 1)/2 interferogram combinations, such that (1) no longer holds. A residual component is observed among each of the three arbitrary multilooked interferograms, that is implying that the phase information of the DS regions is inconsistent among the arbitrary interferograms. Two effects may be observed for the DS phases. 1) E{ φ imk } = 0, reflecting the stochastic noise caused by signal decorrelation. 2) E{ φ imk } = 0, indicating a variant systematic signal among the multilooked interferograms. Here, E shows the statistical expectation of the accompanied random variable.
The former stochastic effect translates to noise in the deformation estimates. More critical is the latter systematic effect that may be present in a subset of the interferograms, and therefore revealed in certain interferogram triplets. The systematic effects introduce a fading signal in the affected multilooked interferograms. If present and not mitigated properly, they are interpreted as deformation and bias the estimates.
The peculiar fading signal is raised by multilooking and absent in single-look observations. Therefore, the discrepancy between the single-look and multilooked observations can be used to investigate first the presence of the fading signals and finally their impact on biasing the deformation estimates. This rationale is adopted throughout this article and expanded in Sections III-B and III-C.
As the target for geodetic application, the atmospheric and surface-deformation signals are consistent within the arbitrary interferograms [8], i.e., φ imk = 0 holds for the phase components corresponding to these signals. The reason is that these physical sources uniformly affect all scatterers within the DS region; therefore, their net residual effect vanishes in E{ φ imk } [8], [13]. To safeguard against the fading signals, one needs to retrieve the consistent component of the interferograms or, in other words, reconstruct the consistency among them.

B. Reconstruction of Consistency
One approach to reconstruct the consistency is to discover the physical source of the fading signals, accurately model their scattering behavior and, thereby, phase contribution, and compensate the corresponding phase from the multilooked interferograms. As one source of the fading signals, we refer the reader to the variation in the moisture content of the scattering media. The effect of moisture on InSAR is studied and modeled in [7]. Based on the proposed moisture model, the corresponding phase contribution is estimated from the residual phase of (2) in [9]. The multilooked interferograms may further be compensated by the estimated moisture-induced phase to form consistent interferograms [9]. The consistency in this case is reconstructed by calibrating the fading signal's phase component.
As our investigations on two distinct regions reveal, the behavior of the fading signal varies significantly across different data sets, land cover, etc. Therefore, the mentioned moisture model does not suffice in explaining the observed inconsistencies (see Section VI), rendering the calibration ineffective. For instance, to explain the observed fading signal over part of our test sites, we propose an alternative model to moisture variation in Section VI. This analysis serves as an example to show that for a an effective calibration, we either need to perform comprehensive research on various possible physical sources of inconsistencies and design case-specific models to explain the interferometric phase of each source or to design a generic model applicable to all possible sources of inconsistencies. In both cases, further studies on various test sites are inevitable. This model-based approach in the reconstruction of consistency is subject to comprehensive research, hence, as of present, still inapt for systematic deformation monitoring.
A practical approach in the reconstruction of consistency is to retrieve the consistent components of the phase within the time series, i.e., by explicitly imposing the condition of φ imk = 0 among all the n(n − 1)/2 interferograms within the time series. This approach is initially designed to reduce the stochastic noise within the interferogram stack [5], [14]. However, as it will be revealed by the analysis of this article, it significantly reduces the effect of systematic phase inconsistencies as well. Referred to as phase linking, the method is explained in Section II-C. Phase linking does not consider a specific model for the fading signal, and it rather gains robustness to such signals by exploring the temporal data redundancy of the interferograms. Successful phase reconstruction by phase linking may be due to the exploitation of sufficiently large number of interferograms within long time series; especially, the long temporal baseline interferograms that are less prone to the short-lived fading signal. However, the examination of this thesis is beyond the scope of this article. Within the scope of this article, we substantiate that phase linking is an efficient and practical solution to reconstruct the consistency.

C. Phase-Linking Techniques
After its pioneering authors, we define phase linking as the estimator that retrieves n − 1 independent common-master interferograms from the partially redundant n(n − 1)/2 multimaster interferograms [14]. Here, we shortly introduce phase linking.
Without loss of generality, we concentrate on one DS region (see [5], [12] for the algorithmic details of the DS region selection). The DS is comprised of an ensemble of the spatially homogeneous region of p pixels in a time series of n synthetic aperture radar (SAR) images, arranged in a matrix Z ∈ C n× p . Based on the central limit theorem, Z follows the zero-mean n-variate Complex Circular Gaussian (CCG) distribution [2]. Under the validity of this distribution, the sample covariance matrix or its normalized version sample correlation matrix (SCM) suffices for the full description of the DS. The SCM, denoted by C, is a Hermitian matrix whose off-diagonal elements pertain to all possible multilooked interferograms I ik within the time series and their corresponding coherence ik , that is Various phase-linking approaches are defined based on different modelings of the SCM [15], [16]. In an earlier work, we proposed a computationally efficient approach to phase linking called Eigen-decomposition-based Maximumlikelihood-estimator of Interferometric phase (EMI) [15]. This proposal decreases the computational cost by reformulating phase estimation into the following Eigen-decomposition problem [15]:φ where v i is an arbitrary complex vector of size n × 1 and • is the Hadamard product. φ is a vector of n wrapped phase values. It contains the consistent interferometric phase component within the exploited interferograms. This phase information is used for the retrieval of the deformation signal. EMI allows the convenient estimation of this phase series by taking the smallest Eigenvector of the matrix C • −1 .
Note that, as common in InSAR, the estimation of absolute phase is ambiguous. An arbitrary image in the time series is selected as a reference scene, and its phase is set to zero and the remaining phases are measured relative to these arbitrary data.

III. METHODOLOGY IN ACCURACY ASSESSMENT
In this section, we explain our experiments to reveal the existence and impact of the phase errors in the multilooked interferograms. The design of these experiments is based on the following rationale.
1) Any systematic discrepancy between the single-look and multilooked observations is an indicator of the fading signals (see Section II-A). 2) Should inconsistent systematic effects exist within the multilooked interferograms, the deformation estimates vary depending on which interferograms are exploited. As the benchmark of the experiments, we perform PSI to retrieve the deformation based on single-look observations. Exploiting various combinations of multilooked interferograms, we further perform multiple processing rounds to investigate the impact on the deformation estimates. In the first processing round, we exploit all possible interferograms within the time series, as explained in Section II-C. In the following processing rounds, we test the ingestion of different subsets of the interferograms; the corresponding estimator for this processing is introduced in Section III-A.
As explained in Section III-B, the resulting deformation estimates from the different processing schemes are evaluated against the benchmark PSI to study the estimation bias. In Section III-C, we expand on our experiments to reveal the existence of the fading signals in multilooked interferograms and quantify their magnitude and impact on the deformation estimates.

A. Enhanced Short Temporal Baseline Subset Algorithm
To allow the ingestion of different subsets of the interferograms and evaluate the impact on the deformation estimates, we designed a variation of the SBAS technique of [4]. The difference with respect to the conventional SBAS is twofold.
1) The baseline constraint is only imposed on the temporal separation between the acquisition pairs. 2) Phase linking is performed on the chosen interferogram subset. Deformation estimation follows based on the wrapped phase estimates. Refraining from unwrapping the interferograms and performing phase linking on the subset instead, the designed approach is less prone to the propagation of the phase unwrapping errors. To reflect these differences with the conventional method, we refer to this approach as Enhanced Short temporal BAseline Subset algorithm (E-StBAS). Note that E-StBAS is composed of a single subset of a fully connected network of interferograms. Contrary to the generic SBAS approach, here, there are no multiple, possibly isolated, interferogram subsets.
The chosen interferogram subset is comprised of bw number of shortest temporal baseline interferograms per acquisition in the time series, such that the total number of exploited interferograms reads as Compared with the conventional phase linking, here, a band matrix will replace the full SCM (see Fig. 1). The bandwidth of the matrix is defined by the parameter bw. The consistent phase based on these interferograms is reconstructed by the following iterative optimization: The iterations can be initialized by the largest Eigenvector of the band matrix C bw , that iŝ In practice, (8) provides a sufficient approximation such that the iteration by (7) is unnecessary. As in Section II-A,φ contains the consistent interferometric phase component within the exploited subset. This phase information is used for the retrieval of the deformation signal.

B. Evaluation of Deformation Bias
Retrieving the consistent phase series for DS using either EMI or E-StBAS, the standard PSI processing [1], [17] is employed on the high SNR DS to mitigate initially the atmosphere and estimate eventually the deformation [5]. Two products may be retrieved from the deformation signal, namely, the relative displacement time series of size n − 1 per DS and, more concisely, the modeled mean displacement velocity as a single parameter per DS.
The intention is to evaluate the accuracy of both DS-derived products. We opt for the PS deformation measurements as the benchmark for validation.
The performance evaluation is conducted as follows: a spatial grid of size 1 km 2 is chosen for downsampling both displacement time series and displacement velocity maps. Atmospheric and ground motion signals are assumed to be stationary within this spatial window. For all PS and DS within the defined reference grid cells, a weighted average of the deformation signal substitutes the sparse estimates. The weighting is based on the a posteriori coherence of PS [17] and DS [5]. For the latter scatterers, the clutter is disregarded using a constant false-alarm-rate detector. Following this approach, the downsampled DS and PS deformations are directly comparable. From this point on, the calculation of the estimation bias in the deformation products is straightforward where . shows the mentioned weighted averaging operator, x, y are the spatial coordinates of the downsampled grid, and d is the evaluated bias. Subscript d can represent either the displacement values in the time series or the displacement velocity. This evaluation will result in a time series of displacement bias as well as the overall displacement velocity bias for the entire time series. Both biases are calculated over the downsampled spatial grid.

C. Evaluation of Interferometric Phase Bias
We intend to track the bias in displacement velocity to the exploited multilooked interferograms. We first introduce a measure to quantify the interferometric phase bias that pertains to the fading signals. Using error propagation, we introduce a second measure to calculate the expected displacement velocity bias from the phase biases.
In the introduction of our first measure, we assume to know the interferometric phase φ , which is free of the systematic and stochastic inconsistencies (see Section II-A). Knowing this phase, the error of the multilooked interferograms can be evaluated for each DS at each interferogram, that is Here, φ shows the multilooked phases over the DS regions, subscripts i, k refer to the master and slave acquisition index, and x, y are the spatial coordinates of the DS. The phase φ can be substituted, for example, by high SNR single-look observations (see Section IV-D for the practical approach to the evaluation of the phase error).
In principle, one can estimate the phase bias from the phase errors by allowing an averaging operator on a chosen ensemble; the averaging is necessary to reduce the stochastic noise. In our bias estimation, we choose a temporal averaging within the time series. The averaging is performed on the calculated phase error of each DS within the interferograms with an identical time lag l separation from their respective master scene. A normalization with respect to the temporal baseline of the interferograms is further considered. Following the temporal averaging and baseline normalization, the intended measure reads as where δ l is the average phase bias of each DS region per time lag l. The temporal normalization allows the direct comparison of the calculated phase biases without the concern for their variable temporal baseline. It as well eases the error propagation from the interferogram level to the displacement velocity, as explained in the following.
We quantify further the expected bias in the displacement velocity given the biases of interferograms by introducing a second measure. To commence, let us simplify the phase reconstruction of (7) by assuming uniform weighting, i.e., ∀i, k ∈ {1, . . . , n} : ik = 1. Furthermore, we assume that the phase errors are small enough that the small-angle approximation holds, i.e., exp j φ ≈ φ . Furthermore, the mentioned temporal normalization of (11) is necessary to propagate the phase bias to the displacement velocity bias. Following the mean propagation law for the linear models [18], the displacement velocity bias reads as where λ is the radar wavelength. To recapitulate, the introduced error propagation assumes the following.
1) The present phase error of the multilooked observations is small. 2) All interferograms are equally weighted in phase estimation. The first simplification holds for short temporal baselines, where the decorrelation noise is negligible enough not to violate the small-angle approximation. The second simplification is more likely to have an effect on the estimation of the phase and displacement velocity bias. Nevertheless, d v provides a mean to verify the evaluated biases in Section III-B with the predicted displacement velocity bias from the phase errors in the multilooked interferograms. In the case of positive verification, we may conclude that the cause of deformation biases stems from the multilooked interferograms and rules out other plausible sources in raising the deformation biases.

IV. ANALYSIS OF BIASES
For the study of the bias, a test site is chosen in the island of Sicily, Italy. The test site is regularly monitored by Sentinel-1 A and B, providing abundant data with a minimum revisit of six days. The data set is composed of 184 acquisitions from October 2014 to September 2018 of ascending.
The geographical location, main cities, and the land-coverclassification map of the region are provided in Fig. 2. The prevailing land cover of the test site is in the form of crop and grassland. The studied area covers approximately 15 000 km 2 . Note that in different analyses of this section, different regions of the test site are covered. These regions are marked in Fig. 2 (top) by the two polygons. We refer to these regions at the caption of each figure to follow. Sicily is monitored by different studies, such as [19], [20], which potentially allows independent performance comparison.
Following the introduced methods in Section III, the biases in the displacement time series and velocity as well as the multilooked interferograms are studied in this section.

A. Comparison Scenarios
The intention is to define a scientifically credible experiment that can isolate the impact of multilooked phases on deformation estimation.
We perform three different DS analyses for deformation retrieval. All steps of the processing and the corresponding latent parameters are kept identical, and the only difference is in consistency reconstruction (see Section II-B). With reference to Sections II-C and III-A, different number of interferograms and different phase-estimation methods are used to retrieve the consistent interferograms, namely, the following: The defined experiment is summarized in Table I.
Identical to the three cases, interferometric wide area processing (IWAP) chain [21] is used for deformation estimation, and the estimation of consistent interferograms is integrated in this chain. As the first step toward DS processing, the statistically homogeneous ensembles surrounding each pixel are detected. The amplitude-based Anderson-Darling  [22]. The latter regions pertain to fast decorrelating scatterers such as water bodies and dense vegetation.
In addition to the abovementioned three DS comparison cases, we perform a conventional PSI [17] and treat the result as the benchmark for our analysis to follow.
Note that in the overall four processing rounds, the reference scene and the reference point are identical. Moreover, the latent parameters are kept identical or chosen in a data-driven fashion to ensure the credibility of the comparisons. Fig. 2 shows the retrieved displacement velocity map of these four described schemes. In the following sections, the results are quantitatively analyzed.

B. Bias in Displacement Velocity
We are interested in the quantitative error of the displacement velocity maps reported in Fig. 3. Following the described method in Section III-B, the PS scheme is taken as the benchmark. According to (9), the discrepancy between the velocity estimated by each of the three DS schemes is evaluated against this benchmark over a downsampled grid. Fig. 4 shows the evaluated d v of each scheme over the test site. Fig. 5 depicts the empirical probability density function (PDF) of the accumulated discrepancies over the entire test site. The first-and second-order moments of these PDFs provide a measure for the overall bias and dispersion of each method in the estimation of the displacement velocity. These performance indicators are summarized in Table I. As revealed by the comparisons, both the bias and dispersion decrease when more interferograms are exploited for phase and, consequently, deformation estimation. The overall performance of the E-StBAS techniques is observed to be worse than the achievable precision of 1-2 mm/year for the Sentinel-1 data stacks [23], [24], while the exploitation of the    (9), reported for the three DS processing schemes. The first μ-order and second σ -order moments of these Empirical PDFs indicate the overall bias and dispersion of the velocity estimates, respectively. The increase in the number of exploited interferograms improves the overall accuracy of deformation estimation.
full covariance matrix helps EMI in retaining this potential performance.

C. Bias in Displacement Time Series
We further extend the bias analysis by evaluating the discrepancy in the displacement time series. The latter time series is the outcome of spatiotemporal phase unwrapping and the removal of the estimated topographic and atmospheric signal components.
Similar to the previous section and following the method in Section III-B, the displacement bias between each DS scheme and the benchmark PSI result is evaluated for each available time epoch. The displacement bias is, therefore, evaluated per acquisition and over the downsampled grid by (9).
For a compact visualization of the temporal behavior, the box-whisker diagrams are chosen here. Each diagram represents the spatially accumulated bias measures for each acquisition in the stack. The diagram provides a robust presentation of the distribution of univariate data. The box represents the interquartile range, i.e., 50% concentration of the data around the median. The whiskers show the extent of the data excluding the outliers. The outliers are considered as samples falling beyond a threshold defined by the 3× interquartile range from the box boarders. Hence, a portion of the data equivalent to 7× interquartile range is represented by the diagram. The box-whisker diagram concisely depicts the following robust statistical measures.  2) Interquartile Range: Visible from the length of the rectangular box. It is a robust measure for the dispersion of the displacement estimates. 3) Skewness: Inferred from the asymmetry between the upper and lower parts of the box-whisker around the median. Skewness is the third moment of data as a measure for the normality of the distribution of displacement error. A skewed distribution indicates a systematic discrepancy between the PS and the DS. Fig. 6 depicts the box-whiskers for the three compared DS schemes. The average displacement bias of each acquisition is

4) A small periodic trend is observed, which is damped
where more interferograms are exploited. 5) Exploiting all multimaster interferograms and applying EMI for the reconstruction of consistency significantly reduces the bias in displacement estimation.

D. Tracking the Bias to Multilooked Interferograms
Up to this point, we observed an increased bias in the displacement velocity as well as a prevailing presence of a linear physical signal in the displacement time series. The observed signals are absent in the result of single-look PS measurements and fade when exploiting all the interferograms within the time series; therefore, they cannot be justified as the surface topography change. In this section, we seek for the cause of this linear trend in the direct multilooked interferograms and attempt to predict the error budget in the displacement velocity from the phase bias of the mentioned interferograms.
Following the method in Section III-C, the focus here is on the quantification of the interferometric phase bias at varied time lags. Note that according to (10) and (11), φ is required for the evaluation of the phase error. From the analysis in the previous sections, we observed a negligible submillimetric discrepancy between the single-look PS measurements and the estimated phases of EMI approach. Therefore, we may approximate the single-look PS phases by the estimated phases based on the full SCM, i.e., the EMI result, and use the latter as the benchmark φ for the evaluation of the phase bias. Having this benchmark, δ l is evaluated for ∀l ∈ {1, . . . , 10}. Fig. 8(a) and (b) depicts an example of the estimated phase biases for l =1 and l = 10. These time lags correspond to the average temporal baseline of 8 and 78 days, respectively (we resorted to the calculation of the average temporal baseline for each time lag due to the irregular temporal acquisition of the images). The empirical PDF of the spatially accumulated δ l measures for variable time lags are provided in Fig. 8(c) (top), and the bottom figure provides the evolution of the average phase bias over time. From these observation, we conclude that the short temporal interferograms, with a smaller time lag l, are more biased than the longer baselines. This observed trend will assist us in the proposal of a simple physical model for the source of phase biases in Section VI.
We further scale δ l by the temporal baseline and convert the resulted phase bias to a slant range to provide the average predicted range bias δr l in millimeter. Fig. 9 depicts this measure. Inspecting the overall trend of the bias per baseline in Fig. 9(c), the range bias shows the same behavior of attenuation with a temporal baseline. However, the bias is seen to increase slightly up to 24 days, and the attenuation starts only from this time. One hypothesis for this trend is that there are two competing scattering mechanisms, which become dominant at different time scales. Initially, the scatterers that cause the bias are dominant. Over time, these scatterers decorrelate, and the initially weaker but long-term coherent, scatterers slowly prevail, yielding a reduction in the bias. The validation of this hypothesis requires dedicated research, as does any other physical explanation in this article.
Moving to the prediction of the displacement velocity bias explained in Section III-C, the phase bias of the multilooked interferograms is propagated to (12). From the integration of the δ l measure and their proper normalization and conversion, d v is further evaluated for the two processing schemes of E-StBAS5 and E-StBAS10. Note that the two schemes are equivalent to exploiting interferograms with an average temporal baseline of up to 30 and 78 days, respectively. Fig. 10 shows the predicted displacement velocity bias based on the observed phase biases. The empirical PDF of the accumulated d v measures in Fig. 10(c) provides the average bias of these processing schemes over the test site. The predicted bias reads as 5.43 and 3.38 mm/year for E-StBAS5 and 10 respectively. Comparing these biases with the empirical values from the bias  estimates of Table I confirms that the measured biases in the deformation estimates are caused by the interferometric phase biases. Note that the discrepancy between the predicted and the empirical displacement velocity bias is due to two approximations: the underlying simplification of the error propagation scheme (explained in Section III-C) and the approximation of φ by the reconstructed phases of EMI. The sign change between the predicted and empirical values only stems from the difference in the convention for setting the direction of positive displacement.

V. INVESTIGATION ON AN ARID REGION
Here, we aim at observing the fading signal on a different climate and land cover compared with Sicily. For this purpose, we chose an area in Qinghai province, China, located in the Tibetan plateau. The area covers approximately 50 000 km 2 , and Fig. 11 (left) provides the geographical extent of the region as well as its land-cover-classification map. Comparing the latter with Fig. 2, here, 83% of the region is bare areas and 15% is vegetated. The type of vegetation is mostly herbaceous. Fig. 12 provides a comparison between the overall  Table I. Note that in these predictions, positive bias is related to phase excess. Therefore, the positive velocity bias corresponds to a motion away from the satellite, i.e., an apparent subsidence. The sign change between the predicted and empirical values stems from this difference in the convention for setting the direction of positive motion. temporal coherence of the two regions by the empirical PDF of the a posteriori coherence of the two regions (for more details on this phase quality measure and its relation to phase stability, see [15], [22]). This pronounced difference in the a posteriori coherence profiles further assures the disparity of the two studied regions. Using Sentinel-1 data for the time span of October 2014 to November 2019, 104 singlelook complexes (SLCs) comprise our data set. Contrary to Sicily, here, only Sentinel-1 A data are available, rendering the increase in the minimum revisit time to 12 days. Table II provides a comparison between the data set of this test site and Sicily, Italy.
We repeat the experiments of Section IV-D for the current test site. With reference to this section, the phase bias δ l is calculated for the time lags of l ∈ {1, . . . , 10}. These time lags correspond to an average temporal baseline of 18-177 days. Note that for the current data set, the average temporal baseline per time lag is increased compared with the Sicily test site, implying that the number of short temporal baseline interferograms is decreased compared with our previous analysis. The empirical PDF of δ l per time lag is reported in Fig. 13. For the sake of brevity, we only provide the bias map of δ 1 here [ Fig. 14(a)]. Comparing Fig. 13 with Fig. 8(c), the following is observed.  II   COMPARISON OF THE STUDIED TEST SITES. DUE TO IRREGULAR  ACQUISITION OF THE DATA WITHIN THE TIME  SERIES, THE AVERAGE TEMPORAL BASELINE  IS CALCULATED FOR EACH TIME  1) Being an arid region, the overall phase bias is lower than Sicily.
2) The overall bias is in the opposite direction, indicating phase deficit, while in Sicily, a prevalent phase excess was observed.
3) The sign change in bias indicates that here a different physical source causes the fading signal compared with Sicily; 4) Despite these observed discrepancies, the general conclusion of attenuation of bias with a temporal baseline is verified in this test site.
In addition to the prevalent phase deficit in this region, Fig. 14(a) reveals the phase excess over the herbaceous vegetated areas. As evident from this figure, a clear signal is present, which is correlated with the ground features. This figure confirms that the phase discrepancies in the multilooked interferograms are systematic. Recall that in Sicily [ Fig. 8(a) and (b)], the phase bias was unidirectional and the current disparity in the sign of the bias was not observed. This duality of the phase bias behavior reaffirms our argument that modeling and calibration of the fading signal are more intricate than mitigating it by proper phase estimation. We deem the in-depth interpretation of the observed fading signals worthwhile and necessary; however, it falls beyond the scope of this article.
Identical to Section IV-D, d v is further evaluated to predict the minimum expected error on the deformation velocity maps for E-StBAS3, E-StBAS5, and E-StBAS10. The three schemes are equivalent to exploiting interferograms with an average temporal baseline of up to 54, 90, and 177 days, respectively. Fig. 14(b) is presented as an example of the predicted error maps for E-StBAS5. Fig. 15 provides the empirical PDF of the predicted bias for the three schemes. From the comparison of these figures with Fig. 10, the following is observed.
1) The dominant bias shows an apparent motion toward the satellite on the contrary to the case of Sicily. 2) As opposed to Sicily, here, two trends with opposite signs are observable, and the sign of the bias is in agreement with Fig. 14(a). 3) Increasing the number of interferograms, the bias is observed to decrease. 4) The bias is in general lower than the case of Sicily; however, note that the temporal baselines are longer in the Tibetan region. In these comparisons, it should be kept in mind that the Tibetan site has fewer short baseline interferograms as the six-day revisit of Sentinel-1 is not applicable here.
Furthermore, coregistering the error map and the land cover, we observe that, on average, the direction of the bias over the bare areas is toward the satellite, while it is in the opposite direction over the herbaceous vegetation (see Fig. 11). The experiment with this test site manifests the complexity of the physical interpretation of the fading signal and the possible variability of the underlying phenomena.

VI. MODELING OF PHASE BIAS
From the analysis of the two test sites, we observe two general behaviors of the fading signal that is correlated with the land cover. Over the vegetated areas of both the Sicily and Tibetan test sites, an apparent delaying effect is observed, which induces a motion away from the satellite. Over the bare Fig. 14. Repeated experiment with the Tibet data set. (a) δ 1 as the average phase bias of lag 1 interferograms, and the average temporal baseline for these interferograms is equivalent to 18 days. (b) dv as the predicted error in displacement velocity reported for E-StBAS5 scheme, exploiting up to 90-day interferograms. These maps are in the radar coordinate system. In comparison with Sicily, here, two general trends with opposite signs are observed, rendering the physical interpretation of the fading signal more intricate. Note that in these predictions, positive bias is related to phase excess; therefore, positive velocity bias corresponds to a motion away from the satellite. areas of Tibet, on average, the observed motion is toward the satellite. The focus here is on providing a physical explanation for these overall trends. The extensive validation of these hypotheses as well as their fit to the data at hand is subject to the dedicated research. Over the bare areas of the Tibetan test site, the moisture model of [7] can justify the sign of the bias. The typical behavior for soils is to get wet fast and dry slowly. According to [8], the nonlinear moisture model in [7] predicts that the cumulative phase change corresponding to a certain moisture change is larger if the change is observed through many small steps compared with that if it happens suddenly or through larger steps. Therefore, slow drying (moving apparently toward the satellite) should prevail over fast wetting. If there is a bias from the moisture cycle, this should appear as a motion toward the satellite. The mentioned model is confirmed by the observations with real data [7], [9]; however, its fit to the Tibetan data set is not pursued here.
Over the vegetated areas of the two test sites, the sign of the bias is incompatible with the moisture model. Contrary to this, the closure phase behavior over dense vegetation was explained in [9] with the same model presented in [7]. However, in the current two test sites, the observed vegetation is mostly in the form of sparse grass and croplands and herbaceous vegetation. If a moisture bias is present here, it is overshadowed by some other biasing signal. An alternative hypothesis for this land cover is that we are observing biomass growth, with an apparent motion away from the satellite, which introduces extra range delay or phase excess. Based on this hypothesis, a model is proposed here. Let us consider the temporal decorrelation model of [25] i,k = γ 0 exp(−|δt i,k |/τ ) + γ ∞ with γ 0 and γ ∞ , respectively, as a short-term decaying and a long-term persistent coherence of an arbitrary interferogram, δt as its temporal baseline, and τ as the duration of signal correlation. To include the observed systematic phase variation, this model is generalized as follows: The model is comprised of three scattering behaviors: a persistent coherence and two decorrelating parts with varied time constants. Introducing a temporally variable phase term by ρ φ , the latter two decorrelating parts explain the vanishing bias with temporal baselines. The phase of this complex model explains the observed phase bias in an arbitrary interferogram pair.
To examine the biomass model, we chose a single burst over Sicily. Excluding the fast decorrelated areas such as water bodies and the stable PS-like scatterers, the average coherence for each temporal baseline is calculated over this burst. Similarly, the observed average phase bias of δ l ∀l ∈ {1, . . . , 70} is evaluated. Using (14), the observed coherence and phase bias are modeled. The modeling yields a short-temporal scattering of τ 1 = 11 days with γ 1 = 0.18, ρ φ1 = 0.03 rad/day, a longer term scattering of τ 2 = 50 days with γ 2 = 0.25, ρ φ2 = 0.002 rad/day, and a persistent coherence of γ ∞ = 0.13. The result of this model fit is depicted in Fig. 16. As seen from the latter, the model successfully approximates the observed phase bias and coherence. Further validation is required for the use of the proposed model in the calibration of the biasing phase.
Note that in [26], as one of the pioneering works in coherence modeling, a correlation duration of 40 days for the C-band data is predicted. Our three-component model with τ 1 = 11, τ 2 = 50, and an implicit τ 3 ≈ ∞ may be approximated by the mentioned model and its correlation duration of 40 days.
Upon comprehensive study and validation of the two moisture and biomass models, they may be used in the model-based reconstruction of phase consistency, as explained in Section II-B.

VII. DISCUSSION AND RECOMMENDATIONS
In this article, we primarily focused on the observation of a peculiar systematic signal in InSAR. We conclude the following. 1) Multilooked interferograms reveal a short-lived, systematic phase component that cannot be attributed to the atmospheric or surface topography variations. This systematic phase component is absent in the single-look phase observations. 2) Phase bias is larger for the shorter temporal baseline, albeit more coherent, interferograms.
3) The sign, magnitude, and temporal behavior of the fading signal vary at different regions and land covers. This variability renders the calibration of the fading signal intricate. 4) The magnitude of the phase bias in each multilooked interferogram may be deemed small, especially compared with the atmospheric perturbations. However, with the attenuation of atmospheric phase and the propagation of the interferogram phase bias within the time series, the corresponding error for the displacement velocity estimation is alarming. 5) The propagation of even small phase biases in long time series compromises the accuracy of the displacement velocity maps from an achievable submillimetric to centimetric per year level. Studying the effect of the observed signal for interferogram time series, we established the following.
1) The redundancy of the multimaster interferograms is decreased as the result of the present systematic signals. 2) Exploiting the temporal data redundancy in large time series yields a degree of robustness of the phase-retrieval algorithms to such phase errors.
3) The major role in the gained robustness is played by the inclusion of the long-term scatterers in long temporal baseline interferograms. Moreover, we shortly discussed the possibility of modeling the phase biases.
1) A complex-valued decorrelation model is proposed as an approximation of the observed bias over vegetated areas.
2) The proposed decorrelation model helps in the interpretation of the physical phenomena behind the systematic phase component. However, as of present, we lack adequate experimentation and validation of this simplified model for the calibration of the biomass-related phase bias. 3) Two regions were studied and a different behavior of the fading signals was evident. Understanding the physical sources of these signals in the multilooked interferograms is subject to comprehensive research. 4) Due to the lack of comprehensive research thus far, the use of model-free phase linking in the mitigation of phase errors is one of the most viable approaches in safeguarding against the fading signals. As hypothesized, the presence of fading signal compromises the accuracy of InSAR in deformation analysis. This problem is highly exacerbated for Big Data processing. Common practice in such processing is to exploit mostly short temporal baseline interferograms [20], [27], [28]. The logic behind such data exploitation is to reduce the number of interferograms by using the most coherent and, thereby, highest SNR observations in the time series. However, the analysis of this article proves that these interferograms are the most affected by the inconsistent fading signals and, therefore, the least reliable for deformation retrieval. This observation challenges the quality of the short temporal interferograms and warns against solely using these observations. Should the small baseline processing scheme be desired for deformation retrieval, we recommend to safeguard against the impact of systematic phase biases by the following. 1) Excluding the short temporal baseline interferograms and using long baselines to decrease the overall phase error. 2) Predicting phase errors according to Section III-C and (11) to exclude the short temporal baseline interferograms with a significant error relative to a desired accuracy.

3) Predicting deformation error budgets according to
Section III-C and (12); choosing the optimum number of long temporal interferograms (similar to bw parameter) for achieving a desired accuracy. 4) The phase series pertaining to single-look PS could be the benchmark for the mentioned error predictions.
However, according to our experiments, the choice of optimum baseline of interferograms for reliable deformation estimates varies for different types of DS. Furthermore, the additional analysis on the optimum number of interferograms increases the computational burden of processing. In principle, the biasprone DS can be discarded from the deformation analysis at the cost of spatial coverage. The introduced δ l measure and a posteriori coherence provide effective features for the robust detection of error prone DS. For reliable deformation monitoring, we advocate the use of phase linking on full SCM, which is fully data-adaptive. We highlight that there are two issues against the use of phase linking: its robustness over fast decorrelating DS and its computational efficiency. To expand on the former, in case there are no sufficient long-term scatterers to provide stability over long temporal baselines, phase linking might effectively exploit only short temporal baselines, and therefore preserve the bias. We expect such scenarios over forest and dense vegetation. To provide a remedy for the computational cost of phase linking, we have previously proposed EMI in [15] and a sequential estimator in [29]. According to our extensive wide area processing experience, partially reported in [30] and [31], EMI provides a viable efficient solution for phase linking and does not pose a challenge for Big Data processing. The sequential estimator further helps to avoid the redundant computations of phase linking in stream processing of large time series [29] (see Fig. 17). This estimator is based on data compression, and it significantly reduces the number of exploited interferograms while retaining the capability to reconstruct the consistency robustly. The accuracy of this sequential algorithm is studied and proven to retain the millimeter-per-year-level target [30].
Our final recommendation for ensuring the accuracy in deformation monitoring is to introduce a new intermediate Efficient exploitation of Big Data stacks using the sequential estimator [29]. With reference to Fig. 1, the sequential estimator divides the data stack into isolated batches. Between the batches, the lost information is retrieved by the formation of the compressed images (shown by the dot in between the two fully connected networks) and the artificial interferograms depicted with the respective arcs and dark shaded in the SCM. The dashed arcs represent the contribution of the isolated batches in the formation of the compressed images. Following this rationale, the full SCM of Fig. 1(b) is recursively approximated to avoid redundant calculations while retaining the accuracy in phase estimation. product level for InSAR, namely, the reconstructed consistent wrapped phase series using EMI and the sequential estimator. The following holds for the envisioned product: 1) contain the consistent physical signal components such as, but not limited to, atmospheric variations and surface displacements; 2) significantly reduce the interferometric phase bias and stochastic noise, thereby enhancing the reliability of InSAR for deformation retrieval; 3) reduce the amount of interferometric data from the n(n −1)/2 pairwise interferograms within the data stack to a time series of n − 1 higher quality and, optionally, downsampled interferograms; 4) provide a unified product for accurate deformation monitoring to the user community. This article was dedicated to the analysis of C-band SAR. The magnitude of the fading signals and, hence, the corresponding errors are expected to increase with the wavelength. Therefore, larger effects are expected and observed for L-band SAR [8], rendering the introduced phase product even more essential to the L-band.