Blind Source Separation for MT-InSAR Analysis With Structural Health Monitoring Applications

Monitoring large areas of the Earth's surface is possible thanks to the availability of remote sensing data obtained by a large collection of diverse satellites orbiting the Earth. Specifically, multitemporal interferometric synthetic aperture radar (MT-InSAR) techniques supply the structural community with time series of line-of-sight displacements of terrain or structures, such as bridges or buildings, resulting from causes such as thermal expansion, contraction, and terrain deformation. The analysis of the different deformation signals observed is crucial in order to identify the different phenomena that cause the deformation. In this article, we explore the possibility of applying blind source separation algorithms to MT-InSAR, with the aim of developing methods towards automatic identification of different deformation patterns both in buildings and structures. We validate the proposed methodology using both synthetically generated datasets and real MT-InSAR data. We also provide a comparison with other similar methods. Our results demonstrate that InSAR time-series analysis can benefit from the use of the proposed blind source separation approach. Furthermore, the proposed technique is robust to the noise which is usually present in MT-InSAR data. This opens the door for monitoring infrastructure at scale over very large areas, helping to monitor civil infrastructures, and providing relevant insights to asset owners regarding the performance of such structures over time.

for agricultural industry, forest monitoring for environmental purposes, soil composition classification for mining exploration, and the study of natural hazards such as earthquakes and volcanoes. By combining remote sensing data together with machine learning techniques, these operations can be more efficiently performed over large areas [1]. The quantity and variety of the available remote sensing data has dramatically increased over the last decade due to the latest technological advances [2]- [4]. This opens the door to new applications specifically within the context of structural health monitoring as the resolutions and frequencies of accessible data increases [5].
Many different sensors can be used for earth observation, formed by employing different sections on the electromagnetic range [6]. For instance, optical sensors provide visually interpretable images by measuring backscattered light energy from the sun. They can be used for discriminating between different land-cover classes even at the chemical level [4], [7]. SAR (synthetic aperture radar) sensors form images which consist of amplitude and phase. These images can be collected irrespective of daylight or cloud, and in all weather conditions due to the active nature of these sensors. The SAR amplitude can be used for applications such as sea ice classification [8], whilst the phase component of SAR can be used in a technique called InSAR (interferometric synthetic aperture radar) [9]- [11]. In-SAR allows users to measure relative surface displacement in the line-of-sight (LOS) of the sensor with millimetre-scale accuracy, and even submillimetre accuracy under some conditions [12]. The consistent application of InSAR techniques at different acquisition times is known as multitemporal interferometric synthetic aperture radar (MT-InSAR). It produces a time-series of displacements. For instance, Fig. 1 depicts typical outputs from MT-InSAR processing. In this case, the time series were obtained after applying the RapidSAR methodology [13] from the commercial provider SatSense Ltd to Sentinel-1 data. Each dot on the map represents the average velocity in the satellite LOS, which can be derived where there is a coherent response of the SAR signal [14], [15].
A number of different SAR constellations are currently available for civil applications, with commercial and open-source options. For instance, TerraSAR-X is a high resolution commercial SAR sensor that allows us to periodically obtain SAR images every 11 days with a ground spatial resolution of approximately 3 m by 3 m down to 0.5 m by 0.5 m [3], [16]. Sentinel-1 images are provided free of charge by the European Commission's Copernicus programme every 6 [13] from the commercial provider SatSense Ltd. (b) Each persistent scatterer has a LOS measurement taken at the time of each acquisition, and the time series measurement plot is displayed for three selected persistent scatterers in different colours in the region of tunnel settlement marked in red in (a). 6 or 12 days in most of the rest of the world [17]. The spatial resolution of Sentinel-1 is lower, at approximately 20 m by 5 m. However, the overall image footprint covers a larger ground area. Despite a lower resolution, the regular and global service provided as part of the European Commission's Copernicus programme provides opportunities for continuous monitoring over wide-scale areas (each image covers a swath 250 km wide).

