On the Use of Weighted Least-Squares Approaches for Differential Interferometric SAR Analyses: The Weighted Adaptive Variable-lEngth (WAVE) Technique

This paper concentrates on the study of the Weighted Least-squares (WLS) approaches for the generation of ground displacement time-series through Differential Interferometric SAR (DInSAR) methods. Usually, within the DInSAR framework, the Weighted Least-squares (WLS) techniques have principally been applied for improving the performance of the phase unwrapping operations as well as for better conveying the inversion of sequences of unwrapped interferograms to generate ground displacement maps. In both cases, the identification of low-coherent areas, where the standard deviation of the phase is high, is requested. In this paper, a WLS method that extends the usability of the Multi-Temporal InSAR (MT-InSAR) Small Baseline Subset (SBAS) algorithm in regions with medium-to-low coherence is presented. In particular, the proposed method relies on the adaptive selection and exploitation, pixel-by-pixel, of the medium-to-high coherent interferograms, only, so as to discard the noisy phase measurements. The selected interferometric phase values are then inverted by solving a WLS optimization problem. Noteworthy, the adopted, pixel-dependent selection of the “good” interferograms to be inverted may lead the available SAR data to be grouped into several disjointed subsets, which are then connected, exploiting the Weighted Singular Value Decomposition (WSVD) method. However, in some critical noisy regions, it may also happen that discarding of the incoherent interferograms may lead to rejecting some SAR acquisitions from the generated ground displacement time-series, at the cost of the reduced temporal sampling of the data measurements. Thus, variable-length ground displacement time-series are generated. The mathematical framework of the developed technique, which is named Weighted Adaptive Variable-lEngth (WAVE), is detailed in the manuscript. The presented experiments have been carried out by applying the WAVE technique to a SAR dataset acquired by the COSMO-SkyMed (CSK) sensors over the Basilicata region, Southern Italy. A cross-comparison analysis between the conventional and the WAVE method has also been provided.


Introduction
Weighted Least-squares (WLS) is a generalization of ordinary least-squares (LS) in which non-negative weights are attached to data points to take into account the varying quality of the
The system of linear Equation (1) defines an LS optimization problem, which is solved as fully detailed in [14]. Precisely, for every identified pixel of the SAR image scene, the topographical error term ∆z(·) as well as a rough estimate of the deformation rate, namely x(·), are first estimated [14,72]. Subsequently, the unwrapped phases are corrected from the topographical errors, thus obtaining: r(P) sin ϑ(P) ·∆z(P) = 4π λ d k (P) + ε k (P), k = 0, 1, . . . , M − 1; ∀P ∈ Az × Rg (2) The system of linear Equation (2) can be synthetically expressed as [14]: Sensors 2020, 20, 1103 5 of 27 where A is the M × N incidence-like matrix of the linear transformation of Equations (2) and (3), Φ = 4π λ D is the vector of the (unknown) deformation values D(P) = [0, d 1 (P), . . . , d N (P)] at the given analysed SAR pixel (where it is implicitly assumed that d 0 (P) ≡ 0), and E is the vector of the error noise contributions. System (3) is solved in the LS sense. Nevertheless, in the more general case that the selected SB InSAR data pairs are arranged to form some independent subsets of data that are separated by large baselines [14], the matrix A turns out to be rank-deficient. In this case, the problem becomes: Problem (4) can be solved by decomposing the matrix A using the singular value decomposition method [67][68][69]. However, a solution that minimizes the norm of Φ is not physically feasible because it is responsible for significant discontinuities in the final deformation result. This problem was circumvented in [14] by manipulating the system of Equations (3) and replacing the phase unknowns Φ(P) = 4π λ D = [0, φ 1 (P), . . . , φ N (P)] with the phase velocities between time-adjacent acquisitions V(P) = [v 1 (P), . . . , v N (P)], which are defined as v k (·) = (φ k (·) − φ k−1 (·))/(t k − t k−1 ), k = 1, 2, . . . , N.
As a result, Problem (4) becomes: Φ : min ||B·V − ∆ Ψ || 2 min ||V || 2 (5) which is solved using the SVD method applied to the matrix B(M ≥ N) (see [14] for its mathematical expression). Mathematically, the matrix B is decomposed as follows: where Q is an (M × M) orthogonal matrix, Z is an orthogonal (N × N) matrix, and S is a diagonal (M × N) matrix, whose diagonal elements represent the singular values [σ 1 , σ 2 , . . . , σ N ] of the matrix B. Accordingly, the solution of Problem (5) is obtained as: with B + being the inverse generalized of B : B + = Z·S + ·Q T , and S + is the diagonal matrix S + = diag 1 σ 1 , 1 σ 2 , . . . , 1 σ N , 0, . . . , 0 ∈ R N×M . The phase values associated with every SAR acquisition are finally obtained by time-integrating the estimated values of the vector V. Once collected, the phases associated with the SAR acquisitions are converted to deformation measurements, and subsequently the Atmospheric Phase Screen (APS) is computed. Finally, the APS-filtered deformation time-series [10,14] are calculated.

