Conference Agenda

This text has been modifed to remove all equations and figures. Please see attached pdf for full abstract.

1. Introduction

Land subsidence in the Netherlands is becoming an increasingly critical issue as it is closely linked with sea level rise, flooding risks and greenhouse gas emissions due to peat oxidation [1,2]. Despite the importance of this issue, it is very difficult to accurately assess subsidence levels across the country. Radar Interferometry (InSAR) is a very promising technique for monitoring land surface motion at large spatial scales with frequent temporal sampling. While InSAR techniques employing stable point scatterers (PS) have been successfully used to monitor subsidence in the Netherlands [3,4,5], these PS points are usually founded at greater depths and the movement of the surrounding landscape has had to be indirectly inferred.

Attempts to directly monitor the peatland surface with distributed scatterer (DS) techniques have encountered significant challenges. One such challenge is the seasonal loss of interferometric coherence every spring, which results in a discontinuous phase time series. Figure 1 illustrates the problem of seasonal coherence loss. Sufficiently coherent interferometric combinations can be made between epochs in the autumn and winter seasons (indicated by the dashed red boxes), allowing time series analysis to be carried out. However, for several months every spring and summer, the observed interferometric coherence is so low that no useful information is likely to be present in any interferogram made with an acquisition during this period.

A further complication is the fact that there are no coherent combinations which can be made between the two individual coherent periods (the NE and SW regions of the matrix), which implies that the two coherent periods are disconnected. We denote this phenomenon with the term loss-of-lock. The disconnect between the two coherent periods means that the gap between them is no longer constrained by integer ambiguities; there exists an unknown real-valued shift between the periods which must be resolved in order to obtain a single, consistent time series. This shift represents the unknown displacement history which cannot be measured that occurred during the incoherent periods.

2. Methodology

2.1 Contextual Data

We assemble a database of combined public cadastral parcel delineations and land cover data, soil maps, and groundwater management zones, available: [6]. These factors play a critical role in either the movement of the land surface, the scattering properties which affect the radar observation, or both. By cross-referencing this data with the SAR imagery, we can assign each pixel to a known parcel ID with known soil, land use, land cover. This ensures that we are processing homogeneous observations which are representative of the same land surface movement phenomena. We multilook the SLC Sentinel-1 observations according to the parcel delineations of the contextual dataset. This is a natural division to make, as the land cover, soil type, and groundwater are approximately consistent within a parcel. The "EMI" method [7] is used to estimate a consistent set of phases in equivalent single master form.

2.2 Segment Identification and Phase Unwrapping

Coherent time series segments are identified by stipulating a minimum number of consecutive epochs in which the daisy-chain coherence exceeds a given threshold. These segments are subsequently unwrapped using the methodology described in [8]. Typical values used to identify a segment are a minimum of five consecutive epochs with gamma > 0.1. At this stage, we are left with a number of unwrapped time series segments, with an unknown displacement between each segment.

2.3 Group Displacement Model Estimation

We postulate that neighbouring parcels with matching land use, land cover, soil type, and groundwater management can be expected to behave in a similar manner, such that we can bridge the incoherent data gap described in Section 1 by combining the coherent observations of several similarly behaving regions to estimate a single set of common displacement model parameters. This model can then be used to estimate the vertical shifts between the time series segments. While the model parameters in X and the shifts Delta z can theoretically be estimated simultaneously, the high degree of correlation between these unknowns can result in a very poor estimation. Instead, we note that the shift is common for all phases in a given segment. Thus by taking the difference in time between phases, the shift term drops out and the model parameters of X can be estimated directly. The shift for a given coherent segment can subsequently be estimated by taking the average difference between the model and the phase time series over the coherent period T.

The selection of the model is an important consideration and can be accomplished by multiple hypothesis testing, which is planned for a future publication. In this abstract, we show the results of an empirical hydrological model based on precipitation and evapotranspiration. Values for these model inputs are provided as daily mean values by the Royal Dutch Meteorological Institute (KNMI).

3. Results and Discussion

The methodology is tested in an area of interest around Zegveld, NL. This area is chosen due to the large peat deposits in the area, and the availability of in-situ validation data. Validation data is provided by extensometer measurements which provide a continuous time series of soil height measurements at one location [9]. The root mean squared error (RMSE) is evaluated between the group median result for the period of overlap (May 2020 - Jan. 2022), giving an RMSE of 6.7 mm. It should be noted that we do not expect these two measurements to match exactly, as the InSAR result is the median of a large spatial extent, while the extensometer data is from a single point.

4. Conclusion

