Conference Agenda

Differential Interferometric Synthetic Aperture Radar (DInSAR) is an advance remote sensing technique that is being widely used in studying and quantifying large-scale displacements due to anthropogenic and natural events, such as earthquakes, landslides, and volcano eruptions (Rosen 2000). More recently, various solutions allowing us to extend the original DInSAR technique, referred to as multi-temporal (MT) DInSAR approaches, have been developed to analyze the temporal evolution of the detected deformation phenomena through the generation of displacement time series. In particular, the basic rationale of the MT-DInSAR approaches is to generate an appropriate set of interferograms by pairing spatially and temporally distributed SAR acquisitions relevant to the same area and to invert such an interferometric stack to retrieve the displacement time series. Among the several MT-DInSAR techniques, Small Baseline Subset (SBAS) is a well-established and consolidated approach that has been widely used for the analysis of deformation events with millimetric accuracy (Berardino et al. 2002; Lanari et al. 2004).

On the other hand, in order to accurately retrieve ground deformation signals through the MT-DInSAR techniques, an effective and robust implementation of Phase Unwrapping (PhU) algorithm is necessary. In this context, several approaches have been proposed to address the PhU problem. Among them, Extended Minimum Cost Flow (EMCF) technique (Pepe and Lanari 2006) is categorized as one of the most effective and commonly used procedure within the SBAS processing chain. The rationale of this approach is to exploit the irrotational property in the temporal/perpendicular baseline plane , where the SAR acquisitions are coupled to generate the multi-temporal interferograms sequence, to spatially unwrap each interferogram of the dataset. More in details, the algorithm firstly estimates the unwrapped phase differences for each spatial arc (spatial gradients), that is created by connecting pixels in azimuth/range plane with the Delaunay algorithm, by using the MCF technique (Costantini and Rosen 1999) in temporal/perpendicular baseline plane. The estimated unwrapped phase differences are then exploited as starting point to perform 2D spatial PhU and solving the network in the azimuth/range plane via the standard MCF approach. The technique is effectively used in SBAS processing chain to carry out the phase unwrapping procedure.

In this work, we present a new PhU procedure developed to improve the EMCF performance by benefiting from the Compressive Sensing (CS) theory. We underline that the CS is an advanced signal processing technique that allows the robust and effective reconstruction of signals from under-sampled noisy measurements. The CS theory has been effectively applied in MT-DInSAR to estimate and mitigate the phase unwrapping errors (Manunta and Muhammad 2022). The presented PhU algorithm follows the same line of action as of the EMCF technique presented in (Pepe and Lanari 2006). Accordingly, the developed procedure consists of three main processing steps (see Fig. 1): (1) networks generation, (2) temporal PhU, and (3) spatial PhU. Each of the involved processing steps is briefly described in the following:

1. in the first step, the algorithm computes two Delaunay triangulation networks referred to as (a) temporal network in the temporal/perpendicular baseline plane, and (b) spatial network in the azimuth/range plane. Note that the temporal network is created by representing the SAR acquisitions in the temporal and spatial baseline plane and by using these acquisitions to create a Delaunay triangulation network, where each edge of the triangles represents an interferogram. Accordingly, this network identifies the sequence of DInSAR interferograms to be unwrapped. To generate the spatial network, we create a mask of coherent pixels, i.e., pixels having triangular coherence (Manunta et al. 2019) higher than certain threshold, that are common to all interferograms. From such a sparse grid of points, a spatial network is retrieved by applying the Delaunay triangulation algorithm; each arc of this network represents a spatial link between two pixels in azimuth/range plane and corresponds to a wrapped phase difference between two pixels in the azimuth/range plane. Further details about this step can be found in (Pepe and Lanari 2006);

