On the Phase Nonclosure of Multilook SAR Interferogram Triplets

This work explores the properties characterizing the phase nonclosure of multilook (ML) synthetic aperture radar (SAR) interferograms. Specifically, we study the implications of ML phase time incongruences on the generation of ground displacement time series through small baseline (SB) multitemporal InSAR (Mt-InSAR) methods. Our research clarifies how these phase inconsistencies can propagate through a time-redundant network of SB interferograms and contribute, along with phase unwrapping (PhU) errors, to the quality of the generated ground displacement products. Moreover, we analyze the effects of short-lived phase bias signals that could happen in sequences of SB interferograms and propose a strategy for their mitigation. The developed methods have been tested using both simulated and real SAR data. The latter were collected by the Sentinel-1 A/B (C-band) sensors over the study areas of Nevada, USA, and Sicily, Italy.

In our work, we study the properties of phase nonclosure among sets of time-redundant networks of ML SAR interferograms to characterize their effects on the ground displacement time series retrieved using Mt-InSAR SB algorithms. They exploit time-redundant, reduced networks of interferograms with small perpendicular and temporal baselines [19], [21], [31], [39], [40], [41]. We show that a set of stable, coherent pixels at the ML scale, can be suitably identified by analyzing the phase triplets obtained from the selected networks of SB interferograms [42], [43]. Furthermore, a method for estimating and reducing the phase biases into sequences of ML interferograms, relying on the exclusive use of phase triplets, is developed. Experiments were carried out on simulated data and two sets of SAR images collected by the European Copernicus Sentinel-1/A-B (S-1/A-B) sensors over Nevada, USA, and Sicily, Italy. The proposed investigation demonstrates the validity of the developed phase bias mitigation method.
This article is organized as follows. Section II presents the theoretical background of the interferograms' phase nonclosures. Section III proposes a method for phase bias reduction in sequences of ML SB interferograms. Section IV shows some simulations. Section V studies the implications This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ of noncompensated phase closures for generating interferometric SAR (InSAR) products through SB-oriented Mt-InSAR processing chains. Experimental results are shown in Section VI. Conclusions are finally addressed in Section VII.

II. PHASE TRIPLET'S CLOSURE PROPERTIES
Let us consider a set of N SAR images collected at the ordered times t = [t 0 , t 1 , . . . , t N −1 ] T and co-registered to a common reference geometry, and let = [φ 0 , φ 1 , . . . , φ N −1 ] T be the (unknown) vector of the full phases (i.e., not restricted to the [−π, π] range) associated with every single SAR image. Given three interferometric SAR data pairs computed from three generic SAR images collected at times t h , t k , t q (see Fig. 1), the following relation holds: where φ n,m = φ m − φ n is the phase difference between the single-look (SL) SAR images at times t m and t n . However, the relation (1) is not applicable when ML interferograms are considered because they involve the estimation of averaged phase values that are independently computed (interferogram by interferogram) over a group of neighboring SAR pixels and they are generally time-inconsistent. In this case, the phase turns out to be a rotational field, see [2] and [44] where φ ML n,m = ϕ ML n,m + 2πU n,m is the phase of the generic (n, m) ML SAR interferogram, ϕ ML n,m is the relevant (wrapped) ML interferometric phase, U n,m is the number of correct (generally unknown) 2π-cycles of the (n, m) unwrapped ML interferogram, and L is the number of looks. Note that j = √ −1 and is the phase extraction operator; also, is the group of averaged SAR pixels at the SL scale used to compute the ML interferograms. The (time) rotational phase field in (2) is called a nonclosure phase triplet. The excess phase that prevents the closure of the triplet in (2) has different contributions related both to random and systematic sources [1], [2], [34], [44], [45], [46]. The nonclosure phase is a time rotational (i.e., a time-inconsistent) field.
Note that given N SAR images, among the whole possible N(N − 1)/2 SAR interferograms that can be generated, only N − 1 are independent. In contrast, the others (N − 1)(N − 2)/2 can be calculated from linear combinations of the previous ones. Multitemporal SB InSAR techniques [18], [47] rely on processing a set of SB ML interferograms M ≤ N(N − 1)/2, which are typically selected by imposing thresholds on the maximum allowed temporal and perpendicular baselines of the SB interferograms. Given such a set of M SB SAR data pairs, a number, say , of phase triplets could then be identified. Specifically, it is worth noting that, considering only a set of SB SAR data pairs, could be noticeably smaller than the maximum number of triplets that can be formed with N SAR images, which is equal to [2]. We observe that the zth ML phase triplet can be expressed as where φ tr z and ϕ tr z = W [φ tr z ] = 0 are the zth unwrapped and wrapped nonclosure phase triplet, respectively, U z is the composite 2π-integer multiple of the zth triplet, and h(z), k(z), and q(z) are the three epochs of the zth triplet. Note also that W (·) is the operator that wraps out the phase into the [−π, π] range. It can be demonstrated that (4) where the last term on the right-hand side of (4) is the spurious phase cycle arising from the observation that the zth nonclosure phase triplet could exceed [−π, π].
Remarkably, the integer terms in (4) are estimated during the space-time PhU operations (see [22]). Nonetheless, the wrapped phase contribution for the considered triplet ϕ tr z persists, even when PhU operations are perfectly accomplished. If not adequately compensated for, the data vector tr = [ϕ tr 0 , ϕ tr 1 , . . . , ϕ tr −1 ] T might influence the quality of the Mt-InSAR results (i.e., the ground deformation time series and the relevant mean ground deformation velocity maps) obtained after inverting the sequence of unwrapped ML interferograms