We demonstrate a new methodology for estimating the ground motion of cultivated peatlands using DS time series InSAR. We show how discontinuities in a decorrelated time series can be bridged by considering the measurements of nearby similarly behaving regions. Our initial results show that the approach is promising, and we have been able to successfully validate our result against the ground truth data we have available with a low degree of error. To our knowledge, this is first accurate multi-year InSAR measurement of peatland surface motion in the Netherlands.

Acknowledgement

This research is part of the Living on Soft Soils (LOSS): Subsidence and Society project, and is supported by the Dutch Research Council (NWO-NWA-ORC), grant no.: NWA.1160.18.259, URL: nwa-loss.nl.

References

[1] G. Erkens, M. J. van der Meulen, and H. Middelkoop, “Double trouble: Subsidence and CO2 respiration due to 1,000 years of Dutch coastal peatlands cultivation,” Hydrogeology Journal, vol. 24, no. 3, pp. 551–568, 2016.

[2] G. Erkens, T. Bucx, Dam, R. D. Lange, and J. G. Lambert, Sinking Cities: An Integrated Approach to Solutions, In: The Making of a Riskier Future: How Our Decisions Are Shaping Future Disaster Risk. World Bank, 2016.

[3] M. Caro Cuenca and R. F. Hanssen, “Subsidence due to peat decomposition in the Netherlands, kinematic observations from radar interferometry,” in Proc. ESA Fringe Workshop, (Frascati, Italy), pp. 1–6, 2008.

[4] M. Caro Cuenca, R. F. Hanssen, A. Hooper, and M. Arikan, “Surface deformation of the whole Netherlands after PSI analysis,” in Proc. ESA Fringe Workshop, (Frascati, Italy), pp. 19–23, 2011.

[5] R. F. Hanssen, F. J. van Leijen, G. Erkens, E. Stouthamer, K. Cohen, and Others, “Land motion service of the Netherlands.” https://bodemdalingskaart.nl/en-us/, 2018.

[6] “Publieke Dienstverlening Op de Kaart (PDOK).” www.pdok.nl.

[7] H. Ansari, F. De Zan, and R. Bamler, “Efficient phase estimation for interferogram stacks,” IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 7, pp. 4109–4125, 2018.

[8] P. Conroy, S. A. N. Van Diepen, S. Van Asselen, G. Erkens, F. J. Van Leijen, and R. F. Hanssen, “Probabilistic estimation of InSAR displacement phase guided by contextual information and artificial intelligence,” IEEE Transactions on Geoscience and Remote Sensing, vol. 60, pp. 1–11, 2022.

[9] S. van Asselen, G. Erkens, and F. de Graaf, “Monitoring shallow subsidence in cultivated peatlands,” Proceedings of the International Association of Hydrological Sciences, vol. 382, pp. 189–194, 2020.

Multitemporal SAR interferometry (MT-InSAR) is one of the most exploited phase-based InSAR techniques, capable of achieving millimeter per year accuracy [3,4,5,7] depending on the number of acquisitions and the spatial scale of the processing at which the displacement rate is measured. It helps to overcome the unwanted effects that may overwhelm the displacement patterns in standard InSAR, in particular the atmospheric delays that often dominate the individual interferograms [9]. Several techniques have been developed to handle the large stacks of SAR data. Some of these techniques, the distributed scatterers and the small baseline methods, employ spatial averaging to reduce the signal-to-noise ratio and to extend the spatial coverage of deformation measurements beyond the persistent scatterers present only in urbanized areas. This averaging in all cases involves a mixture of pixels that are more or less affected by changes of soil moisture and vegetation. It is a complex averaging process and therefore non-linear.

Recently the scientific community realized that, when applying spatial averaging, there was an additional phase delay that biases the deformation recovery in natural areas, which cannot be related to other terms already managed such as the topographic or the atmospheric effects [1]. The amount of bias is small in individual interferograms but its accumulation in time can significantly affect the final estimated velocities. The causes behind this effect remain unclear and debatable [2,6,10]. Researches related this bias to the non-zero closure of multi-looked interferometric phase on specific land covers (vegetation and croplands). Due to the complex spatial averaging, unbiased areas, as roads, will also be affected by the bias, which limits their potential to analyze displacements at mm/yr accuracy.

As part of the broad objective of bias mitigation, we were interested in this study in deepening our understanding of the phase bias in order to be able to limit its effect on unbiased non-natural areas. To achieve our goals, we used a long Sentinel-1 track covering the France territory from South to North, initially processed using the automatic FLATSIM service [8]. The velocity map obtained by automatic processing is, as expected, highly biased in areas where the land cover is dominated by croplands and forests. The bias was observed inversely correlated with the number of unwrapped one year interferograms per pixel. In order to mitigate the bias, we processed again the interferograms, starting with the products provided by FLATSIM service, that is, wrapped interferograms multilooked by a factor 8 in range and 2 in azimuth (hereafter called 2-looks interferograms).