A. Related Works
MT-InSAR techniques supply the structural community with time series measurements of LOS displacements of terrain deformation [18] or structures such as bridges or buildings [19] whose deformation may result from terrain deformation or from causes such as thermal expansion contraction of the asset itself [20]. In this way, MT-InSAR techniques provide relevant insights to asset owners regarding the performance of such structures over time-they are potentially able to detect anomalous movements of the asset, or detect abnormal deformations caused by related factors (such as subsidence effects on buildings or scour effects on bridges). For example, in [21], MT-InSAR image data were used to study subsidence in Mexico City associated with groundwater extraction. In [22], MT-InSAR technology was employed to measure surface ground movements over time, in order to monitor all phases of a tunneling project. In [23], recently excavated subway tunnels were tracked due to very localized subsidence of the above surface along their path and several highways in the Shanghai urban area. In [24], urban assets were monitored in order to support local authorities to detect and evaluate subtle terrain displacements. MT-InSAR can also be used to monitor the behavior of single assets, with opportunities to combine these measurements with structural models and identify unusual behaviors [25]. In [26], the authors demonstrate how MT-InSAR time-series analysis can be used to early detect structural health problems in bridges, such as scour, before the bridge suffered a partial collapse. All these studies demonstrate that MT-InSAR images can be used to monitor infrastructure and obtain useful insights on valuable building and infrastructure assets. However, they require expert knowledge to process the data and interpret the results [19]. The processing is extremely technical and requires the selection of a number of different methods and a large number of parameters over different stages, based on each specific context.

B. Motivations
Expert analysis of MT-InSAR data can provide useful insights for a single asset. However, this approach prevents the application of such techniques at larger scales, e.g., to perform city-scale or larger analysis of structural behaviors and environmental influences, covering large areas and a large number of assets.
To be able to evaluate large areas at scale in an unsupervised way, machine learning algorithms that are able to identify different patterns of displacements over large datasets can be applied. This would open the door for the monitoring of infrastructure assets at scale over wide geographical areas [18], [27].
Machine learning approaches on this topic have been developed in the literature. For instance, in [28], they use independent component analysis (ICA) [29] and principal component analysis (PCA) [30] to discover spatial and temporal surface deformation patterns over San Francisco bay. The results obtained with high resolution TerraSAR-X data demonstrate the ability of these algorithms to identify different components of the motion measured across the region. Some were attributed to continuous landslide deformation motion, other components driven by precipitation-modulated pore pressure changes and some components captured more widespread seasonal deformation, such as groundwater level changes and thermal expansion of buildings.
A similar approach was used in [31], where they also used PCA, ICA, as well as nonnegative matrix factorization (NMF) [32] in order to isolate different latent sources of deformation of volcanoes using time series from MT-InSAR Sentinel-1 data. In [31], the authors introduced the concept of a linear mixture model. This model assumed: 1) the existence of a set of latent sources of deformation on the time domain; and 2) the deformation observed on the different points will be a linear combination of these latent source patterns. The quantities associated to each latent source pattern for each point are the mixtures. This model can be represented by the following expression: where Y ∈ R l×n is a matrix that contains the observed timeseries of the pont scaterers (PS) by columns, with l, representing the number of different acquisitions over time and n the number of different PS. M ∈ R l×p , M = [m 1 , . . . , m p ] is a matrix containing the latent source patterns m i , for i ∈ 1 . . . p by columns and A ∈ R p×n , A = [a 1 , . . . , a n ] contains by columns the mixtures a j , for j ∈ 1 . . . n associated to each latent source for each PS. Finally, N ∈ R l×n is a matrix which represents the noise introduced in the model by the acquisition process, unwrapping errors or other type of noise. In [31], the authors proposed using blind source separation algorithms such as ICA, PCA, or NMF in order to solve this problem and estimated M and A matrices from Y in an unsupervised way, allowing the study of large areas and many different datasets without human intervention.
In the MT-InSAR community, the use of blind source separation algorithms to solve this problem has had limited exploration. However, in the hyperspectral remote sensing community, this is a very well understood problem, known as unmixing. Specifically, linear spectral unmixing has the same definition to that represented in (1). Within the hyperspectral community the latent source patterns contained in M are typically know as the endmembers in the context of unmixing, while the mixtures contained in A matrix are known as the abundances. Due to the nature of signals encountered in hyperspectral imaging, two specific constraints are typically introduced into the model, which are: 1) the abundance nonnegativity constraint (ANC) [33]; and 2) the abundance sum-to-one constraint (ASC) [34], which enforce the abundances to be nonnegative (A i,j ≥ 0, i ∈ 1 . . . p, j ∈ 1 . . . n) and sum to one ( p i=1 A i,j = 1, j ∈ 1 . . . n) for each pixel, respectively. With these two definitions in mind, we can define the concept of purity of a given point. The point will be more pure the closer its maximum abundance is to 1 and the rest to 0, while the closer the abundances to 1/p the more heavily mixed. In the case that we do not apply the ASC, the concept of purity can also be defined as the difference between the maximum abundance and the rest, meaning that if the maximum abundance is similar to other abundances the point is heavily mixed. In the case of the ASC it does not make sense to apply it to MT-InSAR data, as the amplitude of the different time-series displacements is different and independent for each PS. For the case of ANC, it is more of an interpretability issue rather than physical matter as it is in hyperspectral imaging. If we consider negative abundances in MT-InSAR, then, a pattern with consistent positive displacements is going to be the same as another with consistent negative displacements, but with negative abundance. However, if we apply the ANC, then, a consistent positive displacement pattern will be a different from a negative consistent displacement pattern. Another example could be a seasonal periodic pattern in which we observe positive displacements in winter and negative in summer. This pattern would be the same as another periodic inverse pattern, i.e., positive displacements in summer and negative in winter. As positive and negative displacements may be caused by very different events it would make sense to apply the ANC to MT-InSAR data, so we should consider positive and negative patterns as different, rather than being the same pattern with different sign abundances.
The unmixing task is typically split into different stages depending on the approach [35]. The most common stages are: 1) estimation of the number of endmembers p and optional dimensionality reduction of the dataset; 2) estimation of the endmember matrix M; and 3) estimation of the abundances matrix A. These stages may depend on the method employed. There are methods that require a dimensionality reduction of the dataset, and there are other methods that perform stages 2) and 3) within the same stage. For first stage there are several methods proposed in the literature, with the most commonly used being virtual dimensionality [36] and hyperspectral subspace identification by minimum error [37]. For the endmember extraction stage the methods can be divided into two main categories [35]. The first are those that assume that the endmembers are present in the data, i.e., there exist some points in the dataset y k such that y k ≡ m i . The second category does not make this assumption. Finally for the third stage, there are different algorithms depending on the constraints that we consider: If we do not consider any constraint we have unconstrained least squares [34], if we consider the ANC we have nonnegative least squares [33], and if we also consider the ASC we have fully constrained linear spectral unmixing [34].
In this article, we explore the possibility of using traditional hyperspectral unmixing algorithms for MT-InSAR analysis with the aim of providing approaches for automatic identification of different structural behaviors (such as settlement, thermal dilation, heave, etc.). The overarching goal of this work is to automate the measurement, identification, and flagging of specific movement changes of structural assets as measured by regular InSAR time-series monitoring. We, therefore, aim to contribute towards solutions for automatic identification of potential structural problems on a city scale. For example, identifying unexpected subsidence could provide warning to potential damage to structures located in the zone of subsidence. Another example would be identifying a bridge oscillating with thermal dilation and contraction, whose seasonal behavior changes, which could alert asset owners to the fact that the bridge may have articulation or bearing problems [19].
We explore the possibility of discovering different patterns in the time series of displacements. In this way it may be possible to automatically discover abnormal behaviors or other potential structural issues to monitor in an urban environment. This information could be very useful for local authorities, or other private or public organizations for monitoring, planning, and forecasting maintenance activities [38].
The article is organized as follows. Section II describes traditional unmixing algorithms. Section III describes our proposed approach to tackle this problem, specifically adapted to the nature and characteristics of MT-InSAR data. Section IV describes the datasets used for our experiments, compares the traditional unmixing algorithms and our proposed approach in the task of discovering different temporal patterns in MT-InSAR data, and evaluates if traditional hyperspectral unmixing algorithms are suitable for MT-InSAR data analysis. Finally, Section V concludes this article.

