Atmospheric corrections in interferometric synthetic aperture radar surface deformation – a case study of the city of Mendoza , Argentina

Differential interferometry is a remote sensing technique that allows studying crustal deformation produced by several phenomena like earthquakes, landslides, land subsidence and volcanic eruptions. Advanced techniques, like small baseline subsets (SBAS), exploit series of images acquired by synthetic aperture radar (SAR) sensors during a given time span. Phase propagation delay in the atmosphere is the main systematic error of interferometric SAR measurements. It affects differently images acquired at different days or even at different hours of the same day. So, datasets acquired during the same time span from different sensors (or sensor configuration) often give diverging results. Here we processed two datasets acquired from June 2010 to December 2011 by COSMO-SkyMed satellites. One of them is HH-polarized, and the other one is VV-polarized and acquired on different days. As expected, time series computed from these datasets show differences. We attributed them to non-compensated atmospheric artifacts and tried to correct them by using ERAInterim global atmospheric model (GAM) data. With this method, we were able to correct less than 50 % of the scenes, considering an area where no phase unwrapping errors were detected. We conclude that GAM-based corrections are not enough for explaining differences in computed time series, at least in the processed area of interest. We remark that no direct meteorological data for the GAM-based corrections were employed. Further research is needed in order to understand under what conditions this kind of data can be used.


Introduction
Differential synthetic aperture radar interferometry (DIn-SAR) is an advanced remote sensing technique aimed to large-scale surface deformations monitoring with centimeter to millimeter accuracy, by exploiting the round-trip phase components of SAR images relative to an investigated area (Massonnet and Feigl, 1998).Information of ground deformation is associated with the phase difference between two acquisitions, referred to as master and slave images of the interferometric data pair.Both acquisitions must be acquired from relatively close tracks (spatial baseline) and acquisition time (temporal baseline) in order to reduce the temporal and geometric decorrelation phenomena and topographic errors (Berardino et al., 2002).
An effective way for studying and understanding the dynamics of the deformation phenomena and their temporal behavior is the generation of deformation time series.Multitemporal InSAR techniques were developed recently; they are based on combining information obtained from multiple SAR images acquired over a period of time.The most widely used advanced DInSAR algorithms are small baseline subsets (SBAS) approach (Berardino et al., 2002) and persistent scatterer (PS) (Ferretti et al., 2001).SBAS technique is based on an adequate combination of the differential interferograms characterized by a small spatial and temporal separation (spatial and temporal baseline) in order to limit the geometric and temporal decorrelation phenomena and to maximize the number of coherent pixels exploited.Its capability of detecting and investigating long-time deformation phenomena has been already shown in different applications based on exploiting ERS and ENVISAT data of the European Space Agency (ESA) (Pepe et al., 2005;Tizzani et al., 2007).Furthermore, deformation time-series generation capability of the small baseline subset algorithm to the Radarsat-1 data has been demonstrated (Euillades et al., 2009).
Phase propagation delay in the atmosphere is the main systematic error of interferometric SAR measurements.The dominant contribution of the atmospheric phase delay, which may reach tens of centimeters, comes from the temporal variation of the stratified troposphere (Hanssen, 2001).Such variations are related to changes in variables such us temperature, atmospheric pressure and water vapor.In particular, atmospheric water vapor effects represent one of the major limitations for repeat-pass InSAR, and limit the accuracy of deformation rates derived from DInSAR (Li et al., 2006).Due to their dependence on altitude, such variations may produce a spatial correlation between the magnitude of the phase SAR signal and the topographic elevation.
Several methods have been recommended for estimating atmospheric phase delay corrections and isolating displacements from atmospheric artifacts.For example, Onn and Zebker (2006) use zenithal wet delay observations from GPS Network.Li et al. (2006) employ satellite multispectral imagery analysis, and Jolivet et al. (2011) profit from global atmospheric model (GAM) derived data.
The atmospheric contribution becomes potentially greater as the SAR wavelength becomes shorter (Hanssen, 2001).This is an issue when processing data acquired by the new generation X-band satellites like COSMO-SkyMed and TerraSAR-X.Conversely, crustal deformation monitoring can benefit from the high spatial and temporal resolution of these sensors.Furthermore, the fact that several SAR sensors are simultaneously orbiting Earth makes potentially available series of scenes covering a given area during a common time span.As systems are characterized by different polarization and/or acquisition geometry, and the scenes are eventually not acquired on the same days/hours, the emerging question is if the area of interest's deformation history can be accurately characterized by using DInSAR-based techniques.In other words, are time series computed from scenes acquired at different times with different systems or with the same system but different polarization giving the same results?As the underlying deformation is the same, eventual differences could be attributed to not well-compensated atmospheric contributions and/or processing (phase unwrapping) errors.
In this work we present the results of an experiment consistent in computing deformation time series by exploiting HH-and VV-polarized SAR data acquired during a common time span by COSMO-SkyMed satellites.Covered area is a strip of roughly 40 × 50 km located near the city of Mendoza in western Argentina.Obtained results show significant differences among the two elaborations, which motivate us to explore if estimating the atmospheric phase delay from ERA-Interim global atmospheric model (GAM) is useful as an operative correction tool.