EMCF-Based Multi-Temporal Phase Unwrapping
In this sub-section, we present an adaptation of the Extended Minimum Cost Flow (EMCF) PhU algorithm [60]. First, we give some introductory information on the EMCF phase unwrapping algorithm. The EMCF technique [60] represents an extension of the minimum cost flow PhU method [49,50]. It exploits two triangulations; the former is computed in the Temporal/Perpendicular baseline domain. Indeed, every SAR image can be represented in that domain with a point, and the arcs connecting different points represent the selected interferometric SAR data pairs. Those arcs are organized to form a (reduced) triangulation in the Temporal/Perpendicular baseline plane, see Figure 1, obtained by first computing a triangulation starting from the available SAR data and then discarding those arcs that involve at least one large baseline interferogram [60]. The latter triangulation is computed in the {Az × Rg} plane and encompasses the group of coherent SAR pixels that are common to the entire pool of the selected SB interferograms. The unwrapping operation is then performed. First, we identify all the arcs of the computed "spatial" triangles, and, for every arc, we independently carry out a "temporal" PhU Sensors 2020, 20, 1103 6 of 27 step, which is done by applying the MCF technique [49,50] to the temporal/perpendicular baseline grid. The second stage of the EMCF PhU procedure uses, for every arc, the previously computed unwrapped phases as the starting point for the successive spatial phase unwrapping operations. The latter are performed on every single interferogram by using the conventional MCF technique [49,50].
1, obtained by first computing a triangulation starting from the available SAR data and then discarding those arcs that involve at least one large baseline interferogram [60]. The latter triangulation is computed in the {Az × Rg} plane and encompasses the group of coherent SAR pixels that are common to the entire pool of the selected SB interferograms. The unwrapping operation is then performed. First, we identify all the arcs of the computed "spatial" triangles, and, for every arc, we independently carry out a "temporal" PhU step, which is done by applying the MCF technique [49,50] to the temporal/perpendicular baseline grid. The second stage of the EMCF PhU procedure uses, for every arc, the previously computed unwrapped phases as the starting point for the successive spatial phase unwrapping operations. The latter are performed on every single interferogram by using the conventional MCF technique [49,50]. Figure 1. Distribution of the SAR data into the temporal/perpendicular baseline plane for the COSMO-SkyMed dataset of the Basilicata region, used for the experiments shown in Section 4. In particular, the dataset is composed by 50 SAR images. Starting from these data, M = 263 small baseline InSAR pairs are identified. The dates are represented in the "day-month-year" format, and the label "CSK" stands for the satellite mission.
As initially explored in [73], here, we propose to adapt the EMCF procedure to enlarge the pool of the interferograms to be unwrapped and improve the overall PhU performance. To this aim, the selected interferometric SB SAR data pairs are grouped into two different sets, labelled as { 1 } and { 2 } , which are composed of 1 = ( ) and 2 = − 1 interferograms, respectively. The interferograms related to the first set identify the primary network, which forms a (reduced) triangulation in the temporal/perpendicular baseline plane (see Figure 2a) consisting of ( ) data pairs, see [58] for further details, also potentially selected using the sub-optimal procedure described in [22]. Figure 1. Distribution of the SAR data into the temporal/perpendicular baseline plane for the COSMO-SkyMed dataset of the Basilicata region, used for the experiments shown in Section 4. In particular, the dataset is composed by 50 SAR images. Starting from these data, M = 263 small baseline InSAR pairs are identified. The dates are represented in the "day-month-year" format, and the label "CSK" stands for the satellite mission.
As initially explored in [72], here, we propose to adapt the EMCF procedure to enlarge the pool of the interferograms to be unwrapped and improve the overall PhU performance. To this aim, the selected M interferometric SB SAR data pairs are grouped into two different sets, labelled as {S 1 } and {S 2 }, which are composed of M 1 = M (Tr) and M 2 = M − M 1 interferograms, respectively. The interferograms related to the first set identify the primary network, which forms a (reduced) triangulation in the temporal/perpendicular baseline plane (see Figure 2a) consisting of M (Tr) data pairs, see [58] for further details, also potentially selected using the sub-optimal procedure described in [22].
The primary network is ideally suited to the application of the original EMCF PhU T be the set of wrapped phases and be the corresponding set of the unwrapped phases, where P is a generic SAR pixel of the {Az × Rg} spatial grid. The PhU operations are performed over a set of highly-coherent SAR pixels that are in common to the whole set of interferograms of {S 1 }, namely G (Tr) ≡ {P ∈ S 1 }. These pixels can be identified either by computing the average spatial coherence of the interferograms or exploiting the phase triangulation closures (see the appendix of [73]).
( )] be the corresponding set of the unwrapped phases, where P is a generic SAR pixel of the {Az × Rg} spatial grid. The PhU operations are performed over a set of highly-coherent SAR pixels that are in common to the whole set of interferograms of { 1 }, namely ( ) ≡ { ∈ 1 }. These pixels can be identified either by computing the average spatial coherence of the interferograms or exploiting the phase triangulation closures (see the appendix of [74]). For the subsequent operations, the identification of the coherent SAR pixels that are specific for every SAR interferogram represents a very critical step. Over the interferograms of the primary network, as earlier said, a group of high-coherence SAR pixels is selected, and the unwrapping operations are done only on those SAR pixels. After, considering only the interferograms related to the primary network, a rough estimate of the LOS-projected displacement time-series, namely For the subsequent operations, the identification of the coherent SAR pixels that are specific for every SAR interferogram represents a very critical step. Over the interferograms of the primary network, as earlier said, a group of high-coherence SAR pixels is selected, and the unwrapping operations are done only on those SAR pixels. After, considering only the interferograms related to the primary network, a rough estimate of the LOS-projected displacement time-series, namelŷ d (1) N (P) , ∀P ∈ G (Tr) and the DEM errors ∆ẑ(P), ∀P ∈ G (Tr) are obtained by applying the (conventional) un-weighted SBAS procedure (see Section 2.1) or alternative methods. As a result of the removal of some triangles (i.e., those that involve at least one large-baseline interferogram) from the complete triangulation in the temporal/perpendicular baseline plane [60], it is possible that the reduced triangulation may involve N + 1 < N + 1 distinctive SAR images. However, this limitation of the original EMCF PhU technique is straightaway circumvented by complementing the interferograms of the primary network (see Figure 2a) with some additional interferograms (see Figure 2b).
These added interferograms are independently unwrapped following a strategy that is similar to the one initially adopted in [72]. In particular, for every interferometric SAR data pair, the phase unwrapping operation is performed over the relevant set of coherent SAR pixels (which are selected by imposing a threshold on the corresponding spatial coherence map) by using the conventional MCF procedure [50]. First, for every SAR pixel, we consider the local topography errors ∆ẑ(P), ∀P ∈ G (Tr) and the deformation time-seriesd (1) N (P) , ∀P ∈ G (Tr) (including possible atmospheric artefacts) obtained by processing the sequence of the primary network interferograms. The collected rough LOS-projected displacement time series, computed on the unwrapped interferograms of the first network, are temporally low-pass filtered and time-interpolated, and subsequently spatially interpolated to retrieve information on the incoherent pixels. The temporal and spatial interpolation steps are both performed using a simple polynomial interpolator. Thus, the estimated displacement time-series and the residual topographic components are used to generate a model of the (unknown) unwrapped phases. In particular, by considering the k-th interferogram of the group of the added interferograms {S 2 }, the phase model has the following expression: where ∆d (1) k is the LOS-projected relative deformation of the k-th InSAR data pair, ∆ẑ is the DEM error, and b (2) ⊥k is the perpendicular baseline of the k-th InSAR data pair. The modelled phase contributions are then modulo-2π subtracted from the phase of the interferograms, and the residual phase is then unwrapped by the conventional MCF technique [50].
Hence, the whole set of the unwrapped interferograms is obtained. We note that for every SAR interferogram of the primary network, the previously collected unwrapped phases are spatially interpolated over the specific set of coherent pixels of the given InSAR data pair, identified by considering those pixels with a spatial coherence value larger than 0.2.
The diagram flow of the whole process leading to the unwrapping of the complete set of the interferograms is shown in Figure 3.