2. the second step is most innovative part of the algorithm and is referred to as temporal phase unwrapping. It aims to estimate the unwrapped phase gradient of each arc in the azimuth/range plane by exploiting phase closure property in temporal/perpendicular baseline plane. In particular, this step computes the integer ambiguity vector that corresponds to interferograms to be corrected for the arc under consideration, by solving an under-determined system. Since we are considering arcs created by connecting pixels with very high value of triangular coherence (Manunta et al. 2019), it can be reasonably assumed that a physically based solution will be the one having minimum number of corrections, that can also be categorized as sparse vector. As we are in search of sparse solution vector, by following the CS theory (Manunta and Muhammad 2022) we solve the corresponding L₀-norm minimization problem by handling it as L₁-norm minimization problem, which is achieved by applying the properly modified Iterative Reweighted Least Square (IRLS) method as presented in (Manunta and Muhammad 2022). Since the algorithm does not have any constraint on the solution vector to be integer, the retrieved solution is rounded to nearest integer. To evaluate the quality of retrieved solution we compute temporal coherence γ_rnd to measure closeness of the retrieved solution to nearest integer vector. Note that γ_rnd varies from 0 to 1 where the value 1 implies the obtained solution is an integer vector. As the algorithm intends to estimate a sparse vector, we therefore set a threshold th_sparseto ensure the solution satisfies the sparsity constraints. Finally, similarly to EMCF approach, we develop a cost function to evaluate the confidence on the obtained solution. In particular, the cost function assigns very high values to the arcs having γ_rnd and th_sparse higher and smaller than fixed thresholds, respectively. Note also that high cost corresponds to more confidence on the solution and vice versa.

3. in the last step, spatial PhU of each interferogram is carried out by applying conventional MCF approach by taking into account the unwrapped phase gradient of each arc and assigned cost, retrieved in previous step, as an external information. By exploiting the computational efficiency of MCF approach, we perform multiple rounds of spatial phase unwrapping operations. Indeed, in each round we change the cost function by selecting different values for the threshold parameters th_sparse and γ_rnd.The final unwrapped interferometric stack is obtained by computing the weighted average of all the unwrapped solutions.

To evaluate the performance of developed PhU procedure we carry out an experimental analysis on a real SAR dataset acquired by Sentinel-1 descending orbits taken over the area related to Stromboli volcano (Italy), from May 2016 to April 2021. The data consists of 282 SAR acquisitions which are paired in 801 interferograms. Fig. 2(a) and (b) show the mean deformation velocity maps retrieved through the SBAS approach implemented with the developed PhU procedure and the conventional EMCF method, respectively. It is worth to note that in Fig. 2(a) and (b) only coherent points detected by each PhU procedure are shown in the maps. It is quite evident that the developed procedure outperforms in the area subject to strongest deformations, referred to as Sciara del Fuoco. To better explain this, we select a pixel located in Sciara dal Fuoco and present displacement time series retrieved by the developed algorithm (Fig. 2(c)), the conventional EMCF approach (Fig. 2(d)), and their difference (Fig. 2(e)). Note that the pixel has experienced a strong deformation of around 25 cm in 6 months, which was not detected by the EMCF procedure. This shows that the developed CS-based algorithm can significantly improve the quality of SBAS results, even if strong non-linear deformation signals occur.

We have developed a higher-level product data flow for routine reduction of InSAR time series that reduces the number of decorrelated radar pixels and thus yields spatially-comprehensive yet fine-scale observations of displacement. Our approach, combining persistent scattering analysis and small baseline analysis is not new in that each step has been amply developed in the literature [1-3], but including them in a data flow that can be used routinely to process a large amount of data reliably remains a challenge in a system with large throughput. Such methods are crucial for the next generation of InSAR satellites, such as NISAR or ROSE-L, but will be even more important for follow-on missions where data are acquired daily or hourly.