Dataset and study area
We used two COSMO-SkyMed datasets composed of 31 HH-polarized and 27 VV-polarized images respectively.Acquisition geometry is the same for all scenes: stripmap mode (HIMAGE), descending pass with an incidence angle of 38 • .Both datasets cover the period between June 2010 and December 2011, but HH-and VV-polarized scenes were acquired on different days, as can be seen in Fig. 1.
Different land uses are well represented within the selected test site: (1) urban areas (Mendoza city), (2) agricultural coverage around the city, (3) piedmont and bare soil with low vegetation and (4) high relief (Andes Mountains).Area location and main characteristics are presented in Fig. 2.

A brief description of the SBAS algorithm
This brief description is based on Berardino et al. (2002).Considering n + 1 SAR images relative to the same area, acquired at chronological ordered time [t 0 , t 1 , . . .t n ] and assuming that all the images are co-registered, the DInSAR-SBAS algorithm begins with an adequate combination of image pairs, which allows minimizing geometrical and temporal decorrelation.A number m of multilook differential interferograms are generated.Unwrapped differential phase of the k-th interferogram, considered at the generic pixel of azimuth and range coordinates (x, r), is expressed by where k = 1, . . ., m; t A and t B are the acquisition times of master and slave images and ϕ(x, r, t B ) − ϕ(x, r, t A ) represents the phases of the two images involved in the interferogram generation, λ is the radar wavelength, d(x, r, t A ) and d(x, r, t B ) are the radar line-of-sight (LOS) projections of the cumulative surface deformation at the two times t A and t B .ϕ topo k (x, r) takes into account possible topographic artifacts z(x, r) that can be present in the digital elevation model (DEM) used for topographic phase component compensation.The term ϕ atm k (x, r, t A , t B ) accomplishes possible atmospheric disturbances between the acquisitions at times t A and t B , and it is often referred to as a atmospheric phase component.
n k (x, r) accounts for the noise effects.Finally, B ⊥,k represents the perpendicular baseline component and ϑ the SAR sensor look angle.
Considering phase velocities between adjacent acquisitions, Eqs. ( 1) and ( 2) allow defining a system of equations in the n + 1 unknowns: This can be solved by using the single value decomposition (SVD) as B v = δϕ in a matrix form.As a result, the SBAS algorithm gives displacement time series for each processed pixel.Atmospheric artifacts are filtered through the application of low pass filtering step in the two-dimensional spatial domain followed by a temporal high pass filtering.