Weighted Adaptive Variable-lEngth (WAVE) InSAR Method
In this sub-section, we will introduce the adaptation of the methods presented in Section 2.1 to the weighted case. The presented analysis is performed adaptively for every pixel P of the {Az × Rg} domain. As a starting point, let us consider the whole sequence of the unwrapped interferograms, say [∆ψ 0 (P), ∆ψ 1 (P), . . . , ∆ψ M−1 (P)] T , again. As earlier stated, for a given SAR pixel P, only a subset of the available M interferograms has an adequate spatial coherence value (i.e., larger than a prescribed threshold, for instanceγ = 0.2). Differently from [14] (and similarly to [70]), we propose to discard the interferograms that are less coherent (adaptively on a pixel-by-pixel basis). More precisely, for a given SAR pixel P, we propose to discard all the interferograms exhibiting a spatial coherence value less than γ = 0.2 and use the remaining ones for the generation of ground displacement time-series. Accordingly, for every single SAR pixel of the imaged scene, an exclusive network of DInSAR interferograms is selected and, then, inverted to generate the LOS-projected displacement time-series relevant to the selected SAR pixel. As an example, Figure 4 shows two different networks of InSAR data pairs related to two distinctive SAR pixels of the analysed area. More precisely, the DInSAR distributions in Figure 4a

Weighted Adaptive Variable-lEngth (WAVE) InSAR Method
In this sub-section, we will introduce the adaptation of the methods presented in Section 2.1 to the weighted case. The presented analysis is performed adaptively for every pixel of the {Az × Rg} domain. As a starting point, let us consider the whole sequence of the unwrapped interferograms, say [Δ 0 ( ), Δ 1 ( ), … , Δ −1 ( )] , again. As earlier stated, for a given SAR pixel , only a subset of the available interferograms has an adequate spatial coherence value (i.e., larger than a prescribed threshold, for instance ̂= 0.2). Differently from [14] (and similarly to [70]), we propose to discard the interferograms that are less coherent (adaptively on a pixel-by-pixel basis). More precisely, for a given SAR pixel P, we propose to discard all the interferograms exhibiting a spatial coherence value less than ̂= 0.2 and use the remaining ones for the generation of ground displacement time-series. Accordingly, for every single SAR pixel of the imaged scene, an exclusive network of DInSAR interferograms is selected and, then, inverted to generate the LOS-projected displacement time-series relevant to the selected SAR pixel. As an example, Figure 4 shows two different networks of InSAR data pairs related to two distinctive SAR pixels of the analysed area. More precisely, the DInSAR distributions in Figure 4a   Let us now describe the subsequent inversion step of the unwrapped interferograms. Our procedure has some similarities with [63,69] but, in our case, the selection and removal of some interferometric SAR data pairs may lead, as in [14], to an undetermined problem and (potentially) to the generation of variable-length ground displacement time-series.
Before going in detail, let us first spend a few additional moments on the setting of the weights Let us now describe the subsequent inversion step of the unwrapped interferograms. Our procedure has some similarities with [63,69] but, in our case, the selection and removal of some interferometric SAR data pairs may lead, as in [14], to an undetermined problem and (potentially) to the generation of variable-length ground displacement time-series.
Before going in detail, let us first spend a few additional moments on the setting of the weights associated with the phase measurements, in correspondence with the SAR pixels that have a coherence value larger than the considered threshold (e.g., in our experiments,γ = 0.2). To this aim, the variance of the multi-looked interferometric phases is, first, estimated pixel-by-pixel. It is well known that the density probability function (pdf) of the wrapped multi-looked phase φ (ML) is [7,74,75]: where L is the number of the looks, φ is the average multi-looked (true) phase, Γ(·) is the gamma function, and F(·) is the hypergeometric function. Accordingly, the phase variance can be numerically calculated as: However, for the low coherence regions, the estimated coherence values typically result in being biased. Moreover, the phase variance calculation is somewhat computationally expensive, requiring the numerical integration in Equation (10). Alternatively, we can exploit the basic principles of the directional statistics [76]. Indeed, under some simplified assumptions [76], it can be shown that a suitable approximation for the pdf of the wrapped phase (i.e., a directional data vector) is the wrapped normal distribution WN(ϑ; µ, ρ), which is obtained by wrapping onto the circle the normal distribution N µ, σ 2 , that is: where ρ = exp − σ 2 2 is the resultant mean length of the circular phase data [76]. In our case, if we refer to the k-th SAR differential interferogram and the generic (multi-looked) SAR pixel P, the circular data represents the unitary phase vector φ k,0 , φ k,1 , . . . , φ k,L−1 . It is relevant to the single-look phases that averaged together (see Figure 5), one gets the k-th multi-looked phase at the SAR pixel P, where L is the number of averaged phase samples. Accordingly, [76]: where j = √ −1 is the imaginary unit, and: Sensors 2020, 20, 1103

