PSI and DInSAR (1)

Summary:

The three geodetic techniques InSAR, levelling and GNSS are complementary w.r.t. (i) the sensitivity to horizontal and vertical displacement components, (ii) the spatial and temporal resolution of the measurements and (iii) the accuracy of the resulting displacements. We present a strategy to robustly combine the displacement estimates of the three techniques with a focus on long-term, secular movements. In addition, a comparison of the InSAR displacement rates with results obtained from levelling and GPS can serve as a validation of the sensitivity of InSAR.

Database:

Having access to a large amount of geodetic data from InSAR, levelling and GNSS, our study is carried out in the tectonically interesting Upper Rhine Graben (URG) region. The URG is located in the border triangle between Germany, France and Switzerland and is the most prominent part of the European Cenozoic Rift system. In recent years, the area is characterised by moderate tectonic and seismic activity, but on the other hand, carrying a significant probability for strong earthquakes. As the complex fault system of the URG is capable of both extensional and strike slip faulting, vertical as well as horizontal surface movements are expected, albeit on a small scale (< 1 mm/a). We use all available data sets from InSAR (ERS and Envisat), precise levelling and GNSS in order to estimate linear surface motions in the URG area (see attached figure).

SAR data:

SAR data acquired by ERS-1/2 and Envisat covering a period from 1992 to 2000 and 2002 to 2010, resp., both from ascending and descending tracks are analysed using the StaMPS (Stanford Method for Persistent Scatterers) software. Focussing on small displacement rates, we calculate quality indicators at every PS point in order to filter the resulting PS points of each analysed image stack. As a linear displacement rate gets more accurate with a longer time span, we combine the displacement time series of ERS and Envisat for both image geometries. Therefore, spatial interpolation of the PS points of every interferogram is needed, as the locations of ERS and Envisat PS points do not coincide. For the calculation of linear velocities from a combined ERS-Envisat time series, we consider both individual accuracies of a PS point displacement in a specific interferogram as well as temporal correlations between the interferograms induced mainly by the atmospheric filtering. The InSAR time series combination can easily be extended to data acquired from other SAR sensors, such as Sentinel-1.

Levelling data:

The levelling database consists of repeatedly measured height differences along levelling lines in the URG area obtained by the surveying authorities of Germany, France and Switzerland. The lines build a network of measurements with different spatial distributions in each country. As the measurement dates of the levelling campaigns do not coincide even within one country an epoch-wise adjustment of the measurements is not possible. Therefore, we apply a kinematic adjustment on the raw measurements resulting in a linear displacement rate at every levelling benchmark. As levelling is a relative measurement technique the linear rates are relative to a reference point located in a presumedly stable area.

GNSS data:

The transnational cooperation project GURN (GNSS Upper Rhine Graben Network) consists of more than 80 permanently operating GNSS sites in the URG area. At present, we apply a differential processing on a network of sites using GPS observations only. For every site, daily coordinate solutions are estimated, resulting in time series for Easting, Northing and Up components. We use the CATS (Create and Analyse Time Series) software package to handle jumps and periodic variations in the GPS time series finally resulting in linear velocity rates w.r.t. ITRF 2008. The vertical component of GPS is not used in the combination, as it is a factor of 3 worse than the horizontal components and a vertical velocity field with a higher spatial resolution is available from the levelling analysis.

Combination approach:

The data sets from InSAR, levelling and GPS are inhomogeneous in space and time. When combining the displacement results of the techniques, we implicitly assume that (i) InSAR, levelling and GPS measure the same signal, (ii) the tectonic displacements behave linearly in time and (iii) the tectonic signal is smooth in space and does not abruptly change within short distances. The latter requirement is of special importance as we use Kriging to interpolate the different data sets to a common location. Most of the non-linear movements in the area of investigation could be assigned to anthropogenic deformations which have been investigated in several case studies.

As all three techniques are based on relative measurements, they depend on a reference point (levelling), a reference area (InSAR) or a reference frame (GPS), which are located in different regions. We firstly calculate offsets between InSAR and levelling for the vertical component and between InSAR and GPS for the horizontal component in order to shift the rates obtained from InSAR into the datum of levelling and GPS. Therefore, the PS displacements are individually interpolated to the location of levelling benchmarks and GPS sites. From PS velocity rates in ascending and descending image geometries, a mean vertical and horizontal offset between InSAR and the other two techniques is calculated and applied later in the combination step. For the combination, we interpolate the displacement rates of all the three techniques to a common grid. In order to avoid spatial extrapolation, this grid only carries values in the vicinity of existing PS points. At every valid grid point, the linear LOS velocity rates from PS-InSAR (ascending and descending) are combined with vertical velocity rates from levelling and horizontal velocity rates from GPS resulting in a 3D velocity vector (Northing, Easting, Up). A special focus in the least squares estimation of the three displacement components is given to realistic information on the variances and covariances of the input data.