II. METHODS
In this section we describe the methods used in this work. These are classic methods designed for linear spectral unmixing that we are going to apply to the task of MT-InSAR time-series analysis for discovering ground motion patterns. The methods considered are N-FINDR [39], orthogonal subspace projections (OSP) [40], vertex component analysis (VCA) [41], and simplex identification via split augmented Lagrangian (SISAL) [42].

A. N-FINDR
This algorithm searches for the set of points with the largest possible volume by inflating a simplex inside the data. The procedure begins with a random initial selection of points. It then iteratively calculates the resulting volume of replacing every point of the image in each endmember position. If the replacement results in an increase of volume, the point replaces the endmember. This procedure is repeated until there are no more endmember replacements. The mathematical definition of the volume of a simplex formed by a set of endmember candidates is proportional to the determinant of the set augmented by a row of ones. The determinant is only defined in the case where the number of features is p − 1 [43]. Since in hyperspectral data typically l p, a transformation that reduces the size of the input data is required. Often, the PCA transform has been used for this purpose. A possible shortcoming of this algorithm is that different random initializations of N-FINDR may produce different final solutions. In this work, we consider an N-FINDR algorithm implemented in iterative fashion, so that each sequential run was initialized with the previous algorithm solution, until the algorithm converges to a simplex volume that cannot be further maximized.