of 27
Sensors 2020, 20, x FOR PEER REVIEW 12 of 28 Figure 5. Sketch of the pixel grid of the k-th interferogram exploited for spatial multi-look operations.
In the red box are highlighted the single-look phases used for the computation of the variance relative to the generic multi-looked pixel P of radar coordinates (az,rg).
Hence, the phase variances can be easily computed by Equation (13), and the DInSAR interferograms are generated without requiring highly demanding computations as with Equation (10). Alternatively, it is possible to consider the Cramer-Rao bound for the phase variance [6], which is expressed as follows: Once the phase variances are estimated, a suitable way to set the weights that are associated with every single InSAR measurement (i.e., the multi-looked interferometric phase related to every InSAR data pair) is to consider the reciprocal of the phase variances [29], namely Using the matrix formalism, we can also write that: Equation (15) is fundamental for the analysis of the error propagation in the weighted case. To introduce the problem, we refer to the fundamental principles of the uncertainty theory applied to the linear regression problems [79]. We note that, as earlier described, the un-weighted SBAS procedure leads to the solution of the optimization Problem (4), which is done by Equation (7). For a given SAR pixel, the known term ΔΨ, representing the vector of the unwrapped phases, is a random phase vector, and ΔΨ is its related covariance matrix: where the elements on the principal diagonal of the covariance matrix represent the phase variances of the DInSAR interferograms used for the SBAS inversion, and the off-diagonal element of ΔΨ , say , 2 ≠ , represents the covariance between the (unwrapped) phases of the i-th and the j-th InSAR data pair. By applying the basic principles of the linear regression and the uncertainty theory [79] to Equation (7), the covariance matrix of the unknown velocity vector is derived as follows: In the red box are highlighted the single-look phases used for the computation of the variance relative to the generic multi-looked pixel P of radar coordinates (az,rg).
Hence, the phase variances can be easily computed by Equation (13), and the DInSAR interferograms are generated without requiring highly demanding computations as with Equation (10). Alternatively, it is possible to consider the Cramer-Rao bound for the phase variance [6], which is expressed as follows: Once the phase variances are estimated, a suitable way to set the weights that are associated with every single InSAR measurement (i.e., the multi-looked interferometric phase related to every InSAR data pair) is to consider the reciprocal of the phase variances [29], namely W = [ω 0 , ω 1 , . . . , ω M−1 ] with: Using the matrix formalism, we can also write that: Equation (15) is fundamental for the analysis of the error propagation in the weighted case. To introduce the problem, we refer to the fundamental principles of the uncertainty theory applied to the linear regression problems [77]. We note that, as earlier described, the un-weighted SBAS procedure leads to the solution of the optimization Problem (4), which is done by Equation (7). For a given SAR pixel, the known term ∆ Ψ, representing the vector of the unwrapped phases, is a random phase vector, and C ∆ Ψ is its related covariance matrix: where the elements on the principal diagonal of the covariance matrix represent the phase variances of the M DInSAR interferograms used for the SBAS inversion, and the off-diagonal element of C ∆ Ψ , say σ 2 i, j i j, represents the covariance between the (unwrapped) phases of the i-th and the j-th InSAR data pair. By applying the basic principles of the linear regression and the uncertainty theory [77] to Equation (7), the covariance matrix of the unknown velocity vector V is derived as follows: In our treatment, although the computed interferograms have a given grade of dependence (because obtained by the same sequence of SAR data), we assume that the interferometric phases are uncorrelated. This simplified assumption is widely adopted in the literature, and only very few InSAR studies have addressed the problem of the systematic calculation of the correlation among a group of InSAR interferograms computed from the same set of SAR data [78]. Under this simplified hypothesis, Equation (17) becomes: Let us now describe the developed Weighted Adaptive Variable-lEngth (WAVE) DInSAR technique, which permits us to select and adaptively process a set of differential SAR interferograms. Let us refer, again, to a generic SAR pixel P of the multi-looked grid, and let [∆ψ 0 , ∆ψ 1 , . . . , ∆ψ M−1 ] T be the corresponding vector of the unwrapped phases, [γ 0 , γ 1 , . . . , γ M−1 ] T be the relevant spatial coherences, and W = [ω 0 , ω 1 , . . . , ω M−1 ] be the weights. As earlier said, for every SAR pixel, a preliminary selection of the "good" interferograms is made up. Hence, only the group of coherent interferograms survives. Therefore, let I = i 0 , i 1 , . . . , iM −1 be the indexes of the coherent interferograms, withM being the number of selected ones. The extracted vectors of the unwrapped phases and the relevant weights