Results:

Besides the theoretical aspects of our approach to combine surface displacements measured by InSAR, levelling and GNSS, we will present the resulting displacement fields for the URG area. The small horizontal and vertical linear rates recovered from our geodetic analysis are in accordance with geological and geophysical models and indicate that our combination approach is able to resolve tectonic movements at the sub-mm/a level  with a formal error below 0.5 mm/a. Nevertheless, some non-tectonic features are still visible induced by groundwater usage, oil extraction or mining activities. In our presentation, the combined 3D velocities in the URG region will be shown as well as local case studies of man-made surface displacements.

At Fringe 2009 Fabrizio Novali introduced SqueeSAR, a refined approach to deformation analysis developed by Ferretti, Fumagalli, Novali himself, Prati, Rocca and Rucci. The presented comparisons of results generated with SqueeSAR to such generated with PSInSAR exhibited an impressive increase of extracted information. This success is based on a unified spatio-temporal statistical analysis which can be seen to turn DS (distributed scatterers) into PS. Means of this conversion is a maximum likelihood estimator, derived and theoretically investigated by Monti Guarnieri and Tebaldini (2008), that provides the complex signal common to a group of DS pixels obeying a joint circular complex normal distribution. In the framework of SqueeSAR, published 2011 in TGRS, this estimator is made applicable by providing ideas for grouping DS candidate pixels in a neighborhood of an investigated pixel, which is supposed to be statistically homogeneous and to allow for the estimation of the covariance matrix.

Despite the demonstrated capability of this approach to enhance coverage with high quality information considerably and although the use of related ideas in the SAR interferometry community increased in the last years, questions which arise naturally when one tries to put to work this approach are still not answered satisfyingly. These questions concern criteria for grouping pixels, estimators for the covariance matrix and the use of characteristic numbers, which allow recognizing low quality pixels and cutting short time consuming computations for them. Moreover good DS pixels have to be identified in order to be able to decide if they can be included in the start set of pixels on which an initial estimation is performed. This is of particular interest for rural scenes to bridge gaps in the PS net.

For the purpose of exploring these questions we implemented a testing environment for the approach presented in the SqueeSAR paper (2011). DS are processed in a detail of the investigated scene and then connected to a reference result of a PS analysis. Several criteria for grouping pixels and alternatives for estimating the covariance matrix are considered. Different characteristic numbers available at different steps of the algorithm are discussed. Pixels are assessed according to the temporal coherence of their phase differences to nearby PS and versus ground truth. Test data are from a scene in Bavaria, where a LIDAR DEM is available and from the town of Lüneburg, Germany, where we use deformation series from leveling as ground truth. Both stacks are TerraSAR-X high resolution spotlight-mode (300MHz).

Ferretti, A., Fumagali, A., Novali, F., Prati, C., Rocca, F., Rucci, A. (2009): “The Second Generation PSInSAR Approach: SqueeSAR”, Fringe 2009, ESA-ESRIN, Frascati.

Ferretti, A., Fumagali, A., Novali, F., Prati, C., Rocca, F., Rucci, A. (2011): “A New Algorithm for Processing Interferometric Data Stacks: SqueeSAR”, IEEE Transactions on Geoscience and Remote Sensing, Vol. 49, Issue 9.

Monti Guarnieri, A., und Tebaldini, S. (2008), “On the exploitation of targets statistics for SAR interferometry applications”, IEEE Trans. Geosci. Remote Sens., Vol. 46, No. 11, pp. 3436–3443, Nov. 2008.

1. Introduction

Synthetic aperture radar (SAR) interferometry and radargrammetry are two well-known technologies for generating digital surface models (DSMs) from SAR data acquired from satellite or airplane. However, standard SAR interferometry and radargrammetry have strong limits for the derivation of DSMs over urban areas, due to the complexity of urban scenarios and the need for very high precision to reconstruct building shapes. Therefore, DSM or 3D models of urban areas are currently produced practically only based on photogrammetric or Lidar technologies. However, there are important reasons to consider also SAR data for the realization of urban DSMs. In fact, several areas in the world are not accessible by photogrammetric or Lidar technologies, due to frequent cloud coverage, or political or economical reasons. On the contrary, the existence nowadays of several high resolution satellite SAR missions can make available valuable data all over the world without the above mentioned limitations. In particular, the four satellites of the COSMO-SkyMed mission have been acquiring series of interferometric images over all main cities and sites of interest of the world.