B. Orthogonal Subspace Projection (OSP)
This algorithm starts by selecting the point with maximum 2 norm in the scene as the first endmember. Then, it looks for the point with the maximum absolute projection in the space orthogonal to the space linearly spanned by the initial point, and labels that point as the second endmember. A third endmember is found by applying an orthogonal subspace projector to the original image [44], where the signature that has the maximum orthogonal projection in the space orthogonal to the space linearly spanned by the first two endmembers. This procedure is repeated until the desired number of endmembers, p, is found [40]. A shortcoming of this algorithm is its sensitivity to noise, since outliers are good candidates to be selected in the iterative process adopted by OSP. The VCA method discussed in the following section addresses this issue.

C. Vertex Component Analysis (VCA)
This algorithm also makes use of the concept of OSPs. However, as opposed to the OSP algorithm described above, the VCA exploits the fact that the endmembers are the vertices of a simplex, and that the affine transformation of a simplex is also a simplex [41]. As a result, VCA models the data using a positive cone, whose projection onto a properly chosen hyperplane is another simplex whose vertices are the final endmembers. After projecting the data onto the selected hyperplane, the VCA projects all image points to a random direction and uses the pixel with the largest projection as the first endmember. The other endmembers are identified in sequence by iteratively projecting the data onto a direction orthogonal to the subspace spanned by the endmembers already determined, using a procedure that is quite similar to that used by the OSP. The new endmember is then selected as the pixel corresponding to the extreme projection, and the procedure is repeated until a set of p endmembers is found [41].

D. Simplex Identification via Split Augmented Lagrangian (SISAL)
All of the previous described methods assume that the endmembers are present in the data; however, in a scenario where the different endmembers are highly mixed, this might not be the case. The SISAL algorithm does not rely on this assumption. In this case, the algorithms do not need the presence of pure pixels in the dataset in order to generate the endmembers. Most of the techniques in this category adopt a minimum volume strategy aimed at finding the endmember matrix M by minimizing the volume of the simplex defined by its columns and containing the endmembers. This is a nonconvex optimization problem which is much harder than those considered in the previous section in which the endmembers are assumed to belong to the input hyperspectral image.
Craig's seminal work [45] established the concepts regarding the algorithms of minimum volume type. Most of these algorithms formulate the endmember estimation as the NMF of the mixing and abundance matrices [46]- [50], with minimum volume constraint imposed on M. NMF is a hard nonconvex optimization, prone to get stuck within local minima. Aiming at obtaining lighter algorithms with more desirable convergence properties, [42] and [51]- [53] sidestep the matrix factorization by formulating the endmember estimation as an optimization problem with respect to Q = M −1 . The SISAL algorithm implements a robust version of the minimum volume concept. The robustness is introduced by allowing the ANC to be violated. These violations are weighted using a soft constraint given by the hinge loss function (hinge(x) = 0 if x ≥ 0 and −x if x < 0). After reducing the dimensionality of the input data from l to p − 1, SISAL aim at solving the following optimization problem: where Q ≡ M −1 , 1 p , and 1 n are column vectors of ones of sizes p and n, respectively, q m ≡ 1 T p Y −1 p with Y p being any set of linearly independent spectral vectors taken from the hyperspectral data set Y, and λ is a regularization parameter. Here, maximizing log(|det(Q)|) is equivalent to minimizing the volume of M.