DInSAR-SBAS processing
Both datasets were processed separately, and 85 HH and 78 VV differential interferograms were computed according to the SBAS restrictions on temporal and spatial baselines.The maximum perpendicular baseline allowed was 1000 m, which is one-third of the critical baseline for the acquisition geometry.However, the longest one effectively used was 751 m (24 December 2010-25 January 2011) for HH elaboration and 923 m (23 April 2011-9 May 2011) for VV elaboration.The maximum temporal baseline was 150 days.Spatial vs. temporal baseline distribution is shown in Fig. 3, where scenes and interferograms are represented by points and arcs, respectively.The 30 m X-band SRTM DEM (Rabus et al., 2003) was used for the topographic phase removal.HH and VV time series were computed by inverting the interferograms after phase unwrapping process with an extension of the minimum-cost flow (MCF) algorithm (Pepe and Lanari, 2006).Finally, atmospheric contributions are filtered out by using a set of filters in cascade as described in Tizzani et al. (2007).
Mean deformation velocity maps, computed from both elaborations, are shown in Fig. 4 (SAR projection, azimuth in the horizontal direction) where low coherent pixels have been masked out.Time series are referred to a pixel (RPreference point) located in Mendoza city (white triangle in Fig. 4), considered "stable" or point of zero deformation.We considered unreliable those pixels with temporal coherence (Tizzani et al., 2007) lower than 0.8.As expected, coherent pixel density is higher in urban and piedmont areas than in agricultural or high relief ones.
Obtained deformation sequences were co-registered for comparing the results.To do so we computed the azimuth and range shifts of the VV-polarized elaboration master scene As the acquisition times of HH and VV scenes are not coincident, one of the series (e.g., the VV one) must be interpolated for estimating the displacement at the acquisition time of the other one (i.e., the HH series), as illustrated in Fig. 5.For increased robustness we used a mean root mean square obtained by averaging Eq. (4a) of interpolated VV and Eq.(4b) of interpolated HH time series (see Eq. 4c):