A. ML Speckle Noise Model
The mentioned phase triplets' inconsistencies tr could arise, for instance, when in the ML averaging box, different populations of scatterers, characterized by independent phase histories, interfere with one another [9]. Indeed, both the spatial and the temporal inconsistency in the ML averaging window can give rise to nonzero closure phases [38]. To investigate the origin of these phase inconsistencies, let us consider the nth phase triplet, namely ϕ tr n , which involves the three SAR images collected at times (t h , t k , t q ), see Fig. 1. First, let us focus on the (t h , t k ) SAR data pair and consider a scenario with F independent populations of scatterers in the averaging window. The generic f th scatterers' population can be assumed, for the sake of simplicity, to have only two independent phase contributions: 1  can be adopted where the symbol · L denotes the spatial average operation computed over L samples, in which f is the partition of the group of pixels of the averaging window into F independent subgroups, representing the relevant scatterers' populations, with every single set f with L f looks. We remark that is the true, expected signal component related to the f th population, with ρ f being the corresponding coherence value of the scatterer's family. Moreover, n m, f is a multiplicative noise component, N c, f ∼ = 1 − 1/(8L f ), andz L f is defined in [48], where E[n m, f − N c, fz L f ] = 0 and var Of course, in order to describe the characteristics of the scatterers' families involved in the ML window, other speckle noise models, such as the one proposed in [38], might be applied.
Accordingly, considering the speckle noise model described by (5), the ML (expected) phase related to the (h, k) interferometric SAR data pair can be expressed as where ζ h,k is a resulting zero-mean additive noise phase term. As a first approximation, we can assume that the model in (6) depends only on the temporal baseline of the Concerning the temporal decorrelation models that describe how the coherence depends on the interferometric temporal baseline, interested readers are referred to the literature, see, for instance, [49], [50], and [51].
In Section III, we then relax these hypotheses to consider a more general time-variant case. Space-time, physical and statistical properties of local, inherent signals that contribute to the systematic phase biased signal can be found in the literature. For instance, readers are referred to [9], [11], [12], [14], [15], [16], [52], and [53] to have a comprehensive analysis of the models adopted for the characterization of soil moisture content variations and its impact on InSAR investigations. We would like to remark that, in our work, we do not want to discriminate one another the different inherent, local signals φ ML, f h,k , ∀ f that contribute to the phase biased signal nor the time series of the inherent phase contributions. Conversely, we want to estimate and mitigate the effect of the "composite" global systematic phase biased components ϕ bias on the ground deformation products, as obtained using Mt-InSAR algorithms (e.g., [18], [19], [39], [47], [54]).
Experimental results have evidenced that the systematic phase bias ϕ bias is a short-lived signal that rapidly decays as the temporal baseline increases [34], [36], [37], [38], [55]. For small values of t h,k = |t k − t h |, we can locally expand the (full) phase bias related to the (h, k) interferogram φ bias h,k as φ bias where v is a constant decay phase velocity factor and v(t h,k ) is a temporal-baseline-dependent phase velocity difference term. Using (6) and (7), the nth phase triplet ϕ tr n can be expressed as follows: where ζ h,k,q is the zero-mean random noise phase term related to the (t h , t k , t q ) interferometric triplet. From (8), the (wrapped) triplet systematic phase bias can be expressed as , which solely depends on the velocity difference terms v(t i, j ) i, j = h, k, q and is insensitive to the mean constant decay velocity factor v. Because the phase bias is a shortlived signal, with significant phase rates only at very small temporal baselines and negligible rates at medium-to-long baselines, it is reasonably assumed that |φ tr n | ≤ π (see [9]). On the other hand, using polygons of interferograms (see [37], [38]), the probability that the absolute value of the nonclosure phases could exceed π increases as the number of polygon sides increases. This probability is assumed low in [37], and however, in the case that these long phase loops exceed π in moduli, there would be a corrupt estimate and compensation of the bias. On the contrary, in [38], the phaseunwrapping operations are performed on the nonclosure phase loops. In this case, some (unavoidable) PhU errors committed could potentially degrade, to some extent, the reliability of the estimation and compensation of the phase bias.

B. Statistical Properties of the Phase Triplets
Given a set of N SAR images, the selected M SB SAR interferograms can be arranged to form triplets. For every radar pixel, the vector of the (wrapped) phase triplet tr = [ϕ tr 0 , ϕ tr 1 , . . . , ϕ tr −1 ] T is a multisample random circular data vector whose elements have different statistics, depending on the geometrical characteristics of the different families of triplets that could be formed. For instance, if we consider the family of interferometric triplets made by three Sentinel-1 SAR interferograms with temporal baselines of 6, 6, and 12 days, the corresponding systematic phase bias signal φ bias 12. In general, when the family of interferometric triplets with temporal baselines of t h,k = t k,q and t h,q = 2t h,k is concerned, the relevant systematic phase bias equals: Therefore, under the simplified hypothesis that the model in (6) is invariant with time, the set of triplets, namely, Tr ≡ {Tr n } n=1 , can be partitioned as Tr = X χ =1 Tr χ , where X is the number of different homogeneous families of interferometric triplets. 2 Let us now consider the generic family Tr χ , composed of χ elements, and the phase vector tr χ = [ϕ tr χ,0 , ϕ tr χ,1 , . . . , ϕ tr χ, χ −1 ] T . It can be demonstrated that the elements of the vector tr χ are Von Mises-distributed VM(μ χ , κ χ ), with an averaged phase value μ χ and a concentration parameter κ χ [56] where χ is the (sample) mean resultant length of the χth family of phase triplets Equations (9)- (11) show that the systematic bias of the considered family of triplets is the mean direction of the phase triplets' distribution μ χ , whereas the mean resultant length χ gives us a measure of the spread of the phase triplets with respect to μ χ , accounting for the zero-mean random noise contribution ζ h,k,q . The first-order statistics of the amplitude and phase of the directional random data 2 Note that the hypothesis on the homogeneity of the phase triplet families is subsequently relaxed in Section III. Appendix I. It can be demonstrated that the standard deviation of the χth phase triplet family is given by The mean resultant length of the combined multisample data vector tr = [ϕ tr 0 , ϕ tr 1 , . . . , ϕ tr −1 ] T is eventually given by which is referred in this work to as the triangular coherence of the identified network of SB ML interferograms. We remark that triang can be seen as an equivalent coherence [57] and gets a direct measure of the nonclosure (overall) phase triplets' dispersion to its mean value. The triangular coherence can be used to select a group of reliable, coherent SAR pixels at the ML scale. Indeed, the triangular coherence gets a measure of the noise level that affects the selected set of ML interferograms. In particular, the mean resultant length of the phase triplets Mardia and Jupp's book [56]). Accordingly, the variance of this estimator drastically decreases as the number of triplets and the triangular coherence value triang increase. The group of coherent SAR pixels is thus straightforwardly identified by simply imposing a threshold γ triang on the minimum triangular coherence value of the analyzed pixel, as ≡ {P : triang (P) ≥ γ triang }.