III. PROPOSED FORMULATION
In MT-InSAR, the nature of the signals is different from hyperspectral imaging. This leads to substantial differences of the data itself. Two of the major differences between hyperspectral data and MT-InSAR data are: 1) the hyperspectral reflectance data domain is limited to 0-1, so the norm of the vectors is typically very similar or even constant; that is why in hyperspectral imaging it is very common to use the ASC. However, in MT-InSAR this is not the case. In MT-InSAR, the length of the displacement time series may vary a lot from one region of the image to another and from one time-frame to another. Also in MT-InSAR, the domain of the data is not limited to any value. This means that the Euclidean norm of the different points of the image may vary a lot and, therefore, the ASC it is not suitable for MT-InSAR data. 2) The signal to noise ratio (SNR) in MT-InSAR is considerably lower than in hyperspectral images. In hyperspectral images, the SNR is typically uniformly distributed within the same band and concentrated in some very well known bands due to well known physical effects. This is why it is common in hyperspectral imaging to directly remove the bands with low SNR, using only those bands with high SNR. However in MT-InSAR the noise is not equally distributed along all points or time measurements which leads to more difficulties in removing noise patterns and lower SNRs.
Taking into account these two characteristics of MT-InSAR data, we adopted a special case of Tikhonov regularization subject to constraint. We use the Euclidean norm for the misfit and 1 for the penalty factor, although other combinations of different norms for the misfit and penalty factor can be adapted in Tikhonov regularization depending on the application [54]. In our case we proposed the following formulation based on the sparse dictionary learning: where the λ||A|| 1 term aims at obtaining sparse abundances fractions with the λ parameter controlling the tradeoff between the data fidelity term and the sparsity term. That is, in highly noisy scenarios we can increase the importance of sparsity term on the abundances whilst relaxing the data fidelity term, which will help to reduce noise overfitting. Regarding the constraints, we included the nonnegativity constraint and also a normalization constraint on the endmember signatures. The aim of this is twofold: 1) by introducing this constraint on the endmember matrix it is possible to estimate the endmembers together with the abundances; 2) by forcing the norm of the endmembers to one, we aim to obtain the amplitude or size of the displacements in the abundance matrix A, so we can compare the amplitude of the displacements visually between different abundance maps, as the underlying endmembers patterns are normalized.
To implement the proposed formulation we employed the Python library "SParse Optimisation Research COde" (SPORCO) [55]. Specifically, we solved the proposed formulation via the basis pursuit denoising [56] with dictionary learning and nonnegativity constraint, which is implemented through an equivalent formulation using the alternating direction method of multipliers approach [57]. This implementation is very efficient, modular, and nevertheless provides very accurate results. It consists of an iterative and alternative optimization of each variable after splitting the variables for each term through an augmented Lagrangian procedure. In the following we will refer to this algorithm as blind sparse source separation (B-SSS).

IV. EXPERIMENTAL RESULTS
In this section, we perform an assessment of different spectral unmixing algorithms in the task of identifying different time-series patterns obtained from InSAR data. Our goal is to evaluate the applicability of these techniques for structural health monitoring purposes. This section is structured as follows. First, we describe how we produced several datasets created specifically for this work. After this we perform a quantitative and a qualitative assessment of the different algorithms analyzed using both the synthetic dataset generated and also a real dataset over the region of "The City" in central London.