Stratified atmospheric contribution correction
We computed a second set of HH and VV deformation time series by previously correcting the differential interferograms for stratified atmospheric contributions.Employed SBAS settings were exactly the same as already described in the previous section.
For computing the corrections, we used the Python-based Atmospheric Phase Screen mitigation library (PyAPS) (Jolivet et al., 2012), which uses meteorological information derived from global atmospheric models (GAMs) to estimate the dry and wet components of the atmospheric delay.In particular, we employed results from the ERA-Interim 2011 reanalysis computed by ECMWF (Dee et al., 2011).
ERA-Interim is the latest global atmospheric reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF).It is based on 4D-Var assimilation of global surface and satellite meteorological data (Dee et al., 2011).This re-analysis provides several meteorological parameters on a global ∼ 75 km grid from 1989 to the present.Such parameters are computed daily at 12:00, 06:00, 00:00 and 18:00 UTC.Vertical stratification is described on 37 pressure levels, densely spaced at low elevation (25 hPa), with the highest level around 50 km (1 hPa).The pressure levels located under the local elevation of grid nodes are obtained by extrapolation (Jolivet et al., 2011).
COSMO-SkyMed SAR images were acquired at ∼ 21:54 UTC.For each scene date, we got the closest ERA-Interim data (http://rda.ucar.edu/datasets/ds627.0/), which are those computed at 00:00 GMT.Input files (GRIB extension) corresponded to a World Meteorological Organization (WMO) standard for distributing gridded data.We used PyAPS for interpolating the temperature, water vapor and dry air partial pressure provided at each pressure level.This allows estimating the atmospheric delay as a function of elevation by using a mixed Clausius-Clapeyron law (Jolivet et al., 2012) on each ERA-Interim grid point located in the vicinity of the radar image.Next, a bilinear interpolation in the horizontal dimensions and a spline interpolation along altitude are applied for producing the tropospheric delay maps in radar geometry (full resolution).High-resolution tropospheric delay maps were multilooked to the COSMO-SkyMed scene processing resolution.Furthermore, we referenced them to the time of the first acquisition and to the reference point (RP).Subsequently, we generated the corrections for each one of the 2 π-wrapped differential interferograms and applied them to the original ones.
In Fig. 6    suggesting that they are due to topography-dependent water vapor effects.
The whole procedure employed for generating uncorrected and corrected time series is summarized in the flow chart displayed in Fig. 7.

Results and discussion
Mean deformation velocity maps computed as described in Sect.2.2 are shown in Fig. 4. The maps show significant differences.Near the reference point RP, in the urban area of Mendoza city, the velocity is almost zero in both elaborations.As we move away from the RP, noticeable differences in magnitude and sign of the LOS (line-of-sight) mean deformation velocities can be observed.In general, deformation velocity of the whole area, for the time span processed, is lower than 1.3 cm yr −1 in absolute value.However, localscale subsidence patterns of about 2 cm yr −1 were detected in both HH and VV elaborations (see Fig. 4) near Potrerillos Dam (a and b), and in Mendoza city (c and d).VV elaboration reports velocities systematically higher than HH ones in these areas.
As the mean deformation velocity gives only a quick look of the underlying data, which are the deformation time series, we based our analysis on them.We computed the average RMS map by using Eq. ( 4c).An amplitude scene, DEM, and average RMS map are shown in Fig. 8A-C, respectively.From the figure, a correlation between RMS map and topography is apparent: lower areas (around the Mendoza city) show RMS of 0 mm, which increases up to ∼ 35 mm in higher areas.One can detect roughly three areas of different Correlation between average RMS and topography, and the fact that high RMS values are due to differences between time series at some dates (e.g., 24 December 2010 to 14 March 2011), at least in the green RMS zone, motivated us to hypothesize that differences in computed time series are related to non-compensated stratified tropospheric phase components.In order to prove this idea, we corrected the differential interferograms and re-computed the deformation time series with the procedure already described in Sect.2.3.Obtained results are presented in Fig. 9, where RMS maps computed by Eqs.(4a) and (4b) for non-corrected-HH, -VV and corrected-HH, -VV time series are displayed.Clearly, RMSs for corrected time series are higher than RMSs for non-corrected ones in areas topographically higher, whereas it diminishes in lower areas, particularly those located in the near range: see highlighted box and time series A, B, C and D in the figure .In order to understand why the time series, and consequently the RMS metric, are heavily affected in the mountainous region (see VV RMS map and E, F time series in Fig. 9), we analyzed the unwrapping residuals.Unwrapping residuals compute the difference between differential interferograms and synthetic differential interferograms reconstructed from phase time series after SVD inversion (more details can be seen in Tizzani et al., 2007).If both interferograms are equal, which means that phase unwrapping was correct, residuals are near zero.Sometimes non-zero residual are found in relatively isolated areas in terms of coherence, thus indicating non-reliable unwrapping there.Figure 10 presents significant residual maps computed after the ERA-Interim VV-and HH-corrected elaborations.Note that non-zero residuals are found in areas that correlate well with areas of RMS worsening.This could explain the observed effect: high RMS is not necessarily the result of failed corrections but of processing errors in the phase unwrapping step.
At this point, it is interesting to check if the correlation between unwrapped phase and topography before and after correction systematically diminished.For that comparison we do not used the unwrapped interferograms but the unwrapped phase obtained after the SBAS inversion via single value decomposition.Note that the unwrapped phase map relevant for each acquisition can be understood as a synthetic interferogram computed between it and the first acquisition (in time) of the series.We computed the Pearson correlation coefficient, which improved in 8 HH and 13 VV scenes after correction, whereas it became worse in 23 HH and 14 VV scenes.If, instead of computing correlation in the whole area, we restrict it to near the range of low RMS zone (white box in Fig. 9), the results are as follows: Pearson's coefficient improved in 17 HH and 9 VV scenes after correction, whereas it become worse in 14 HH and 18 VV scenes.The correction seems to work in the HH time series, but fails in the VV one.
In summary, many interferograms show topography correlated fringes, which are an indication of atmospheric stratification components in the differential phase.Corrections computed from ERA-Interim GAM remove the correlation between phase and topography in some scenes, but not in all of them.We also detected unwrapping errors affecting the high relief areas, so they should not be taken into account for atmospheric correction statistics.But inside the area without unwrapping errors, corrections perform well only in the HH elaboration.Possibly, corrections derived from a GAM in an area without enough direct meteorological data do not accurately represent the atmospheric state.Errors could be also derived from using GAM data computed for a time that differs in a few hours from the SAR passing time.

Conclusions
We computed deformation time series from HH-and VVpolarized COSMO-SkyMed SAR images acquired during a common time span.Algorithm used was DInSAR-SBAS.Results show noticeable differences between both elaborations, which is clearly a processing artifact given that the underlying deformation is the same.As scenes were acquired on different days, we explore using atmospheric data for correcting them.
We detected correlation between unwrapped phase and topography in both elaborations.Furthermore, differences between time series, quantified with a RMS metric, also show correlation with elevation: differences increase with height.
We applied corrections for stratified atmospheric components, computed by using independent data provided by ERA-Interim GAM.From the obtained results, we conclude that they are not enough for minimizing difference between the time series.In fact, the correction improves some of the interferograms (by removing topography-correlated fringes) but worsens others (increasing the correlation).This could suggest that GAM models fail to represent accurately the atmospheric state in this region, at least for some of the acquisition dates.Lack of direct meteorological data could explain this behavior, but further work for verifying this hypothesis needs to be done.

Fig. 1 .
Fig. 1.Distribution of HH and VV data acquisitions.Most of the acquisitions were alternately selected in HH and VV polarization mode in 8 days.

Fig. 3 .
Fig. 3. COSMO-SkyMed SAR data representation in the temporal/perpendicular baseline plane for HH (A) and VV (B) acquisitions.In both, the arcs correspond to the generated differential interferograms.

Fig. 4 .
Fig. 4. HH and VV mean deformation velocity maps superimposed on the amplitude images of the study area of the city of Mendoza.(RP) identifies the reference point of the start the unwrapping process.(A) and (C) represent a singular area in HH results, and (B) and (D) in VV results.deformation time series: (A) and (B) near the Potrerillos Dam; (C) and (D) in the city of Mendoza.
we show an example of the computed atmospheric corrections for the differential interferogram composed by acquisitions 15 June 2010 (master) and 1 July 2010 (slave).Note the topography-correlated fringes are indicated by two arrows in Fig. 6a and the stratified delay maps predicted by ERA-Interim in Fig. 6c.Finally note in Fig. 6d how these signals were significantly reduced after correction,

Fig. 5 .
Fig. 5. Example of deformation time series (DTS).White triangle represents HH displacement, asterisk VV deformation and black triangle VV DTS interpolated at HH series.

Fig. 6 .
Fig. 6.Example of a COSMO-SkyMed interferogram and atmospheric correction over the area of interest using the ERA-I global atmospheric model.All images are in radar geometry.(A) Wrapped interferogram from SAR acquisitions on 15 June 2010 (master) and 1 July 2010 (slave).Perpendicular baseline = 184 m., temporal baseline = 16 days.(B) Black dots: pixel phase values as a function of elevation.(C) Corresponding stratified delay map predicted by ERA-Interim.(D) 2 π -wrapped differential interferogram corrected from the topography-dependent atmospheric effects.

Fig. 7 .
Fig. 7. Flow chart of the procedure employed for generating uncorrected and corrected deformation time series (DTS).

Fig. 8 .
Fig. 8. Correlation between RMS map and topography: (A) amplitude of COSMO-SkyMed SAR image on 25 January 2011 (master).(B) Digital elevation model in radar geometry generated from X-Band SRTM (Shuttle Radar Topography Mission).(C) Average RMS map.(a)-(f) represent HH and VV deformation time series.