III. PHASE BIAS ESTIMATION AND MITIGATION
This section presents a method to estimate the phase bias affecting a sequence of SB ML interferograms characterized by the maximum temporal baseline t max .

A. Time-Invariant Case
Let us first assume that the model of (6) Under this hypothesis, the expansion given by (7) is suitable. We remark that the temporal baseline of an interferogram is a multiple of the sampled temporal revisit time of the used SAR sensor, namely, δ, which, for instance, is equal to six days for the twin constellations of Sentinel-1 A/B sensors. Accordingly, if we consider two generic interferometric SAR data pairs with temporal baselines (λ − 1)δ and λδ, with λ ∈ Z, (7) particularizes as Comparing (14) and (15), the following iterative relation is derived:

Equation (8) defines a system of linear equations with respect to the (unknown) phase bias velocity differences
that can be solved in the least-squares (LS) sense as V = Z † · tr , where Z † is the pseudoinverse of the matrix Z; the symbol · stands for the matrix multiplication (rows by columns) operator. Then, the estimates V are used to iteratively compute the phase biases at the different temporal baselines t λ = λδ, ∀λ = 1, 2, . . . , t max /δ through (16) using the initial condition φ bias (t max ) = 0. As claimed in several independent investigations [34], [36], [37], [38], the phase bias is a signal that rapidly decays as the temporal baseline increases. However, it is not guaranteed that the maximum temporal baseline t max of the selected set of SB interferograms is large enough to assume that φ bias (t max ) = 0. A strategy to understand whether the maximum temporal baseline t max is adequate for the phase bias correction is to compute, and for every single SAR pixel of the scene, If ϒ is larger than a given tolerance ϒ toll (e.g., ϒ ≥ 10 −4 ), some additional interferograms with longer baselines must then be added. It is worth remarking that the additional long-baseline ML interferograms are exclusively used to compute the phase bias. However, they are not exploited to generate the ground deformation time series using an SB-oriented algorithm [18], [19], [39], [47], [54].

B. Time-Variant Case
Let us now assume that the model of (6) . In this case, the strategy described in Section III-A can be specialized by locally applying it to single time windows, by dynamically selecting and using only subgroups of triplets' families that encompass the selected time window to estimate the phase bias. To describe the developed method, let us focus on the generic i th ML SAR interferogram that spans the time window between the times t h and t k of duration m h,k δ. The phase bias estimate related to this i th interferogram is carried out by using (17) considering only the triplets {α, β, γ } (see Fig. 2) that satisfy both the following conditions. 1) At least one of the three arcs {α, β, γ } of the given triplet must wholly be included or include the reference time window [t h , t k ]. 2) For all the three arcs: the generic arc of the triplet is the reference time window [t h , t k ] or its duration is different from the reference time window duration.
The block diagram of the proposed phase bias correction method is shown in Fig. 3.

C. Compensation of the Random Phase Signal Components
Once the systematic phase bias components are estimated and compensated for, however, some uncompensated timeinconsistent random phase noise contributions, which lead to phase nonclosures, can persist. They can effectively be compensated for using the methods described in [1], [3], [4], [46], [58], [59], [60], and [61]. In particular, the application of the first step of the extended minimum cost flow (EMCF)-based processing chain [3], which is fully detailed in [57] and here referred to as enhanced multitemporal noisefiltering (E-Mt-InSAR) algorithm, allows one to obtain a set of optimized, fully time-consistent set of ML SAR interferograms. More specifically, for every SAR pixel of the ML grid, the E-Mt-InSAR method is based on searching for the (unknown) phase vector of the wrapped phases related to the available N SAR acquisitions that minimize the (weighted) circular variance of the random phase vector representative of the difference between the original and the optimized interferograms reconstructed from the computed (wrapped) phases associated with every SAR acquisition (i.e., the residual phases). It is worth remarking that the (weighted) circular variance gets a measure of the dispersion of the residual phases about their (weighted) mean direction. However, no constraint is imposed about the mean (weighted) direction of the residual phases, and the estimator is not able to adequately discriminate between the short-lived systematic phase contributions and the zero-mean random phases associated with the ML SB interferograms, see also the experimental results shown in Appendix II. An effective strategy that could be adopted is, first, to compensate/mitigate the systematic phase bias components using the algorithms described in Sections III-A and III-B and then apply the E-Mt-InSAR noise-filtering algorithm to the set of compensated ML interferograms {φ ML i } M i=1 . Accordingly, using this strategy, both the systematic and random noise phase contributions can be adjusted.