respectively.
Similarly to [14], because of the adaptive selection of the interferograms to be used and the discarding of the remaining ones, the SAR data can generally be organized in disjointed subsets, which have to be properly linked. In our approach, we also check the temporal overlapping of the produced data subsets. From our analyses, we discard those SAR pixels for which there is no time overlap between the data subsets. Note also that the extracted interferograms now do not necessarily involve all the available SAR images. Indeed, only a group of SAR acquisitions are preserved, let φ 0 ,φ 1 , . . . ,φN T be the unknown phase relevant to thoseN + 1 "survived" SAR images, which are those acquired at ordered times t 0 ,t 1 , . . . ,tN T . As for the un-weighted SBAS case, the problem is still formulated concerning the vector of the unknown velocities between contiguous SAR imagesV = v 1 , v 2 , . . . , vN The weighted inversion problem to be solved is: where the matrix B is also adaptively obtained for the given SAR pixel taking into account the distribution of the "survived" SAR acquisitions and the relevant InSAR data pairs, following the same lines of [14]. The norm indicated in Equation (19) is defined as in [79], namely ||X || W = X T ·W·X 1 2 .
Accordingly, Problem (19) leads to the minimization of the following norm: subjected to a minimum-norm constraint on the velocity. Note that, in our case, W = diag ω 0 , ω 1 , . . . , ωM −1 , see Equation (15), and the vector e represents the phase residuals between the unwrapped phases and those reconstructed from the retrieved solution, namely e = B·V − ∆Ψ.

of 27
The Problem (20) represents the particularization of the following more general weighted least-squares problem: where the matrixes W RM ×M and T RN ×N put weights on the row and column of the matrix B RM ×N , respectively; its solution is unique. It is obtained, see [79], by first calculating the W-T-singular decomposition of matrix B [79], assuming thatM ≥N, which decomposes B such that U −1 ·B·Z = diag µ 1 , µ 2 , . . . , µN , where µ i , i = 1, . . . ,N are the singular values of the decomposition, U is orthogonal to W, and Z is orthogonal to T.
Note that, the W-T decomposition of matrix B is obtained [79] as follows. First, the matrixes of the constraints W and T are written by the Cholesky factorization as: W = L·L T , T = K·K T . Then, the SVD decomposition of the matrix C = L T ·B·K −T is obtained, such as Q T ·C·X = diag µ 1 , µ 2 , . . . , µN . Finally, the decomposition of the matrix B is done by setting U = L −T ·Q and Z = K −T ·X.
The final solution is eventually obtained as follows: where: Z ∈ RN ×N , D + B = [1/µ 1 , 1/µ 2 , . . . , 1/µ r , 0, 0, .., 0] ∈ RN ×M with r ≤N, and U ∈ RM ×M . In our case, matrix T is the identify matrix of RN ×N and, by Equation (5), the weights are set . Accordingly, the matrixes W and T are diagonal and: As a consequence, we have: Of course, the solution (23) is now pixel-dependent, because both the selected InSAR data pairs and the associated weights are different from pixel to pixel. In critical regions, it may also happen that the pre-selection of the "good" interferograms, for the given SAR pixel, leads some SAR acquisitions to be discarded. In this latter case, the problem in (19)-(21) is mathematically defined in the same way as in the previous sub-section. Still, the vectors of the unknown velocities have a reduced length (i.e.,N is less than N). As a consequence, once the estimated velocity values are time-integrated to recover the corresponding LOS-projected deformation time-series, namelyφ = φ 0 ,φ 1 , . . . ,φN −1 , the latter have a reduced length (i.e., less than N + 1).
The diagram flow chart of the WAVE method is shown in Figure 6. We would like to remark that although the proposed WAVE method has been developed in the framework of the SB approaches, its main characteristics allow a simple extension to the more general non-SB context.
Finally, let us spend a few moments on the propagation of errors for the weighted case. By considering Equation (22), it is easy to prove that: whereB + = Ω·W 1/2 , and C ∆Ψ is the covariance matrix of the selected unwrapped phases. Under the hypothesis that the interferograms are uncorrelated, namely C ∆Ψ = W −1 , Equation (24) simplifies as: We would like to remark that although the proposed WAVE method has been developed in the framework of the SB approaches, its main characteristics allow a simple extension to the more general non-SB context.
Finally, let us spend a few moments on the propagation of errors for the weighted case. By considering Equation (22), it is easy to prove that: where ̂+ = Ω • 1 2 ⁄ , and ΔΨ is the covariance matrix of the selected unwrapped phases. Under the hypothesis that the interferograms are uncorrelated, namely ΔΨ = −1 , Equation (24) simplifies as: To check the quality of the obtained variable-length LOS displacement time-series, a weighted version of the temporal coherence factor [60], which also takes into account the weights [ 0 , 1 , … ,̂− 1 ] of the selected interferograms, has also been introduced: where ̂( ) is the number of above-coherence-threshold interferograms at the given pixel , and and are the indexes of the master and slave time acquisitions of the k-th InSAR data pair. , the co-variance matrix of the velocity terms is constant for every element of the (estimated) velocity vector, thus implying that errors propagate in the computed deformation time-series uniformly over time.
To check the quality of the obtained variable-length LOS displacement time-series, a weighted version of the temporal coherence factor [60], which also takes into account the weights ω i 0 , ω i 1 , . . . , ω iM −1 T of the selected interferograms, has also been introduced: whereM(Q) is the number of above-coherence-threshold interferograms at the given pixel Q, and IM k and IS k are the indexes of the master and slave time acquisitions of the k-th InSAR data pair.