In this work, we propose the realization of urban DSMs from high resolution satellite SAR data: advanced interferometry techniques can provide a cloud of points (persistent scatterers – PS) localized with metric precision in the 3D space; then surface reconstruction techniques applied to this point cloud are used to calculate the urban surface model. For the first step it is important to use advanced interferometric techniques able to provide a cloud of points as dense as possible. The second step remains critical because the PS points are in any case very sparse and not distributed homogeneously. However, a priori information on the distribution of the PS and on the structure of the urban scenario can be exploited to help the surface reconstruction.

2. Interferometry techniques

Persistent scatterer interferometry, founded on the idea (first introduced in [1]) of analyzing long series of full resolution SAR images, allows detecting and locate in the 3D space a set of sparse points that remain correlated at different acquisitions, named persistent scatterers (PS). More recently, we developed the persistent scatterer pair (PSP) method [2-4], based on the idea of both identifying and analyzing PS working only with pairs of points
(“arcs”), which makes the methods insensitive to spatially correlated signals such as atmospheric or orbital artifacts. In particular, the enhanced PSP algorithm searches for an optimized set of arcs connecting pair of points in order to fully test all the image pixels as potential PS, and significantly improves the already very good performance of the standard PSP approach. This makes possible to identify (and locate in the 3D space with metric precision) a large number of PS, corresponding to practically all SAR image pixels where you could expect a coherent signal, even when there are not strongly scattering structures and the sensed signal is low. The method is capable to select PS also in parts of the structures where standard techniques are not able to extract information (see Figure 1), and is very useful to build dense point clouds from which to reconstruct urban DSMs.

In order to assess the precision of the PSP measurements from COSMO-SkyMed data, we evaluated the statistics of the height measurements obtained from of 30 COSMO-SkyMed ascending Stripmap SAR images acquired from February 2011 to June 2013 in an area near Novara, Italy, with baselines spanning 1590 m on the flat roofs of two industrial buildings (see Figure 2). Results show that the standard deviation of the height measurements is about 0.7 m.

It is worth mentioning that another approach to consider in the problem of collecting point clouds for urban DSM is SAR tomography [5-6], a technique that allows detecting (and locate in the 3D space) PS points by processing the whole complex information of the SAR images. The tomographic technique requires a preliminary calibration of the data (that can be performed for example also with the PSP method), in particular to remove atmospheric components from the phase. Tomography is particularly important to analyze buildings in detail, since allows detecting different scattering mechanisms mixed in a single pixel of the SAR image, a frequent circumstance in urban areas.

3. surface reconstruction

The problem of surface reconstruction from point clouds has several applications (e.g., medical scanners, laser range finders, vision techniques, Lidar). The general problem is known as surface reconstruction from unorganized points, and several methods have been proposed to solve it [7], [8].

We have performed several experiments, on simulated and real data, for obtaining a 3D model from SAR-derived cloud of points, testing different surface reconstruction algorithms among which, Screened Poisson [9], Power Crust [10], Hoppe et al. [11]. Ball pivoting is a simple but effective algorithm that we mostly used for our preliminary tests. It computes a triangle mesh interpolating a given points cloud. The principle is very simple: three points form a triangle if a ball of a user-specified radius r touches them without containing any other point. Starting with a seed triangle, the ball pivots around an edge (i.e. it revolves around the edge while keeping in contact with the edge’s endpoints) until it touches another point, forming another triangle. The process continues until all reachable edges have been tried, and then starts from another seed triangle, until all points have been considered. The relatively small amount of memory required by the BPA, its time efficiency, and the quality of the results obtained compare favorably with existing techniques.

In general, for all types of surface reconstruction algorithms, the surface reconstruction from the PS point clouds generated by interferometry techniques is very difficult, because the PS points are very sparse and not distributed homogeneously. In order to improve the results, it is important to exploit the highest SAR resolution data and it is useful to exploit some a priori information on the distribution of PS (e.g., their tendency to accumulate along the borders of roofs and roads), and on the structure of the surface to be reconstructed (mainly vertical and almost horizontal surfaces).

4. First tests