D. Role of PhU Errors
The presented analyses do not consider the effects of biased time-inconsistent PhU errors. As a matter of fact, once the ML interferograms are unwrapped, some of the observed discrepancies in the generated InSAR products can also be due to time-polarized PhU errors, e.g., 2π-multiples PhU errors that are superimposed on phase triplets and are responsible for time incongruences that can propagate through the selected network of SB interferograms. These effects arise when the SB interferograms are unwrapped independently [62], [63], [64], [65] and/or using some hybrid 2-D + 2-D space-time PhU methods [22], [26], [28] that do not ensure that the PhU solution is time-irrotational. These effects are appropriately considered by the temporal coherence factor, which is computed after the SB inversion of the unwrapped interferograms, as detailed in Section V. Hence, the temporal coherence value is used to detect those pixels that are more affected by time-inconsistent phase artifacts and exclude them from the subsequent analyses.

IV. SIMULATION
The developed phase bias correction method was first tested in a controlled environment by running some simulations. Specifically, we have considered the following cases: 1) the model (6) is time-invariant (see Section III-A) and 2) the model (6) is time-variant where the nonstationary phase correction method described in Section III-B is applied. For both cases, we have also considered the effects of the decorrelation noise that has been introduced in the test environment.

A. Time-Invariant Case
In this first case, the assumption is retained that the phase bias depends only on the time span of the considered InSAR where t h and t k are the two generic time acquisitions. The simulation has been carried out considering the time distribution and the SB interferometric network settings related to the Nevada case-study area (see Section V). We considered N = 65 SAR acquisitions and M = 895 SB InSAR data pairs selected by imposing a maximum temporal baseline t max = 96 days. The interferograms' temporal baselines are sampled as: = 6305 InSAR triplets have been identified. By referring to the phase model described by (6), we assumed the presence of F = 2 independent populations of scatterers in the averaging ML spatial window and for the sake of simplicity that ϕ ML,defo h,k = 0 ∀h, k. First, a free-noise scenario was assumed, i.e., ζ h,k = 0 ∀h, k. In this case, the adopted InSAR phase model (6) particularizes as: For the first family of scatterers, we set φ ML,1 h,k = 0, ρ 1 = 1, and ψ 1 = 1. Differently, for the second family of scatterers, the following values are set: φ ML,2 h,k = ξ t h,k , ρ 2 = e −α t h,k , and ψ 2 = 1, where α = (11/2)(1/t max ) and ξ = qπ(1/t max ). Fig. 4 shows the plots of simulated inherent (independent) phase term for the second family of scatterers versus the interferograms temporal baseline, considering the two cases with q = 1 and 4. In this round of simulations, we have assumed that the average phase of the second family varies linearly with time and ρ 2 = e −α t h,k in accordance with the temporal decorrelation models proposed in the literature (see [50]). In such a way, the made hypothesis that model (6) is time-invariant is valid. Based on the parameters listed above, the adopted InSAR phase model becomes: The simulated phases, which are plotted versus the temporal baseline (blue lines) in Fig. 5(a) and (b), correspond to a maximum (simulated) phase bias velocity of roughly 2.1 and 8.7 cm/year, related to InSAR data pairs with a temporal baseline of t = 6 days and assuming a wavelength of 5.546 cm (i.e., that of Sentinel-1 A/B), for the simulations with q = 1 and 4, respectively. We have applied the phase bias estimation method described in Section III-A to simulated phases. We want to remark that, even with q = 4, the simulated phases do not lead to phase triplet ambiguities, that is, φ tr ∈ [−π, π]. For instance, if we consider a triplet with side lengths of 6, 6, and 12 days, the excess phase 2φ 6 days − φ 12 days is of about 0.1787 rad. The same does not happen if we use for phase estimation, instead of triplets, closed loops forming polygons in the temporal/perpendicular baseline plane with several SB arcs. Indeed, if we consider a polygon formed by A = 10 arcs of side length 6 days and one arc of side 6 A days, the excess phase Aφ 6 days − φ 60 days is of about 3.1909 rad, which is outside the range [−π, π]. In this case, the polygons must first be unwrapped, and this operation could introduce some undesired artifacts in the phase bias estimates.
The estimated phases are shown with red triangles in Fig. 5(a) and (b). The results show that when the model in (6) is time-invariant, and in the absence of noise, the method perfectly reconstructs the InSAR biased phases, also when they exhibit a sign change [see Fig. 5(b)]. Fig. 5(c) and (d) shows the results of the simulations obtained by adding to simulated phases the phase noise components ζ h,k = 0, for the two cases with q = 1 and 4. ML noise signals have been simulated (e.g., see the statistics shown in [48]) considering a coherence value of 0.35 and an equivalent look number (ELN) of 80. The results demonstrate that the proposed method is robust to decorrelation noise artifacts.