We demonstrate our system by reducing a set of Sentinel-1 data acquired over two areas: a portion of the San Andreas fault (SAF) in California and an aquifer system in Texas, not so much as to claim that we have produced a better map but mainly because the higher-level data product can be delivered routinely with little human intervention. Our results are shown in figs. 1 and 2, which display the average radar line of sight distance rate in radians per day. When the SAF data are converted to an assumed right lateral fault displacement the rate peaks at 3 cm/yr, similar to that observed using GPS [4]. The Texas deformation field is highly corrupted by decorrelation that proves a challenge for the phase unwrapping step, yet “blind” application of our algorithm yields a clear subsidence and uplift signal.

To yield these figures, we processed multiple years’ data acquisition (SAF: two years, 72 scenes, resulting in 670 interferograms, Texas: 5 years, 123 scenes, 662 interferograms) from Sentinel-1 A/B from level-0 raw data products, without any monitoring of intermediate results as will be necessary when the volumes of data produced daily by future systems exceed the capabilities of human oversight. We have found that our flow is usually quite reliable in the sense that the products are geophysically feasible, there is little adjustment of processing parameters on a site to site basis, and the computation is efficient. It thus can serve as a template for large-scale systems to deliver high level products that are readily usable by non-radar-specialist users.

Our data flow is as follows. First, we reduce level-0 measurements directly to geocoded single look complex (GSLC) images so that the data are in ingestible form. We stack these automatically coregistered products and compute an initial set of persistent scattering (PS) pixels using MLE methods. These candidate PS pixels are each tested using a cosine similarity criterion to remove the non-PS pixels from the candidates. We then test every remaining pixel against the new PS set to fill in as many non-PS pixels as possible. The next step is to interpolate the residual holes in every interferogram with a spiral interpolator to obtain complete coverage. Finally, the filled interferograms are processed as SBAS time series to obtain the displacement history of the region. Figs. 1 and 2 above are those series presented as an average line of sight rate.

The only parameter that we currently alter by hand is the maximum temporal temporal baseline in the SBAS analysis. As of this writing we do not have a reliable algorithm to automatically select this quantity. Choosing too small a value does not allow enough averaging to produce the mm-level accuracies we desire, and choosing too large a value results in aliasing that leads to underestimation of the displacement rates [5]. Nonetheless, our approach has shown itself reliable for many different types of terrains, ranging from the example shown here to volcanos or hydrological uplift and subsidence.

[1] Hooper, A., H. Zebker, P. Segall, and B. Kampes, A new method for measuring defor- mation on volcanoes and other natural terrains using InSAR persistent scatterers, Geophysical Research Letters, 31 (23), 5, doi:10.1029/2004GL021737, 2004.

[2] Berardino, P., G. Fornaro, R. Lanari, and E. Sansosti, A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms, IEEE Transactions on Geoscience and Remote Sensing, 40 (11), 2375 – 83, 2002.

[3] K. Wang and J. Chen, "Accurate Persistent Scatterer Identification Based on Phase Similarity of Radar Pixels," in IEEE Transactions on Geoscience and Remote Sensing, vol. 60, pp. 1-13, 2022, Art no. 5118513, doi: 10.1109/TGRS.2022.3210868.

[4] Tong, X., Sandwell, D. T., and Smith-Konter, B. (2013), High-resolution interseismic velocity data along the San Andreas Fault from GPS and InSAR, J. Geophys. Res. Solid Earth, 118, 369– 389, doi:10.1029/2012JB009442.

[5] Pepin, K., On the use of interferometric synthetic aperture radar for characterizing the response of reservoirs to fluid extraction and injection at wells, PhD thesis, Stanford University, Chapter 6, 2022.

Synthetic Aperture Radar (SAR) systems are powerful tools to monitor the climate change and to support adaptation actions. SAR systems provide important temporal information for monitoring changes over time. In the case of SAR systems, the first source of temporal information is the multitemporal series of the direct radar observables. For repeat pass configurations, and due to the interferometric capabilities of SAR systems, additional and complementary temporal information is captured by the complex interferometric coherence, also referred to as temporal coherence, constructed from pairs of SAR images [1], exploiting the temporal dependence between SAR images.