In order to verify the potential of the proposed method some preliminary tests were performed using COSMO-SkyMed SAR data. The results (see Figure 3), obtained using stripmap data (3 m ground resolution) acquired over Rome, Italy, demonstrate the high potential of the proposed method for building urban DSM in areas of the world not easily accessible by photogrammetric or Lidar technologies (e.g. due to frequent cloud coverage, or political or economical reasons). More complete tests are in progress and will be included in the final paper and in the presentation. In particular, the use of very high resolution COSMO-SkyMed spotlight SAR data (1 m ground resolution) will greatly improve the quality of the obtained urban DSM. Moreover, it is expected that the performance of surface reconstruction techniques can be improved by introducing in the processing a priori information on the distribution of PS and the typical structure of the urban scenario.

4. References

[1] A. Ferretti, C. Prati, and F. Rocca, “Non-linear subsidence rate estimation using permanent scatterers in differential SAR interferometry,” IEEE Trans. Geosci. Remote Sensing, vol. 38, pp. 2202–2212, Sept. 2000.

[2] M. Costantini, S. Falco, F. Malvarosa, F. Minati, “A new method for identification and analysis of persistent scatterers in series of SAR images” in Proc. Int. Geosci. Remote Sensing Symp. (IGARSS), Boston MA, USA, pp. 449-452, 2008.

[3] Mario Costantini, Federico Minati, Francesco Trillo, Francesco Vecchioli, “Enhanced PSP SAR Interferometry for Analysis of Weak Scatterers and High Definition Monitoring of Deformations over Structures and Natural Terrains”, in Proc. Int. Geosci. Remote Sensing Symp. (IGARSS), Melbourne, Australia, 2013.

[4] M. Costantini, S. Falco, F. Malvarosa, F. Minati, F. Trillo, F. Vecchioli, “Persistent scatterer pair interferometry: approach and application to COSMO-SkyMed SAR data,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 7, no. 7, July 2014.

[5] F. Lombardini, “Differential tomography: a new framework for SAR interferometry,” IEEE Transactions on Geoscience and Remote Sensing,  (Volume:43 ,  Issue: 1 ), pp. 37-44, 2005.

[6] D. Reale, G. Fornaro, A. Pauciullo, X.Zhu, and R. Bamler, “Tomographic Imaging and Monitoring of Buildings With Very High Resolution SAR Data,” IEEE Geosci. Remote Sens. Letters, vol. 8, no.4, pp. 661–665, July 2011

[7] J. Hrádek, “Methods of surface reconstruction from scattered data,” Technical Report No. DCSE/TR-2003-02, University of West Bohemia in Pilsen, March 2003.

[8] M. Berger, A. Tagliasacchi, L. M. Seversky, P. Alliez, J. A. Levine, A. Sharf, C. T. Silva, “State of the Art in Surface Reconstruction from Point Clouds,” Eurographics 2014.

[9] M. Kazhdan, M. Bolitho, H. Hoppe, “Poisson Surface Reconstruction,” in Proceedings of the fourth Eurographics symposium on Geometry processing, 2006.

[10] N. Amenta, S. Choi, R. K. Kolluri, “The Power Crust, Union of Balls and Medial Axis Transform,” Computational Geometry, 19(2) (2001), 127-153.

[11] H. Hoppe, T. DeRose, T. Duchampy, J. McDonaldz, and W. Stuetzlez, “Surface reconstruction from unorganized points,” in Proceedings of the 19th annual conference on Computer graphics and interactive techniques, 1992.

The paper will describe a new approach to Persistent Scatterer Interferometry (PSI) data processing and analysis, which is implemented in the PSI chain of the Geomatics (PSIG) Division of CTTC. The PSIG chain includes three main processing blocks. In the first one a set of correctly unwrapped and temporally ordered phases are derived, which are computed on Persistent Scatterers (PSs) that cover homogeneously the area of interest. The key element of this block is given by the so-called Cousin PSs (CPSs), which are PSs characterized by a moderate spatial phase variation that ensures a correct phase unwrapping. This block makes use of flexible tools to check the consistency of phase unwrapping and guarantee a uniform CPS coverage. In the second block the above phases are used to estimate the atmospheric phase screen. The third block is used to derive the PS deformation velocity and time series. The key tool of this block is a new 2+1D phase unwrapping algorithm.

The inputs of the PSIG include a stack of N co-registered SAR images, the amplitude dispersion (DA) and M wrapped interferograms, with M>>N. The main processing steps are briefly described below.

1) Candidate Cousin PS (CPS) selection. In this step, a set of PSs with phases characterized by a moderate spatial variation is sought. This is accomplished by using at least a seed PS and searching for its “cousins”, i.e. PSs with similar characteristics.