B. Time-Variant Case
At this stage, we study what happens in the more general case that the model (6) The phase estimation method described in Section III-B is applied in this case. Fig. 6(a) [Fig. 6(b)] shows the plots of the estimated interferograms biased phases (biased phase velocity) versus the InSAR temporal baseline for six selected temporal windows (dashed lines), which corresponds to the InSAR data pairs listed in Table I. In this first case, we have considered ξ = qπ(1/t max ) with q = 1. The same has been repeated considering q = 4. The plots of the estimated biased phases (biased velocity phases) versus the InSAR temporal baselines are shown in Fig. 6(c) [Fig. 6(d)]. Black solid lines in the plots of Fig. 6 represent the estimated phases retrieved using  the time-invariant algorithm of Section III-A. We can observe that the time-invariant method cannot follow the time-variant fluctuations of the adopted model. Fig. 7(a) and (b) shows the comparison between the simulated and the estimated biased phases for both the scenarios with ξ = qπ(1/t max ) using q = 1 and 4. Simulated biased phases (blue squares) accounting for the time-variant η h,k terms have been shown as well as the relevant estimates (red triangles) obtained by applying the method described in Section III-B. The achieved results evidence that the biased phases have correctly been estimated using the developed method, which can track the different fluctuations of the phase bias in the interferograms that belong to every group of temporal baselines, especially for those with small and very small baselines. Also, in this case, we have evaluated the effects of the noise on the phase estimates. The results are shown in Fig. 7(c) and (d).

V. SB MT-INSAR METHODS
In this work, we refer to a unified representation of the SB algorithms. Indeed, the different implementations of the SB methods proposed in the literature have individual peculiarities; however, they can almost be unified [47], [66] considering that they solve a linear optimization problem that relates the vector of the (known) unwrapped ML SB interferograms, namely, = [φ 0 , φ 1 , . . . , φ M−1 ] T , to a model of Q unknown parameters of the ground deformation M d = [M d,0 , M d,1 , . . . , M d,Q−1 ] T . The adopted unified representation of the SB linear transformation is given as follows: where B ∈ R M×Q is the design matrix of the considered linear transformation. Precisely, for the implementation of the Small-BAseline Subset (SBAS) algorithm [18], the model parameters represent the velocities between time-consecutive SAR acquisitions, namely, )] T , and the mathematical expression of the design matrix B is that detailed in [18] and [66]. Once (18) is solved, the obtained model vectorv is time-conservative, that is to say, there exists a unique vectorˆ that satisfies the following relation B ·v = A ·ˆ and the Euclidean twonorm A ·ˆ − 2 is minimal, where A is the incidence matrix of the SB network nonplanar graph (which is equivalent to the discrete gradient operator) andˆ is the vector of the (unwrapped) phases associated with every single SAR image [66]. Here, we explicitly address some critical properties of the vector . First, we observe that the mth (wrapped) phase ϕ m can be decomposed into the sum of a time-consistent ϕ m and a time-inconsistent ϕ m phase component [67]. Therefore, where U m is the correct (unknown) ambiguity number associated with the mth ML interferogram, which also considers the ambiguity number that could arise from the sum of the wrapped phase terms ϕ m + ϕ m . Of great relevance for the generation of the ground displacement products through SB methods are the time-inconsistent phase terms φ m , which (as said) are responsible for nonclosure phase triplets (see Section III). This section shows how the residuals of the LS solution of (18), namely, res = B ·M d − , are related to the nonclosure phase triplets. Moreover, we summarize some basic properties of the temporal coherence quality factor [22] to provide a quantitative estimate of the PhU errors committed after the SBAS inversion [3], [18], [22], [68]. We can express it as where φ res m is the mth element of the residual vector res = (B ·M d − ) ∈ R M . Accordingly, the temporal coherence temp is the mean resultant length of the circular phase data vector res . We can observe that where ∈ Z ×M is the discrete curl operator related to the SB network graph, which is a sparse matrix composed of (number of triplets) rows and M (number of interferograms) columns. Considering (19), (21) can be manipulated as follows: · res = · A ·ˆ − · + + 2πÛ (22) whereÛ = U is the vector of ambiguity numbers estimated through PhU operations, which could generally differ from the correct (unknown) vector U. Note also that · A ·ˆ = 0 and · = 0 being the former the discrete counterpart of the ∇ × ∇ operator and a time-conservative vector. Furthermore, if the estimated ambiguity number vector is not correct, namely,Û = U + U + U , where U and U , respectively, are the time-correlated and time-uncorrelated ambiguity number errors, the following relation holds: because the (true) phase ambiguity cycles and the time correlated ambiguity number errors are time-conservative, namely, 2π( · U) = 0 and 2π( · U ) = 0. If we assume that no PhU errors were committed, i.e., U = 0, (22) particularizes as · res = − tr (24) where tr = W ( · ) ∈ R is the vector of temporal nonclosure phase triplets, see (3). In (24), we have also assumed that phase triplets are within [−π, π]. Equation (24) relates the LS residuals after SB inversion res and the nonclosure phase triplets tr , making it evident that noncompensated phase triplets correspond to errors after the SB inversion and have an impact on the temporal coherence. Finally, if we relax the hypothesis that no PhU errors were committed, we have · res = − tr −2π ·U . Therefore, the phase residuals increase due to PhU errors and the corresponding values of the temporal coherence decrease. Note that even with a perfect compensation of the systematic and random (wrapped) nonclosure phases, some significant phase residuals could still arise after applying the SB inversion to unwrapped interferograms if substantial PhU errors are present.