In classical interferometry, coherence is modelled as the product of several decorrelation terms [2], but these terms degrade the quality of the information that can be extracted from coherence. A group of these factors depends on the SAR system and the imaging process, such as the geometric decorrelation terms. Other decorrelation terms, such as the volume decorrelation [3], depend on the scatterer characteristics. For multitemporal SAR acquisitions, coherence depends also on the temporal decorrelation term, which accounts for changes in the scatterer between acquisitions. Therefore, scenarios characterized by a low temporal coherence, for instance, long temporal baseline configurations or datasets affected by strong weather effects, prevent the extraction of temporal information from coherence.

This work analyzes the extraction of common temporal information from multitemporal SAR and Polarimetric SAR (PolSAR) datasets in low temporal coherence scenarios, and establishes the relation with the temporal information extracted from direct radar observables. The basis of this study is on different previous works on coherence and change detection analysis. On the one hand, coherence has been proved to be helpful for land cover and vegetation mapping [4], where shorter temporal baselines, leading to larger coherence values, have been shown to perform better [1]. On the other hand, polarimetry has been employed for optimized change detection [5] and crop phenological monitoring based on the optimized polarimetric contrast ratio [6].

We propose a reinterpretation of the temporal coherence, linking the previous works under a common theory, demonstrating that coherence can be separated into two terms: a first symmetric term accounting for coherent changes and a second asymmetric term accounting for radiometric changes. For low temporal coherence scenarios, the symmetric term presents low values, preventing the use of coherence. Nevertheless, the information provided by the asymmetric term can be used in these circumstances to exploit common information between SAR acquisitions. We propose the use of this information, as an alternative to coherence, for information retrieval in low temporal coherence scenarios.

To prove the usefulness of this approach, we consider different classification strategies on Sentinel-1 (C-band), Radarsat-2 (C-band) and UAVSAR (L-band) where the improvements of the classification overall accuracy may range between 20% and 50%, compared to classification based on coherence. Our results demonstrate that the proposed approach can provide accurate information retrieval even in scenarios with low temporal coherence. This is important for monitoring climate change and supporting adaptation actions, as SAR systems can provide valuable insights into changes in the Earth's surface over time.

[1] A. W. Jacob, F. Vicente-Guijalba, C. Lopez-Martinez, J. M. Lopez-Sanchez, M. Litzinger, H. Kristen, A. Mestre-Quereda, D. Ziolkowski, M. Lavalle, C. Notarnicola, et al., “Sentinel-1 InSAR coherence for land cover mapping: A comparison of multiple feature-based classifiers,” IEEE J-STARS, vol. 13, pp. 535–552, 2020

[2] H. Zebker and J. Villasenor, “Decorrelation in interferometric radar echoes,” IEEE Trans. Geosci. Remote Sens., vol. 30, pp. 950–959, Sep 1992

[3] S. R. Cloude and K. P. Papathanassiou, “Polarimetric SAR interferometry,” IEEE Trans. Geosci. Remote Sens., vol. 36, pp. 1551–1565, Sept. 1998

[4] N. Joshi, E. T. Mitchard, N. Woo, J. Torres, J. Moll-Rocek, A. Ehammer, M. Collins, M. R. Jepsen, and R. Fensholt, “Mapping dynamics of deforestation and forest degradation in tropical forests using radar satellite data,” Environ. Res. Lett., vol. 10, no. 3, p. 034014, 2015

[5] A. Marino and I. Hajnsek, “A change detector based on an optimization with polarimetric SAR imagery,” IEEE Trans. Geosci. Remote. Sens., vol. 52, no. 8, pp. 4781–4798, 2013

[6] A. Alonso González, C. López Martínez, K. Papathanassiou, and I. Hajnsek, “Polarimetric SAR time series change analysis over agricultural areas,” IEEE Trans. Geosci. Remote. Sens., vol. 58, no. 10, pp. 7317– 7330, 2020