A. Datasets Used
We generated synthetic data in order to being able to perform a quantitative evaluation of the different methods. We also consider a real dataset from Sentinel-1 located over "The City" in central London.
The procedure used to generate the synthetic images was based on the procedure used in [58] to generate synthetic hyperspectral images. However, in this case, instead of spectral signatures we use temporal patterns extracted from real MT-InSAR data and also perform a random sampling in the spatial domain to simulate the effects of noncoherence. Fig. 2 shows a block diagram of the procedure with the different steps numbered from 1 to 7. The procedures are separated in two columns, on the left hand column are the procedures that work on the spatial domain, on the right the procedures that work on the time domain, and on the middle the ones that work both on the spatial and space domains.
The first step is to generate a grey-scale fractal image with the diamond-square algorithm [59]. Several natural objects can be approximated by fractals to a certain degree, including clouds, mountain ranges, coastlines, etc., thus providing a baseline for simulating spatial patterns often found in nature. This image is then segmented using the k-means algorithm, such that each cluster corresponds to a different temporal pattern. A first version of the abundances matrix A is generated using this segmentation as baseline. The different rows for the matrix A correspond to the masks of the different classes, with ones if the pixel belongs to a particular class and zeros otherwise. The third procedure works with this initial abundances and tries to simulate natural mixtures by performing a Gaussian convolution on the spatial domain, such that, the points at the transition areas between two classes will result in more mixed (the maximum abundance value is lower) than the inner points of each segment, which will result in less mixed or more pure (the maximum abundance value is higher or close to 1). At the end of this procedure, we  obtain the final ground truth abundances A. Fig. 3(b) and (c) show examples of the output of this procedure.
Procedure number 4 selects a set of different temporal patterns corresponding to the M matrix that we are going to introduce in the synthetic dataset as endmembers. In order to select these patterns we ran the VCA endmember extraction algorithm over a real MT-InSAR dataset obtained with the Sentinel-1 sensor over a time-frame that ranges from 2015 to 2020, processed using the RapidSAR methodology [13]. Fig. 3(a) displays the different patterns considered as the endmembers (M) for this simulation. Once we have our A and M matrices, in step 5 we generate our ground truth data according to the model in (1). Finally in order to simulate noncoherence and environmental noise typically present in real scenarios we perform a random sampling on the spatial domain to simulate noncoherence in step 6 and a noise addition procedure to simulate environmental noise in step 7.
We created three different datasets with this procedure. The size of all datasets is 200 × 200 pixels. The number of patterns introduced in all cases is p = 16. Two different scenarios are generated depending on the kernel of the Gaussian convolution performed on step 3. For the first scenario generated with code-name "Syn-1," the kernel size and standard deviation was chosen to produce an scenario where the maximum abundance of the points is close to 1, or in other words, a scenario in which there are present the pure patterns on the dataset. For the second scenario with code-name "Syn-2," we manually tuned the kernel size and the standard deviation to produce a highly mixed scenario where the maximum abundance for the points varies from 0.56 to 0.91. For the noise addition procedure zero-mean Gaussian noise (N) was added to the synthetic scenes in different signal to noise ratios (SNRs)-from 5 dB to 10 dB-to simulate contributions from ambient and instrumental sources, following the procedure described in [37]. In the random sampling procedure to simulate the effects of noncoherence, a set of 10% of the the pixels were randomly selected from the original 200 × 200 pixel image resulting in a 40 000 points dataset.
We also generated a semisynthetic dataset using the Syn-2 image, but in this case, instead of Gaussian noise we added real MT-InSAR noise in step 7. In order to extract the noise from a real MT-InSAR data we applied the PCA procedure, retaining as noise the latest components with eigenvalues close to 0. In the following we will refer to this dataset as "Syn-real".
For illustrative purposes, Fig. 3(a) shows the endmember patterns used in the simulation. The abundance maps associated to each endmember pattern for both synthetic images are also displayed in Fig. 3(b) and (c). An example of an unwrapped interferogram corresponding to the date of 13th Jan 2017 in the case of the Syn-1 dataset with 5 dB of SNR is displayed in (d), while the corresponding randomly selected points are displayed in (e).
In addition to the synthetic images we also used real data from Sentinel-1 corresponding to the region of "The City" in central London. The acquisitions were collected from an ascending orbit at six-day intervals for a time-frame that ranges from 2015 to 2020. The time series were obtained after applying the Rapid-SAR methodology [13] from the commercial provider SatSense Ltd. RapidSAR is a time series algorithm developed by Spaans and Hooper [13]. The algorithm is designed to rapidly ingest new data and to handle pixels that have intermittent coherence. One key feature is that coherence is estimated from nearby pixels with similar scattering characteristics. In [60], the authors compared the results from RapidSAR with other processing algorithms. Fig. 4 shows some of the different displacement patterns that are observable in this real dataset. The PS marked in blue come from an area around Bank Underground metro station where there were tunneling construction activities, which result in very abrupt subsidence of the ground in that area. Panel 4(b) shows a noticeable downwards displacement starting in late 2017. The points in the two red regions show an overall trend of positive displacements toward the satellite. The reason is not known, but could be attributed to construction works. Panel 4(c) shows a light uptrend starting on late 2017. The green regions mark the Blackfriars road and rail bridges. Panel 4(d) shows a clear seasonal pattern of dilation and contraction. Finally, the purple region marks pier structures used by small boats and ferries on the river. Panel 4(e) shows that these structures also display a seasonal pattern, but, inverse to that present in the Blackfriars bridges.
Due to the lack of ground truth available for these datasets, we perform a qualitative assessment by comparing the output results of the different algorithms with regards to the reference patterns shown on Fig. 4(b)-(e).