A. Nevada, USA
The area is located in the Monte Cristo Range, approximately 38-km west-northwest of Tonopah. On 15 May 2020, a severe 6.5-Mw earthquake struck the west side of the selected area, fortunately without casualties. The 2020 Monte Cristo Range earthquake represents one of the strongest earthquakes in Nevada in the past 100 years, precisely the strongest since 1954 [69]. Starting from the available SAR, ML SB interferograms were generated, by considering a temporal baseline threshold of 96 days. For the interferogram generation, we adopted an ML factor of 4 and 20 pixels for the azimuth and range directions, respectively. The one-arcsec shuttle radar topography mission (SRTM) digital elevation model (DEM) of the scene and precise orbits of the Sentinel-1 A/B satellites was used to compute the topographic phase and flatten the interferograms. Fig. 8 shows the triangular coherence map of the investigated area, calculated considering = 6305 phase triplets computed over a selected network of 895 SB ML interferograms. Subsequently, we have imposed the threshold γ triang = 0.50 and exclusively processed the selected coherent pixels, which are almost three million. The phase related to the group of coherent SAR pixels was unwrapped through the minimum cost flow (MCF) solver [62], and then, the SBAS [18], [68] inversion was applied. As a result, the line-ofsight (LOS) mean displacement ground velocity map of the area was obtained, which is shown in Fig. 9, where the effects of the strong rupture ascribable to the Mw 6.5 Monte Cristo Range earthquake are evident. Indeed, on the left part of the scene, the areas straddling the fault line, showing opposite LOS deformation trends, are evidenced with an absolute value of maximum LOS mean displacement velocity of about 20 cm/year. We specify that the mean displacement velocity map, shown in Fig. 9, has been computed from the interferometric network composed of 895 ML InSAR pairs with a maximum temporal baseline of 96 days, without having applied any phase bias correction method.
Next, we have applied the phase bias estimation methods described in Sections III-A and III-B (for the timeinvariant and time-variant cases) to the Nevada SAR data set using networks of SB interferograms, with a progressively reduced maximum temporal baseline (i.e., from 96 to 12 days). The selected SB ML interferograms were independently  corrected, unwrapped, and inverted through the SBAS technique [18]. Then, for every SAR pixel, the values of the mean ground displacement velocity and the temporal coherence have been computed. Note that, for the subsequent analyses, the atmospheric phase screen (APS) and the residual topography components were not compensated for. Fig. 10 shows the plot of the (average) absolute values of the difference between the computed ground deformation velocities using SB networks at given temporal baseline thresholds with respect to those achieved by considering as a reference a maximum temporal baseline of 96 days. Three groups of coherent SAR pixels have been identified with temporal coherence values greater than 0.7 in Fig. 10(a), 0.90 in Fig. 10(b), and 0.98 in Fig. 10(c). For every group, we plotted the absolute velocity biases (mm/year) for the original ML interferograms (blue line), the ML interferograms compensated using the time-invariant phase bias estimation method (red line), and those compensated with the time-variant method (black line). The results show that with high temporal coherence values, the developed phase bias estimation methods reveal effective, and the time-variant algorithm has a better performance than the time-invariant one, especially at very SBs. Both the phase bias estimation methods can reduce the effects of the ground displacement velocity biases of the ML interferograms with respect to what happens using the original, uncompensated interferograms. Fig. 11 shows the maps of the ground deformation velocity differences between the results at temporal baselines of 12 and 96 days, where only pixels larger than the given values of the Fig. 12. Average absolute values, related to the Sicily Island test case, of the difference between the computed ground deformation velocities using SB networks at given temporal baseline thresholds with respect to those achieved considering a maximum temporal baseline of 96 days. No correction (blue), time-invariant correction (red), and the time-variant correction (black). Temporal coherence is greater than 0.95. temporal coherence are mapped. More specifically, Fig. 11(a)-(c) shows the ground deformation velocity bias map for the original interferograms, those compensated using the time-invariant phase bias estimation method, and finally those obtained using the time-variant ones, respectively. The maps portray only SAR pixels with temporal coherence values larger than 0.7. Fig. 11(d)-(f) is the same as in Fig. 11(a)-(c) but shows only SAR pixels with a temporal coherence larger than 0.9, and Fig. 11(g)-(i) shows only those with a temporal coherence larger than 0.98. As demonstrated in Section V, the temporal coherence is sensitive both to time-inconsistent phase terms and PhU errors. No substantial ground deformation velocity biases can be appreciated in the maps of Fig. 11(i) as well as in the relevant plot (black line) of Fig. 10(c), whereas with the same temporal coherence threshold, they were evident. This demonstrates that the developed methods are effective. The areas with lower temporal coherence values are more affected by time-inconsistent PhU errors and the latter can be responsible for some of the observed ground deformation velocity differences, see, for instance, the area highlighted by the black rectangle (dashed line) in Fig. 11(c).