Interferometric synthetic aperture radar (InSAR) has shown great capability in detecting and measuring earth's surface deformation caused by different driving mechanisms. Despite this successful application, the efforts to describe the quality of InSAR deformation time series in terms of precision are limited in scope. Firstly, the exploitation of complete error propagation from raw SAR images to final deformation time series has received less attention due to the complexity of the processing steps, incomparable strategies and algorithms among different InSAR methodologies, and the large volume of spatio-temporal InSAR observations. Secondly, most of the previous studies on InSAR quality description and noise modeling have mainly focused on describing the characteristics of different noise components (e.g., scattering noise and atmospheric effects) in a single interferogram or subset of them, while the effect of time series processing and atmospheric filtering steps on the spatio-temporal noise variability in final time series is overlooked. Note that InSAR processing steps, mainly the spatio-temporal atmospheric filtering, affect InSAR deformation results, and consequently they alter the spatio-temporal noise structure of InSAR deformation estimates. Furthermore, as the filtering step and its setting varies from case to case, it is difficult to derive a generic stochastic model for InSAR deformation time series.

To cope with these limitations, the objective of this study is to develop a methodology to propagate all the error sources from raw InSAR observations to final deformation estimates in order to describe the quality and uncertainty of the final results. An efficient error propagation scheme through all the processing steps is designed. The main focus is on the second central moment of errors, or the error covariance matrix. Note that the primary output of TInSAR algorithms is a 3D spatio-temporal dataset, i.e., deformation time series for a large number (usually hundreds of thousands) of persistent or distributed targets. The final covariance matrix for all these observations is therefore a very large matrix, not possible to work with in practice. The goal is to derive an analytical closed-form expression to reconstruct the variances and covariances for any given spatio-temporal deformation estimate.

In this regard, we first exploit the existing body of knowledge about InSAR error sources to propagate errors through the interferogram generation steps and the time series phase retrieval (i.e., from statistics of SAR data to statistics of interferometric phase time series), and then we focus to a large extent on the further error propagation in the atmospheric effect mitigation step. For the latter, we formulate the filtering step in a mathematical framework based on the prediction theory (similar to least-squares collocation or Wiener filtering). This mathematical formulation is: i) generic enough to cover different existing methodologies and algorithms, ii) flexible enough to digest deterministic and stochastic assumptions about spatiotemporal behavior of different signal and noise components, and last but not least iii) simple enough to allow the application of the linear error propagation concept. Using this formulation, we present an analytical closed-form expression for the final covariance matrix of deformation time series. This covariance matrix not only comprises the information on the noise dispersion of individual time series, but also includes information about the correlation among noise components both in the space and in the time domain. The easy-to-use closed-form expression captures the majority of the assumptions and settings in a mathematical sense and can be applied to infer the contribution of different factors and processing settings. For instance, the proposed model comprises the initial noise structure of the data derived from the coherence matrix of targets and atmosphere noise models, applied multilooking factors, the number of acquisitions and time-interval between them, the type of deformation mechanism and nonlinear behavior of deformation, the temporal or spatial kernels used for atmospheric mitigation, and the distribution and location of targets used for spatial filtering/interpolation of atmospheric effects. It is shown that in InSAR deformation time series, we encounter with a highly correlated (or smooth) noise both in the space and in the time domain. The proposed model shows how much the characteristics of the smooth/correlated noise depend on the processing settings or on the initial error structure of input SAR data. The presence of such a spatio-temporally correlated noise can sometimes be misinterpreted as a real deformation. We demonstrate that how the proposed model helps to avoid such a misinterpretation. Finally, the validity of the proposed model and the derived closed-form expression is tested in different simulation scenarios, and also its application on a case study over a subsidence field in southern Tehran (Iran) based on Sentinel-1 data is demonstrated.