In a first step, we started by analyzing the bias to build a good proxy for biased or unbiased pixels in 2-looks radar geometry. That for, we constructed a time series using all the shortest baseline unfiltered 8-looks interferograms. The resulting velocity map was compared with the high-resolution THEIA land cover map. Then, averaged phase time-series were computed for each land cover allowing the understanding of bias accumulation and evolution through time. These displacement time-series confirmed that the urban areas are stable over time. Rice is found to be the cropland with the highest bias while vineyards only suffer from moderate bias. We observed that for the croplands, the bias is mainly accumulated during the period of plants growth and stabilizes during the harvest period. This common behavior of almost all the croplands indicates that the observed bias might be related to physical properties of plants during the growth season (size, humidity, etc.). However, the complete loss of coherence on vegetated pixels during harvest prevent any return to zero. This contrasts with the seasonal behavior of forests characterized with cyclic seasonal motion, where we also note a different behavior between broad-leaved and coniferous forests. In parallel we also compared available information in 2-looks, such as the temporal coherence, amplitude dispersion, and interferogram amplitude variability, to the land cover. This allowed the extraction of useful statistical properties of each land cover to distinguish its pixels based on SAR data only.

In a second step, we based on these results to propose a proxy for biased pixels in 2-looks, and a methodology to unmix the reliable unbiased pixels from the biased ones. Therefore, a map of unmixing coefficients was built providing a confidence indicator for each pixel based on their statistical properties. This map, which gives high values to stable unbiased pixels and low ones to biased pixels, will replace the amplitude of 2-looks wrapped interferograms to be used as a weight when multilooking into 8-looks, prior to filtering, unwrapping and time series inversion. New time series are then computed again using the “unmixed” smallest baseline interferograms. The resulted velocity map is still affected by strong bias in crop areas, however, roads, isolated farms, etc. are now devoid of bias. The use of weighted moving average filter with the calculated unmixing coefficient as a weight, allowed us to keep track of the isolated unbiased pixels that are likely to be mixed and hidden in the surrounding bias. We performed a statistical comparison of velocities before and after unmixing, as a function of the type of land cover and as a function of a multi-looked version of the unmixing coefficient. We show that the proxy used for bias is relevant for isolating bias-prone pixels. Such kind of methodology is important, for practical reasons, for services computing massive numbers of interferograms such as ARIA, FLATSIM or LICS, which only provide multi-looked interferograms.

[1] H. Ansari, F. De Zan, and A. Parizzi. Study of systematic bias in measuring surface deformation with sar interferometry. IEEE Transactions on Geoscience and Remote Sensing, 59(2):1285–1301, 2021.

[2] F. De Zan, M. Zonno and P. López-Dekker, "Phase Inconsistencies and Multiple Scattering in SAR Interferometry," in IEEE Transactions on Geoscience and Remote Sensing, vol. 53, no. 12, pp. 6608-6616, Dec 2015.

[3] M.P. Doin, F. Lodge, S. Guillaso, R. Jolivet, C. Lasserre, G. Ducret, R. Grandin, E. Pathier, and V. Pinel. Presentation of the small baseline nsbas processing chain on a case example: the etna deformation monitoring from 2003 to 2010 using envisat data. Proceedings of the ESA Fringe 2011 Workshop, Frascati, Italy, (19-23 September 2011), 2011:19–23, 2011.

[4] A. Ferretti, C. Prati, and F. Rocca. Permanent scatterers in sar interferometry. IEEE Transactions on Geoscience and Remote Sensing, 39(1):8–20, 2001.

[5] A. Gabriel, R. Goldstein, and H. Zebker. Mapping small elevation changes over large areas: Differential radar interferometry. Journal of Geophysical Research: Solid Earth, 94(B7):9183–9191, 1989.

[6] Y.E. Molan, Z. Lu, and J.W. Kim. Influence of the statistical properties of phase and intensity on closure phase. IEEE Transactions on Geoscience and Remote Sensing, 58 (10):7346–7354, 2020.

[7] P. Rosen, E. Gurrola, G.F. Sacco, and H. Zebker. The insar scientific computing environment. In EUSAR 2012; 9th European Conference on Synthetic Aperture Radar, pages 730–733, 2012.