B. Sicily, Italy
As a second case study, we have considered a dataset of Sentinel-1 SAR data gathered over the area of Sicily, Italy. The investigated area is home to one of the most important active volcanoes of the world, the Mount Etna volcano, which has been continuously active for the last three decades. Several InSAR investigations have been carried out in recent years for the characterization of the magmatic sources beneath the volcano cones and for the study of the superficial effects of volcanic eruptions and relevant local seismicity [42], [70], [71]. First, starting from the available SAR images, we have generated 1000 SB ML interferograms with a maximum temporal baseline of 96 days. The interferograms were flattened using a three-arcsec DEM of the Sicily Island and precise orbital information. The interferograms were multilooked with a box of 4 × 20 (azimuth and range) SAR pixels and unwrapped [62]. We have identified different networks of SB InSAR data pairs by imposing progressively decreased maximum temporal baseline thresholds, from 96 to 6 days, and we cross-compared the values of the ground displacement velocities taking as a reference the results achieved using a maximum InSAR temporal baseline of 96 days. Also, in this case, to exclude their effects on the cross-comparison analysis of the InSAR products, we did not filter out the APS and the residual topographic components from the generated InSAR products. It is worth remarking that in practical cases, other strategies could also be used to identify subsets of SB ML interferograms and use them to generate the ground displacement products [3], [54], [72]. In this work, we want to emphasize the performance of the developed methods in the general case that conventional ML SB interferograms are used, which are straightforwardly selected by imposing a threshold on their maximum allowed temporal baseline. Fig. 12 shows the plot of the (absolute) SBAS-driven ground displacement velocity differences versus the maximum temporal baseline of the used network of SB interferograms, taken as a reference the ground displacement velocity obtained using 96 days as the maximum temporal baseline. The three plots are related to the original interferograms (blue line), the interferograms compensated using a time-invariant model (red line), and those obtained considering a time-variant model (black line). To exclude that (at most) the ground deformation velocities difference could be due to time-inconsistent PhU errors, 3 only SAR pixels with a very high value of the temporal coherence (i.e., larger than 0.95) have been considered. The results show that the developed phase bias estimation methods are capable of considerably reducing the effects related to systematic phase biases in the generated ground displacement maps. Note that in this case, we have also considered the extreme case when only six-day interferograms are used to compute the ground displacement time-series. We want to remark that any potential PhU error on such six-day interferograms has a severe impact because it is time-integrated and appears as an undesired jump in the obtained ground deformation time-series. Fig. 13 shows the maps of mean displacement velocity differences between the results obtained from 96 to 6 days, for the original one in Fig. 13(a)-(c), and those corrected by applying the developed phase estimate timevariant method in Fig. 13(d)-(f). The maps are obtained by considering different temporal coherence thresholds to make it evident the spatial distribution of SAR pixels that are more affected by PhU errors.

VII. CONCLUSION
The properties of nonclosure phase triplets in sequences of ML SB interferograms have been addressed in this research study. We have shown that phase triplets and phase residuals, arising from the LS inversion of the unwrapped SB interferograms, are linked to one another. We have discussed the implications of ML nonclosure phase triplets on the claimed fading signals [34] into the SB-driven mean ground displacement maps obtained using very SB SAR interferograms. Two methods that could be adopted to reduce the effects of noncompensated systematic phase biases in the generated ground displacement time series have been proposed, considering both a time-invariant and time-variant model for the systematic phase bias signals. The results obtained by applying the developed algorithms to simulated and real SAR data have also been presented and discussed. The methods rely on the exclusive use of triplets of ML interferograms, instead of using polygons of ML SB interferograms, and assume that the systematic phase bias rapidly decays over time, and it is negligible at large temporal baselines. We also proposed a recursive strategy to understand at which temporal baseline we could reasonably assume that the phase bias is almost zero. It consists of adding long-baseline ML interferograms that are, however, only used to implement the phase bias correction step, and they are subsequently discarded for the generation of ground displacement time series. Yet, we have shown that nonclosure phase triplets can be used to effectively identify a group of coherent SAR pixels at the ML scale over which ground deformation InSAR products can be computed.
Globally, our research study demonstrates that, under proper hypotheses, with the preliminarily compensation/reduction of the systematic phase errors in the sequence of original wrapped interferograms, the SB Mt-InSAR methods reveal effective for the generation of ground displacement products, also with new-generation SAR systems characterized by frequent repetition times (on a weekly basis or less). Nevertheless, further developments are still required to understand the ultimate consequences of the very small baselines and nonclosure phase triplets for the detection and analysis of soil properties and other "hidden" signals that can be extracted from sequences of SB ML InSAR interferograms.
As a final remark, we want to highlight that for the presented investigations, we have intentionally not used 3-D [73], [74], [75] or hybrid 2-D + 2-D PhU methods [22], [26], [28]. Conversely, we have unwrapped every single interferogram independently [62]. Also, we have selected the SB interferograms to be inverted by simply imposing a threshold on the maximum allowed temporal baseline of the interferograms. We followed this strategy precisely to emphasize the role played by the systematic phase artifacts and PhU errors.