Dataset and Study Area
The study area where we have successfully applied the WAVE methodology encompasses the Basilicata region, Southern Italy (see Figure 7). The territory of Basilicata is predominantly mountainous. It is surrounded by the Mediterranean Sea, with a relatively small flatland zone near the coastal area overlooking the Ionian Sea, which is one of the two outlets to the sea, where the other ones are on the Tyrrhenian side. Moreover, the region is crossed by the Apennine Mountains that extend from the Ligurian Alps to the Northern Sicily mounts, with famous massifs, such as Mount Sirino (2005 m), Mount Pollino (2267 m), Mount Volturino (1835 m), and many other ones. The hills are almost entirely composed of clay soil, which is subjected to erosion processes caused by the frequent rainfalls, mostly in the winter season. These processes are the starting trigger for landslides, which have afflicted and still afflict several cities [80][81][82][83][84], such as Montemurro, Latronico, Pomarico, Maratea, Craco, and Montescaglioso. It is also important to underline the high seismic risk of the entire area, as confirmed by the numerous earthquakes [85][86][87] that have occurred during the last 200 years, which, in many cases, have destroyed several small cities [88], where thousands of people tragically died. The intertwining of these issues, including also the hydrogeological instability, contributes to the overall instability of zone. That results in a problematic scenario for correctly applying the conventional interferometric SAR techniques. The available SAR dataset is composed of 50 SAR acquisitions, collected by the COSMO-SkyMed (CSK) sensors of Italian Space Agency (ASI), from which a set of 263 differential SAR interferograms were generated.
Sirino (2005 m), Mount Pollino (2267 m), Mount Volturino (1835 m), and many other ones. The hills are almost entirely composed of clay soil, which is subjected to erosion processes caused by the frequent rainfalls, mostly in the winter season. These processes are the starting trigger for landslides, which have afflicted and still afflict several cities [82][83][84][85][86], such as Montemurro, Latronico, Pomarico, Maratea, Craco, and Montescaglioso. It is also important to underline the high seismic risk of the entire area, as confirmed by the numerous earthquakes [87][88][89] that have occurred during the last 200 years, which, in many cases, have destroyed several small cities [90], where thousands of people tragically died. The intertwining of these issues, including also the hydrogeological instability, contributes to the overall instability of zone. That results in a problematic scenario for correctly applying the conventional interferometric SAR techniques. The available SAR dataset is composed of 50 SAR acquisitions, collected by the COSMO-SkyMed (CSK) sensors of Italian Space Agency (ASI), from which a set of 263 differential SAR interferograms were generated.  For our purpose, we have used HH polarized StripMap (HIMAGE) images, collected from February 2012 to November 2018 along ascending orbit, each with a right-looking radar direction and viewing angles of approximately 34 • . Every raw SAR acquisition has been first focused [90,91] to obtain the focused In-Phase and Quadrature terms for each acquisition, e.g., the Single Look Complex (SLC) images. The information of the entire SAR dataset is summarized in Table 1, where for every SAR image, we list: (i) the acquisition dates, (ii) the (average) Doppler Centroid (DC) and (iii) the orthogonal baseline (calculated w.r.t the common, reference master image acquired on 1 October 2014). Every SAR image covers a small area of the Basilicata region, between the meridians ranging from 15.6 • to 16.0 • and the parallels ranging from 40.0 • to 40.4 • . In this zone, the Pertusillo dam, the Lagonegrese, and the Val D'Agri territory, as well as the Vallo di Diano territory, are included. In Figure 8, an optical image of the area is shown, and some locations of particular interest are highlighted. orthogonal baseline (calculated w.r.t the common, reference master image acquired on 1 October 2014). Every SAR image covers a small area of the Basilicata region, between the meridians ranging from 15.6° to 16.0° and the parallels ranging from 40.0° to 40.4°. In this zone, the Pertusillo dam, the Lagonegrese, and the Val D'Agri territory, as well as the Vallo di Diano territory, are included. In Figure 8, an optical image of the area is shown, and some locations of particular interest are highlighted.