[8] F. Thollard, D. Clesse, M.P. Doin, J. Donadieu, P. Durand, R. Grandin, C. Lasserre, C. Laurent, E. Deschamps-Ostanciaux, E. Pathier, E. Pointal, C. Proy, and B. Specht. Flatsim: The form@ter large-scale multi-temporal sentinel-1 interferometry service. Remote Sensing, 13(18), 2021.

[9] J. Wasowski and F. Bovenga. Investigating landslides and unstable slopes with satellite multi temporal interferometry: Current issues and future perspectives. Engineering Geology, 174:103–138, 2014.

[10] Y. Zheng, H. Fattahi, P. Agram, M. Simons, and P. Rosen. On closure phase and systematic bias in multilooked sar interferometry. IEEE Transactions on Geoscience and Remote Sensing, 60:1–11, 2022.

SAR interferometry is a well-established technology, recently applied even at continental scale [1], for monitoring ground motions with millimeter-scale precision from time series of satellite SAR acquisitions. A key step of the technology is the identification, among billions of pixels, of the points (typically corresponding to man-made structures, rocks, or bare soil) that exhibit interferometric phase coherence for all the stack acquisitions. We will call here these points persistent scatterers (PS), regardless of the scattering mechanism that can be point-like or distributed. PS identification is not trivial (e.g., due to atmospheric and other systematic disturbances affecting the phase), and several techniques have been developed for the identification of PSs, based on statistics of the stack image amplitudes (amplitude dispersion, signal-to-clutter ratio) and/or phases in the spatial and/or temporal domains.

In this work, a novel algorithm, which we will call point coherence estimation (PCE), is presented for identifying PSs in a clean and simple way, without the need for spatial averages, amplitude/phase calibrations, or critical assumptions/approximations. PS selection is based on a novel technique we devised to estimate the temporal coherence (related to the phase noise) of each single point of the considered interferometric data stack from the coherences between pair of points, which can be directly calculated.