APPENDIX I
Let us consider the χ -length random directional data vector [ϕ tr χ 0 , ϕ tr χ 1 , . . . , ϕ tr χ χ −1 ] T which are assumed to be Von-Mises distributed VM(μ χ , κ χ ). From the work [56], we can consider its resultant vector, which is the following random variable: where R χ and χ are the mean resultant length and the mean direction of the considered directional data vector, respectively. Let us also ρ χ and μ χ assume be the average (sampled) values of the mean resultant length and the mean direction, respectively. For large samples, e.g., χ 1, it can be shown [56] that the following relations hold: Interested readers can find some approximate relationships for the variance of χ and R χ in [56, pp. 76-78]. For instance, the variance of the phase is better described with where κ χ is the concentration parameter of the Von Mises distribution. Equations (I.2)-(I.4) show that the precision of the mean direction and mean resultant length measurement critically depends on the number of elements constituting the directional data vector population χ .
APPENDIX II In this section, we show some experimental results obtained by applying the E-Mt-InSAR algorithm to conventional (nonphase-bias-corrected) ML interferograms. The method is detailed in [57]. Here, we want to point out that, for every analyzed SAR pixel, the optimized set of phases associated with the available N SAR acquisitions are exploited to reconstruct the set of M optimized, wholly time-consistent ML SB interferograms, namely, {φ ML opt,i } M−1 i=0 , which markedly have a zero phase closure. 4 Accordingly, once the E-Mt-InSAR noisefiltering method is applied, any subset of the optimized ML SB interferograms that could be extracted from and that involves all the available N SAR acquisitions is, 4 Note that the E-Mt-InSAR noise-filtering method also involves a subsequent nonlinear combination step between the original interferograms and those reconstructed from the optimized set of acquisition phase. It was proposed in [3, eq. (4)] to preserve the spatial coherence of the very coherent interferograms. This nonlinear operation, however, reintroduces some nonclosure phase triplets. by design, time-consistent. Therefore, the ground deformation time series that one could generate by inverting these subsets of optimized ML SB interferograms using any SBoriented Mt-InSAR method [18], [19], [20], [47], [54] are necessarily coincident. Of course, this statement is true when we exclude the effects of localized PhU errors and the role of some specific processing stages, such as those for the estimation of the residual topographic phase components and the compensation of the APS, whose results could depend on the selected subset of the optimized ML SB interferograms. Although the ground deformation time series obtained by inverting these different subsets of M optimized interferograms extracted from {φ ML opt,i } M−1 i=0 are coincident, however, this does not mean that the E-Mt-InSAR noise-filtering method [57] is capable of wholly compensating for the effects of the systematic phase biases in the ground displacement InSAR products, which could still be present even with coincident time series if the phase bias compensation was made incorrectly.
To prove the validity of these statements, we performed an experiment using the conventional SB ML interferograms related to the Nevada case-study area. We independently applied the E-Mt-InSAR noise-filtering method [57] to the set of SB interferograms selected with a maximum temporal baseline of 96 days [see Fig. 14(a)] and 12 days [see Fig. 14(b)]. The optimized noise-filtered interferograms were independently unwrapped [62] and inverted by the SBAS algorithm [18]. Fig. 14(c) shows the difference of the mean ground displacement velocity of the study area as obtained using the SB networks with maximum temporal baselines of 12 days [see Fig. 14(a)] and 96 days [see Fig. 14(b)].
The map only depicts SAR pixels with temporal coherence values greater than 0.95 to exclude that some of the observed differences could be ascribed to PhU errors.
The achieved results show that some differences in the ground displacement velocity maps persist. As a matter of fact, the E-Mt-InSAR method [57] was designed to mitigate the decorrelation noise in sequences of ML interferograms by exploiting their temporal relationships, but it did not consider specifically the presence of short-lived systematic phase contributions that, as said, occur at small/very-small temporal baselines. Nonetheless, depending on the used network of SB interferograms, the E-Mt-InSAR method allows obtaining time-consistent filtered ML interferograms that are characterized by perfectly compensated phase closures. We also remark that the method does not require that the SB interferograms form a planar triangulation in the temporal/perpendicular baseline plane. 5 To further prove the validity of the above-mentioned statements, we have carried out a second experiment. We have applied the E-Mt-InSAR algorithm to the network of SB interferograms shown in Fig. 14(a), characterized by a temporal baseline threshold of 96 days. Subsequently, from the generated group of 895 filtered interferograms, we have extracted a subset of them consisting of 125 filtered interferograms that form the SB network shown in Fig. 14(b) with a maximum temporal baseline of 12 days. The two sets of interferograms were unwrapped and the relevant ground displacement velocity maps were compared. Fig. 14(d) shows the mean ground deformation velocity difference between the results obtained with the SB networks at 96 days in Fig. 14(a) and 12 days in Fig. 14(b) but using for both the same set of optimized noise-filtered ML interferograms, as above clarified. The results show that the ground displacement velocity difference is almost zero everywhere but some localized zones with tiny PhU errors. Also, in this case, we only portray SAR pixels with temporal coherence values larger than 0.95. As said, this outcome does not automatically mean that the systematic phase contributions were perfectly compensated for. Conversely, it does mean that these phase artifacts cannot be easily discriminated by the (true) ground deformation signals because they are time-consistent and, hence, they do not lead to time discrepancies that could be straightforwardly visible by comparing the obtained ground displacement velocity maps.
We hope that the results of our analyses could be useful to unravel some of the (potential) reasons that are behind the apparent disagreement existing between the main outcomes of the investigations presented in [34] and [76].
Open Access funding provided by 'Consiglio Nazionale delle Ricerche-CARI-CARE' within the CRUI CARE Agreement Francesco Falabella (Graduate Student Member, IEEE) received the B.Sc. degree in computer science and information technology and the M.Sc. degree (cum laude) in information technology and telecommunications engineering from the University of Basilicata, Potenza, Italy, in 2016 and 2019, respectively, where he is currently pursuing the Ph.D. degree in engineering for innovation and sustainable development.
In 2019, he joined the Institute for Electromagnetic Sensing of the Environment (IREA), Naples, Italy, and the Institute of Methodologies for Environmental Analysis (IMAA), National Research Council (CNR) of Italy, Tito, Italy, as a Research Associate. His research interests include synthetic aperture radar (SAR) data processing, ground-based SAR interferometry, advanced multitemporal interferometric SAR (InSAR) techniques, high-performance computing (HPC), and environmental remote sensing and applications.
Mr Italy. In 2018, he was an Adjunct Professor of mobile telecommunication systems with the University of Naples Federico II. His main research interests include the development of advanced differential synthetic aperture radar interferometry (DInSAR) algorithms for the monitoring of surface deformation phenomena induced by subsidence, volcano activities, and earthquakes, with a particular interest toward the phase unwrapping problems. More recently, he has developed research activities for the generation of DInSAR products through multiplatform/multiangle and the new generation synthetic aperture radar (SAR) instruments, the generation of hybrid scan SAR to stripmap DInSAR analyses, the integration of SAR and optical images, and the analysis of land changes in flooded areas.
Dr. Pepe was a recipient of the 2014 Best Reviewer Mention of the IEEE GEOSCIENCE AND REMOTE SENSING LETTERS and the 2017 Best Reviewer Mention of Remote Sensing (MDPI) journal. He acts as a reviewer for several peer-reviewed international journals.