Experimental Results
In this section, we show the results achieved by applying the WAVE method to the set of SAR data presented in Section 3. The experiments were conducted by first selecting from the available N = 50 SAR acquisitions a group of InSAR data pairs, selected by imposing a maximum perpendicular separation of 800 m and a maximum time span of 2 years, respectively. As a result, we have identified M = 263 InSAR data pairs. Once chosen, the M differential SAR interferograms have been generated. To this purpose, precise satellite orbital information and a three-arc Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) of the area have been used to remove the topographic phase components. The computed differential SAR interferograms have also been independently multi-looked (10 azimuth and 10 range looks, respectively) and pre-filtered by using the approach in [92]. As explained in Section 2, the interferometric SAR data pairs are organized into two independent groups. The former corresponds to the primary interferometric data pairs network shown in Figure 2a, and the latter represents the set of additional interferograms (see Figure 2b). These two sets of interferograms are then unwrapped by using the strategy outlined in Section 2.2.
To investigate the potential of the WAVE approach, we first processed the sequence of the M unwrapped interferogram using the conventional, un-weighted SBAS approach. Figure 9 shows the mean displacement velocity map of the investigated area, overlaid to an amplitude SAR image of the zone. Only the SAR pixels with a temporal coherence greater than 0.6 have been depicted on the map. We also show in Figure 9b-e the LOS-projected displacement time-series of some selected SAR pixels, which are located in areas subjected to sensible deformation signals. The shown displacement time-series are obtained after the estimation and removal of the Atmospheric Phase Screen (APS). The latter operation is carried out by the cascade of a High-Pass temporal filter and a Low-Pass spatial filter [10] applied to the achieved LOS-projected, un-filtered deformation time-series.
zone. Only the SAR pixels with a temporal coherence greater than 0.6 have been depicted on the map. We also show in Figure 9b-e the LOS-projected displacement time-series of some selected SAR pixels, which are located in areas subjected to sensible deformation signals. The shown displacement timeseries are obtained after the estimation and removal of the Atmospheric Phase Screen (APS). The latter operation is carried out by the cascade of a High-Pass temporal filter and a Low-Pass spatial filter [10] applied to the achieved LOS-projected, un-filtered deformation time-series. Subsequently, we have applied the WAVE method to the same set of unwrapped interferograms. We remark that, even though we pre-selected a set of SB interferograms, all the possible DInSAR interferograms could be potentially processed. However, this would turn out to be extremely time-consuming. Similarly to [14], however, the removal of the low-coherent interferograms may lead to the generation of multiple disjointed subsets, which are linked by using the weighted version of the SVD strategy, as detailed in Section 2. Furthermore, in some noisy areas, the possibility of having LOS-projected displacement time-series with variable length also arises. The application of the WAVE technique to the whole set of SAR pixels has led to the generation of the corresponding LOS-projected displacement time-series. To check the quality of the achieved results, we have calculated the weighted temporal coherence; see Equation (26), namely {Γ( )} ∈{ × } . With respect to the temporal coherence factor initially introduced in [60], in our case, we only consider the set of DInSAR interferograms that has actually been used during the inversion step for the generation of the surface displacement time-series. Accordingly, the weighted temporal coherence factor turns Subsequently, we have applied the WAVE method to the same set of unwrapped interferograms. We remark that, even though we pre-selected a set of SB interferograms, all the possible DInSAR interferograms could be potentially processed. However, this would turn out to be extremely time-consuming. Similarly to [14], however, the removal of the low-coherent interferograms may lead to the generation of multiple disjointed subsets, which are linked by using the weighted version of the SVD strategy, as detailed in Section 2. Furthermore, in some noisy areas, the possibility of having LOS-projected displacement time-series with variable length also arises. The application of the WAVE technique to the whole set of SAR pixels has led to the generation of the corresponding LOS-projected displacement time-series. To check the quality of the achieved results, we have calculated the weighted temporal coherence; see Equation (26), namely Γ(P) P∈{Az×Rg} . With respect to the temporal coherence factor initially introduced in [60], in our case, we only consider the set of DInSAR interferograms that has actually been used during the inversion step for the generation of the surface displacement time-series. Accordingly, the weighted temporal coherence factor turns out to be more sensitive to the unwrapping errors in the medium-to-high coherent regions. Conversely, the weighted temporal coherence does not take into account the decorrelation artefacts that were present in the interferograms that were adaptively discarded, because of the low values of their spatial coherence. Additional constraints have also been imposed to discard some extra SAR pixels with low quality. To do that, we have generated three spatial maps containing, for every SAR pixel, the information related to: (i) the number of interferograms actually used during the inversion of the unwrapped phases, namely N I (P) P∈{Az×Rg} ; (ii) the number of time acquisitions of the computed variable-length displacement time-series, namely N D (P) P∈{Az×Rg} ; and (iii) the number of independent data subsets that are linked by the WSVD decomposition described in Section 2.3. Figure 10 shows the map of the number of used interferograms, as obtained in our case. Finally, the set of the well-processed SAR pixels are those that satisfy the following multiple conditions: whereγ, ι and δ are the thresholds on the weighted temporal coherence, the number of used interferograms, and the number of used data, respectively.
that were present in the interferograms that were adaptively discarded, because of the low values of their spatial coherence. Additional constraints have also been imposed to discard some extra SAR pixels with low quality. To do that, we have generated three spatial maps containing, for every SAR pixel, the information related to: (i) the number of interferograms actually used during the inversion of the unwrapped phases, namely { ( )} ∈{ × } ; (ii) the number of time acquisitions of the computed variable-length displacement time-series, namely { ( )} ∈{ × } ; and (iii) the number of independent data subsets that are linked by the WSVD decomposition described in Section 2.3. Figure  10 shows the map of the number of used interferograms, as obtained in our case. Finally, the set of the well-processed SAR pixels are those that satisfy the following multiple conditions: where γ , ι and δ are the thresholds on the weighted temporal coherence, the number of used interferograms, and the number of used data, respectively. Figure 10. A spatial map of the investigated area in SAR coordinates related to the number of interferograms used for each pixel by the WSVD time-series inversion. Figure 11 shows the map of the mean displacement rate of the area as obtained by using the WAVE approach, where only the well-processed SAR pixels satisfying the multiple Conditions (27) have been portrayed. As a result, we got a total number of 696,898 well-detected SAR pixels against the 132,399 corresponding values obtained by processing the same set of unwrapped interferograms using the un-weighting SBAS procedure, with an improvement of roughly 526%. We also show in Figure 11b-e four examples of the obtained LOS-projected displacement time-series. In particular, Figure 11b,c show the plots related to two SAR pixels (see their location in the map of Figure 11a) characterized by one single data subset and with full-length displacement time-series (i.e., consisting of 50 time samples), and Figure 11d,e show the plots related to two SAR pixels with variable-length deformation time-series. We also show in Figure 12a-d the relevant networks of InSAR Figure 10. A spatial map of the investigated area in SAR coordinates related to the number of interferograms used for each pixel by the WSVD time-series inversion. Figure 11 shows the map of the mean displacement rate of the area as obtained by using the WAVE approach, where only the well-processed SAR pixels satisfying the multiple Conditions (27) have been portrayed. As a result, we got a total number of 696,898 well-detected SAR pixels against the 132,399 corresponding values obtained by processing the same set of unwrapped interferograms using the un-weighting SBAS procedure, with an improvement of roughly 526%. We also show in Figure 11b-e four examples of the obtained LOS-projected displacement time-series. In particular, Figure 11b,c show the plots related to two SAR pixels (see their location in the map of Figure 11a) characterized by one single data subset and with full-length displacement time-series (i.e., consisting of 50 time samples), and Figure 11d,e show the plots related to two SAR pixels with variable-length deformation time-series. We also show in Figure 12a-d the relevant networks of InSAR interferograms, as represented in the temporal/perpendicular baseline plane, and related to the SAR pixels whose LOS displacement time-series are shown in Figure 11b interferograms, as represented in the temporal/perpendicular baseline plane, and related to the SAR pixels whose LOS displacement time-series are shown in Figure 11b-e. Finally, the LOS-projected displacement time-series have been compared to available external GPS measurements. GPS data were downloaded from the website [95] and projected along the sensor line-of-sight. Figure 13 shows the results of this cross-comparison analysis in correspondence to the seven GPS stations available in the area. The red stars in Figure 13 highlight the location of the GPS sites in Figure. Table 2 shows the values of the root mean square error (RMSE) of the difference between the obtained LOS-projected displacement time-series and the corresponding GPS measurements.  Finally, the LOS-projected displacement time-series have been compared to available external GPS measurements. GPS data were downloaded from the website [93] and projected along the sensor line-of-sight. Figure 13 shows the results of this cross-comparison analysis in correspondence to the seven GPS stations available in the area. The red stars in Figure 13 highlight the location of the GPS sites in Figure. Table 2 shows the values of the root mean square error (RMSE) of the difference between the obtained LOS-projected displacement time-series and the corresponding GPS measurements.