2) 2D phase unwrapping. 2D phase unwrapping is performed on the candidate CPSs using the redundant set of M interferograms.

3) Phase unwrapping consistency check. This check is based on a least squares estimation, followed by the analysis of the so-called residuals. The final set of CPSs is selected at this stage.

4) APS estimation and removal. Using the selected CPSs, the APS is estimated and then removed from the original interferograms, obtaining a set of M APS-free interferograms.

5) Estimation of deformation velocity and residual topographic error (RTE). The deformation velocity and RTE are computed over a dense set of PSs (much denser than the selected CPSs), from the M wrapped APS-free interferograms, using the method of the periodogram. Optionally, an extension of the two-parameter model can be used to account for the thermal expansion.

6) RTE removal. The RTE phase component is removed from the wrapped APS-free interferograms. The linear deformation component can optionally be removed and then, in a later stage, added back to the deformation time series. The same procedure can be done with the thermal expansion component.

7) 2+1D phase unwrapping. A 2+1D phase unwrapping is executed on the set of M APS- and RTE-free interferograms to obtain the final deformation phase time series, a quality index for each time series and other parameters related to the detection and correction of unwrapping errors.

The performance of the procedure will be illustrated using X-band data and, if possible, with the Sentinel-1 data.

The Wide Area Product (WAP) is a new interferometric product developed to provide measurement over large regions. Persistent Scatterers Interferometry (PSI) has largely proved their robust and precise performance in measuring ground surface deformation in different application domains. This technology was validated within a PSI certification activity included in the Terrafirma project, funded by the European Space Agency, where PSI results produced by different processing chains were validated against them and ground data. The inter-comparison reported that precisions below 1mm/year are achieved and the Root Mean Square Error (RMSE) ranges from 0.8 to 1.5 mm/year, based on ERS and ENVISAT data and ground survey. This validation demonstrates the great accuracy of PSI technique for measuring ground surface deformation in relatively local (20km) and pilot test areas (urban and rural areas with non-pronounced topography).  

In this context, however, the accurate displacement estimation over large-scale areas (more than 10.000km²) characterized by low magnitude motion gradients (3-5 mm/year), such as the ones induced by inter-seismic or Earth tidal effects, still remains an open issue. The main reason for that is the inclusion of low quality and more distant persistent scatterers in order to bridge low-quality areas, such as water bodies, crop areas and forested regions. This fact yields to spatial propagation errors on PSI integration process, poor estimation and compensation of the Atmospheric Phase Screen (APS) and the difficult to face residual long-wavelength phase patterns originated by orbit state vectors inaccuracies.  

Research work for generating a Wide Area Product of ground motion in preparation for the Sentinel-1 mission has been conducted in the last stages of Terrafirma as well as in other research programs. These developments propose technological updates for keeping the precision over large scale PSI analysis. Some of the updates are based on the use of external information, like meteorological models, and the employment of GNSS data for an improved calibration of large measurements.

The main characteristic of a Wide Area processor shall be the capacity of concatenating parallel orbit tracks and consecutive image frames to avoid any constrain in terms of ground coverage. This is a key factor and represents the novelty of the product for particular applications like tectonic studies where maps at regional scale should be provided. Despite covering wide regions, the retrieval of large scale motions can, however, not be possible depending on the processing approach. In terms of data processing, the tests performed showed that several steps may condition the final result and should thus be carefully taken into account:  the spatial resolution of the interferograms, the point network selection, the integration from a unique seed or based on several distributed seeds, and the APS estimation and calibration between different tracks among others. Depending on the configuration and performances of the different PSI processing steps, large trends or local episodes may be underestimated; therefore a standard processing is yet to be defined.

The present work wants to contribute to this open topic by exposing the experience of the authors in PSI Wide Area processing over different thematic fields. The lessons learnt and the main limitations along the way will be described and analyzed. Furthermore, the different WAP strategies developed by the authors will be discussed in order to maximise the trade-off between large-scale motion and local phenomena monitoring. Conclusions and recommendations will be extended in terms of the Sentinel-1 data. The default mode of Sentinel-1, the Interferometric Wide Swath Mode, as well as the regular acquisitions plan over tectonic belts foreseen by ESA and other initiatives like COMET+, highlight the significance of this topic. Consequently, it is important to set the state of the art, identify requirements and limitations and propose solutions to them.

During the round table, seed questions proposed by the chairs will be discussed with the audience.