Let us start considering the phase differences between neighboring (within a few tens or hundreds of meters) points. The temporal coherences of these pairs-of-points can be determined since, as well known, the spatially correlated components (such as atmospheric and orbital artifacts, large scale motions) cancel out in the phase differences, whereas the temporally correlated components (i.e. the differences between the elevations and velocities – or higher order motion models – between the two points are estimated by maximizing the temporal coherence. Hence, the temporal coherence of each pair-of-points mainly depends only on the phase noises (e.g. temporal, spectral, geometric decorrelations, thermal noises) of the two points.

In the hypothesis that the phase noises in the two neighboring points of each considered pair are statistically independent (which might require to exclude pairing the nearest neighboring pixels if the images are oversampled), it can be easily demonstrated that the expected value of the temporal coherence for each pair is equal to the product of the temporal coherence expected values for the two paired points (analogous relations can be obtained also considering a finite number of samples instead of the theoretical expected values). Then, taking the logarithm of the obtained equations, an overdetermined system of linear equations is obtained. The overdetermined system can be solved by means of existing efficient solvers, with the solution corresponding to minimize, typically according to the L₁ or L₂ norm, the residuals of the equations. A reliable and consistent estimate of the temporal coherence of each single point is then obtained, based on which the PSs can finally be identified.

It is worth noting that the method does not need any assumption about the probability distribution of the phase noise. However, when considering a Gaussian probability distribution, the above system of equations states that the noise variance of each phase difference between a point pair is the sum of the noise variances of the two points.

The PCE method we propose can be applied to full-resolution data as well as to data with degraded resolution for a previous processing such as a multi-look or distributed scattering processing. Moreover, it is important to note that the method is quite stable, in the sense that applying the algorithm to the whole set or to a subset of the points produces similar results. In fact, it is possible, and it can be convenient in some cases, to iteratively apply the method to the previously selected points. In addition, this stability makes it possible also to apply our technique to a set of preselected candidate PSs, obtained for example by the AD and SCR methods with very relaxed thresholds, or by other quick techniques. This would reduce the computational time, which in any case is absolutely affordable. In fact, in addition to the calculation of the temporal coherence, which is common to all PS methods, the proposed algorithm requires the solution of an overdetermined system of linear equations, for which very efficient solvers can be used.

In all the tests performed, the method has proved to be very effective, providing for each single point a reliable measure of the temporal coherence, and of the related phase noise variance, therefore making it possible to detect a very large number of coherent points, i.e. PSs, with very few false detections.

We describe in the following some tests performed on two stacks of Sentinel-1 interferometric SAR images acquired over a pre-alpine area in Piemonte, Italy (86 acquisitions from January 2020 till October 2022) and over an area between Sicily and Calabria, Italy, including the Etna volcano (120 acquisitions from January 2020 till December 2022), respectively. The analyzed areas are affected by different kinds of displacement phenomena associated to natural and anthropic activities, and include different types of land cover, among which continuous and discontinuous urban fabric, transport infrastructures, bare soil, agricultural fields, mountains, and a big volcano.

The quality of the obtained results can be clearly appreciated by visual inspection of the selected PSs vs a very high resolution optical image of the ground. Moreover, we show some comparisons with the PSs identified by the classical, although basic, methods of the amplitude dispersion (AD) and/or signal-to-clutter ratio (SCR). Choosing thresholds such that the three methods approximately have the same PS false detection rate, our method provides significant improvements not only to the PS density, but also to the PS coverage of the ground, i.e. more areas and objects are covered by PS measurements.

In order to better clarify the difference between the three considered PS identification methods, we also computed 2D histograms relating AD, SCR and the temporal coherence estimated by our method. Their analysis shows that our algorithm is able to identify also PSs characterized by low SCR and high AD, confirming the effectiveness of the method proposed in this work.

One of the major products of interferometric synthetic aperture radar are displacement time-series which are of high importance in various applications including but not limited to monitoring the dynamics of volcanic activity, landslides, subsidence, earthquakes, ice and water. Different techniques have been developed and used to obtain displacement time-series from a stack of SAR images based on the type of scatterers which includes the Permanent Scatterers (PS) and Distributed Scatterers (DS). The most common techniques to estimate displacement time-series over DS pixels are based on classic Small BAseline Subsets (SBAS) and Phase linking methods. The former uses a subset of interferometric pairs and the latter uses all possible interferometric pairs to estimate displacement time-series. Recent studies have shown that Phase linking algorithms with full covariance matrix results in unbiased or less biased estimates of ground displacement compared to the SBAS algorithm. However, estimating displacement time-series from all possible interferometric pairs (full covariance matrix) is computationally expensive. A sequential estimator proposed by Ansari et al (2017) provides an efficient algorithm for processing the full covariance matrix in batches.

The current and future availability of dense Synthetic Aperture Radar (SAR) data from Sentinel-1 and upcoming NISAR missions has sparked the need to efficiently produce unbiased ground displacement time-series at fine resolution and in near real time. With the unprecedented InSAR big-data, producing displacement estimates for latest acquisitions with short latency (e.g., 24-72 hours from the acquisition time) requires novel algorithms to update the archived displacement time-series in contrast to reprocessing the entire archive. Although the sequential estimator is big-data friendly and potentially allows to update existing time-series with new acquisitions, it imposes a long latency of a few months to update existing time-series.

In this study, we propose an algorithm to update the InSAR displacement time-series with very short latency (few hours from the acquisition of new SAR data) without reprocessing the whole stack of the data. The algorithm is based on the phase linking approach and modifies the sequential estimator to meet a short latency of a few hours. In this algorithm we define a ministack as a subset of subsequent images with a size of N in which it may grow up to 2N-1 with new acquisitions, after which the latest minstack shrinks back to N. At each shrinking stage of the ministacks a compressed SLC which is a linear transformation of all SLCs in that latest ministack is estimated and used to form interferograms between the actual and compressed SLCs, i.e., ensuring the contribution of the long temporal baseline interferograms into the estimation of displacement at each acquisition. With this technique, only a limited amount of data will be pulled for the analysis and that includes the previous compressed images and the growing ministack with the size varying from N to 2N, therefore the computational efficiency improves.

In order to verify the near real time time-series algorithm, we simulate a displacement time-series with different decorrelation scenarios including long-term coherent, long-term decorrelated, light seasonal decorrelated and strong seasonal decorrelated, and we calculate the residuals obtained from near real time InSAR time-series technique and compare with the traditional sequential estimator. By comparing the estimated displacement time-series with the simulated displacement, the simulation results indicate low residuals for long term coherent as well as light and strong seasonal decorrelation. For the long-term decorrelated scenario where targets lose coherence rapidly over time, both real time and traditional sequential estimators show large residuals.

We also apply the above mentioned techniques to a stack of real data over a small region near Bristol dry lake in California which is known for the systematic closure phase bias when processing with conventional small baseline approach. We compare the performance of the sequential EMI method with the near real time algorithm. We also evaluate the performance of the time-series estimation when the stack is divided into two parts such that the first half is processed with traditional sequential EMI and the second half processed with the near real time estimation algorithm.

The results from real data demonstrate that the displacement time-series from the near real time algorithm is comparable with the traditional sequential EMI and significantly less biased compared to conventional SBAS algorithm.