B. Comparison of Different Methods for MT-InSAR Blind Source Separation
In this section, we will compare several algorithms in the task of blind source separation. In particular, we compare traditional algorithms used in the literature for this task, such as PCA or ICA proposed in [31] and [61], with other traditional unmixing algorithms used in hyperspectral imaging such as, VCA, N-FINDR, automatic target generation process (ATGP), or SISAL. We compare the results of these algorithms with the proposed B-SSS algorithm in Section III. The number of estimated endmembers extracted for each algorithm was p = 16, which corresponds to the ground truth number of endmembers. For all algorithms the parameters were optimized between a reasonable range of values for better performance. In the case of SISAL and B-SSS, several experiments with different values of λ were performed, choosing the λ value that provided the better result in each case.
To perform a quantitative assessment, we used the cosine similarity, defined in (4), between the ground truth endmember matrix used for generating the synthetic images M and the estimated endmember patterns for the different algorithms for each dataset M. This metric measured in degrees, measures the similarity in the shape of the different endmember patterns independently of its norm, and lower values indicate better results First, we will study the impact of the λ parameters on the results. Table I shows the averaged estimated cosine distances in the case of B-SSS for the different SNR ratios and different values of λ considered. In bold is marked the best result indicating the optimum value of the λ parameter in each case. The table shows that the higher the noise the higher the optimum value of the parameter, which means that the higher the noise the higher the impact of the 1 sparse regularization to improve the result and deal with the noise, as expected. Table II shows the averaged estimated cosine distances for the different considered algorithms in the case of the Syn-1 dataset for different SNR. The B-SSS formulation provides the best result in all the different levels of noise in this case. VCA, SISAL, and ATGP also provide good results, with values lower than 10 degrees in cases of higher SNR. The N-FINDR algorithm performs a little worse, while the PCA and ICA algorithms do a poor job of finding the different endmember patterns. This is because these algorithms are specialized in finding bases of the subspaces to better represent the data, considering thenegative coefficients also. The unmixing algorithms, on the other hand, are specialized in finding the more extreme pixels such that the abundances A are positive, which is how the patterns were generated. In other words, ICA and PCA are not looking for   Table III shows the averaged estimated cosine distances in degrees for the different considered algorithms in the case of the Syn-2 dataset for different SNR and for the "Syn-real" dataset. In this case the proposed B-SSS formulation also provides the best result in the cases of high noise: 5 dB, 10 dB, and real noise. This is because of the 1 norm sparsity term works well in preventing noise overfitting and generalizing the underlying patterns. The 1 sparsity term also works well in the case of real noise, in which case the noise is not uniformly distributed along time and space. This will be the case in real MT-InSAR images where the level of noise and its distribution is very different from hyperspectral imaging. This is the main reason why we propose this formulation. However, in the case of lower noise conditions and highly mixed patterns, the SISAL algorithm provides the best result. In this case, the algorithms: VCA, N-FINDR, and ATGP, which assume the presence of pure pixels in the image, do not provide a result as good as in the Syn-1 case. Also in this case, the PCA and ICA do not provide a good result for the same reasons as in the Syn-1 case.
For illustrative purposes, Fig. 5 shows an estimated endmember by the B-SSS algorithm (yellow dotted line), its corresponding ground truth endmember (purple dotted line) and the corresponding purest pixel with and without noise (blue dotted line and red line, respectively) for different noise scenarios. The example endmember in each case was chosen such that the error of the estimation was the closer to the average of all the endmembers for each case; in other words, the most representative example of the average error for all endmembers. The result provided by the B-SSS algorithm is very close to the ground truth endmember; even in the scenarios of high noise (5 dB), the result provided by B-SSS is almost noise-free, demonstrating the ability of the 1 term regularization to avoid noise overfitting and providing very accurate results. In the case of the Syn-2 dataset, the estimated endmember is closer to the ground truth than the purest pixel. This shows that although the proposed approach has not been able to estimate perfectly the ground truth endmember for this highly mixed scenario, the estimation is still purer than the purest pixel in the image. This demonstrates the ability of the proposed formulation to estimate endmembers in highly mixed conditions to some extent. Finally, Panel 5(f) shows a comparative of the estimated endmembers by the different algorithms VCA, N-FINDR, ATGP, SISAL, and B-SSS. The chosen endmember for this comparison corresponds to the first endmember, which has the highest norm and is the easiest endmember to identify by all algorithms. The figure shows the cumulative squared error between the chosen endmember and the estimated corresponding endmember for the different algorithms, so the lower the line, the better. The B-SSS provides the result with the lowest error.