Figure 12. (a)-(d)
InSAR distribution in the temporal/perpendicular baseline plane related to the SAR pixels labelled as I, II, III, and IV in Figure 11.

Discussion
We now discuss the main topics that emerged in this work, taking into account the attained experimental results. We have developed a hybrid InSAR method that complements the adaptive selection of a suitable set of interferograms with a weighted least squares inversion strategy, to improve our capability to monitor the evolution of surface displacement, especially in the low-coherent regions, by considerably increasing the number of measurements points. This is done, in correspondence with some SAR pixels, at partial expenses of the time sampling of the measurements (i.e., the number of sampling SAR acquisition times). Therefore, we have produced the so-called variable-length LOS-projected displacement time-series. The benefit of this approach with respect to the conventional, un-weighted SBAS procedure is that it preserves all the principal features of the SBAS technique (and potentially of other alternative SB approaches) but complements it with a spatial adaptive selection procedure, which allows working always with moderate-to-high coherent interferograms, thus mitigating the impact of noise decorrelation and phase unwrapping artefacts. Even though the WAVE method represents an extension of the SBAS procedure, it does not impose, in principle, any pre-selection of the small baseline InSAR data pairs to be used. For reducing the computational burden of the algorithm, not all possible InSAR interferograms are generated, but only a fraction of them. Because of this, the WAVE approach does not necessarily belong to the class of the SB methods. Even in the large baseline interferograms, some isolated pixels still preserve their coherence, and they can be profitably used for the generation of ground displacement time-series. Therefore, our approach avoids the time-consuming inspection (and in some cases, the potential manual removal) of the InSAR pairs to be used for the generation of the deformation time-series. The use of the multiple-criterion (27) for the identification of the well-processed SAR pixels also allows quantifying the improvement in terms of analysable SAR pixels. The weights used during the interferograms inversion procedure are essential to reduce the impact of noise artefacts in the obtained DInSAR products and, as demonstrated by Equation (25), they are set in such a way to have LOS-projected displacement time-series with uniformly time-distributed uncertainty values. The WAVE method can be applied to any set of multi-looked SAR interferograms, no matter the InSAR tools or the pre-processing (e.g., the noise-filtering) methods that have been used to generate and process the interferograms. In addition, the WAVE method can be straightforward implemented in concurrent and/or parallel computing architectures, thus only requiring independent, pixel-by-pixel analyses. Further efforts are required to study how to set these weights in the most general case that the interferograms are correlated (as they are). Also, the WAVE method has to be validated in wider geophysical contexts, especially in areas with moderate-to-low coherence, where the un-weighted SBAS approach (and other alternative SB methods) still suffers from significant decorrelation noise and phase unwrapping disturbances. All these issues are matters for future investigations.

Conclusions
This study has addressed the problem of generating variable-length LOS-projected displacement time-series in moderate-to-low coherent areas using a hybrid InSAR method, which is named as Weighted Adaptive Variable-lEngth (WAVE) technique. The paper provides an introduction of the weighted least-squares (WLS) InSAR methods, with a specific focus on those applied in the small baseline InSAR framework. The mathematical and statistical fundamentals of the WAVE method have been detailed in the work. The pros of the WAVE method are essentially related to the possibility to extend the number of detectable pixels, but at the potential costs of the measurement time sampling. However, this represents a reasonable compromise, especially in moderate-to-high-vegetated areas, such as for the Basilicata region, where the decorrelation noise and the phase unwrapping mistakes may lead to unsatisfactory results. The experimental results have shown the validity of the proposed method in the selected study area. The new WAVE method requires further extensive studies to demonstrate its general applicability for the retrieval of deformation signals in various geophysical contexts. As earlier stressed, we would like to remark that the WAVE technique naturally fits with the requirements of the SBAS technique as well.