C. Comparison of the Different Methods in the Real Scenario
In this section, we perform a qualitative assessment of the different algorithms in the case of the real MT-InSAR dataset, previously described over the region of "The City" in central London (Fig. 4). The number of endmembers p = 16 was manually selected after a qualitative analysis of the eigenvalues of the data correlation matrix. As for the synthetic data, all algorithm parameters were optimized through a qualitative analysis for all the algorithms and the endmembers that better matched the phenomena described in Fig. 4 were manually selected. Figs. 6-9 display the different abundances fractions and the endmember patterns associated to them corresponding to the different events displayed in Fig. 4. More specifically, the selected endmember patterns displayed correspond to the most similar patterns with regards to the average pattern shown on the Fig. 4(b)-(e), respectively.
Regarding the settlement due to tunneling at the Bank underground metro station, Fig. 6 shows the abundances fractions and the endmember pattern associated to this specific event for all the different algorithms. For PCA and ICA the pattern has been very poorly identified. For instance, in the case of PCA, the pattern related to this event is a generic positive pattern, which characterizes the event with high negative abundance. In the case of ICA, we observe a similar result to that provided by PCA but in this case the pattern provided is more noisy and less clear. However, in the other cases we observe that due to  defined than the previous cases with zero overall velocity from 2016 to mid 2017 and a more noticeable step in the second half of 2017. Also, on the abundances fractions maps, we observe that for B-SSS the pattern is more localized on the Bank area and the rest of the image has lower abundances than the rest of the algorithms. This is due to the sparsity constraint that we are enforcing, which promotes zero abundances on the pixels where the pattern is not present, removing these residual low abundances and concentrating the pattern on those regions where the pattern is present.   of the B-SSS algorithm, the result is again a sudden displacement upwards during 2018. This disparity of results shows that there is not a clear single source of these displacements in this area and that the displacements are actually a result of mixtures of different causes and events. This can explain why the different algorithms detect different patterns that are mixed in different proportions, and there is not a clear single pattern identified in that area. Fig. 8 shows the abundance fractions and the endmember pattern for the Blackfriars bridges. For PCA and ICA, all the periodic events are grouped in the same component, one with positive abundances and the other one with negative. This is not the case for the other algorithms, which, due to the ANC, identify two different endmembers, one for the Blackfriars bridges and another for the piers and platforms on the river. In this specific case, the algorithms ATGP, SISAL, and B-SSS provide the best results with a clearer temporal pattern around the bridges. For VCA and N-FINDR, the pattern is not so clear and is partially mixed or noisy.
Finally, Fig. 9 shows the abundance fractions and the endmember patterns for the periodic displacements associated to the pier structures on the river. In this case, the PCA and ICA results have been omitted as the result is the same as for the previous case, but with opposite sign on the abundance. For the other algorithms the pattern is more clearly identified by VCA, N-FINDR, and B-SSS, whilst it is partially mixed and noisy for ATGP and SISAL. Again, we observe that the pattern is more well defined and clearly localized in the case of the B-SSS algorithm. This is again due to the introduction of the sparsity constraint. Fig. 10 shows the reconstruction error in terms of normalized mean squared error (NMSE) between the original data Y and the estimated reconstructed data M A by the different algorithms. This metric is defined as follows: The figure shows that the algorithms ICA and N-FINDR provide the worse results in terms of reconstruction errors. The proposed algorithm B-SSS shows reconstruction errors in the same order as other methods such as PCA, VCA, or ATGP, while the best algorithm in terms of reconstruction error is SISAL. In general, the regions with higher reconstruction errors are those regions where there are not any significant displacements. On the other hand, the regions studied in this work show low reconstruction errors in general. This is expected, as the NMSE is an error measurement relative to the power of the signal. Thus, in the regions where the displacements are low, the relative power of the noise is higher than in those regions where there is a powerful displacement signal.

V. CONCLUSION
We have explored the use of traditional hyperspectral unmixing techniques applied to MT-InSAR time-series analysis in order to detect and characterize different temporal patterns, and their location in space and time. We also proposed a B-SSS algorithm, specifically designed for MT-InSAR data. The results show that traditional hyperspectral unmixing algorithms can be applied over MT-InSAR data in order to find patterns of particular interest over urban areas-such as subsidence-and different seasonal patterns in structures over urban areas. The traditional hyperspectral unmixing algorithms can even provide better insights than other traditional techniques typically used in MT-InSAR field, such as PCA or ICA. The proposed formulation generally provides better results in the case of MT-InSAR data, due to the fact that it is more robust to the high noise conditions often present in these type of data. The 1 sparse regularization term helps to deal with high noise conditions and provides better results in general, as it prevents over fitting. Using our B-SSS we demonstrate that it is possible to characterize different events occurring over large areas in an unsupervised way, both in terms of how they mix and also in the size of the displacements. This is an important step towards the automatic analysis of how different causes of displacements interact between themselves, and how they contribute to the performance of the structures. This approach may help stakeholders to better understand the performance of particular assets without the need to set up expensive equipment in situ. Additionally, it can be used to detect risks over large urban areas with little human intervention, and alert stakeholders where to carry out in situ monitoring at an early stage.
In the future, we will further explore the use of these techniques applied over very large datasets for identification of different sources of deformation, as well as methods for the automatic evaluation of the causes underlying the different deformation patterns. A semisupervised approach to identify these causes after the proposed analysis could help to perform this task, without human intervention helping to fully automate the processing and identification of possible problems and rising warnings when unexpected behaviors happen. Last but not least, in this article, we have focused on the application of the blind source separation method proposed on to infrastructures; however, the same methodology could be applied to discover latent deformation signals with other applications, such as, regional deformation analysis related to landslides or volcanoes. This is something that can be explored in the future.
VI. DATA AVAILABILITY Sentinel-1 data were obtained via the Copernicus Program, a joint initiative of the European Commission (EC) and the European Space Agency (ESA). RapidSAR data were processed by SatSense Ltd. TerraSAR-X satellite data were purchased by the Satellite Application Catapult through grant funding provided by Innovate U.K., supplied by Airbus for academic use.