Land Subsidence Monitoring with Differential SAR Interferome

Land Subsidence Monitoring with Differential SAR InterferometryTazlo Strozzi, Urs Wegmiiller, Luigl Tosl, Gabrlele Bitelli, and Volker Spreckels
The potential of differential SAR interferometry for land subsidence monitoring is reported on. The principle of the technique and the approach to be used on a specific case are first presented. Then significant results using SAR data from the ERS satellites for various sites in Germany, Mexico, and Italy, representing fast (mlyear) to slow (mmlyear) deformation velocities, are discussed. The SAR interferometric dis~lacement maps are'validated with available leveling data. he accuracy of the subsidence maps produced, the huge SAR data archive starting in 1991, the expected continued availability of SAR data, and the maturity of the required processing techniques lead to the conclusion that differential SAR interferometry is suitable for operational monitoring of land subsidence. Land subsidence (land surface sinking) occurs in many parts of the world, particularly in densely populated regions on deltas (Poland, 1984;Barends et al., 1995). Subsidence is the result of the natural compaction of sediments; of extraction of ground water, geothermal fluids, oil, gas, coal, and other solids through mining; and of underground construction. Most of the major subsidence areas around the world have developed in the past half-century at accelerated rates due to the rapidly increasing use of ground water, oil, and gas. Even if the hazards associated with subsidence are different from those caused by sudden and catastrophic natural events like floods and earthquakes, because surface sinking is a slow event, expansive damage can occur. In particular, many areas of known subsidence are along coasts where the phenomenon becomes obvious when the ocean or lake waters start coming further up on the shore. In addition, spatially heterogeneous subsidence produces damage in buildings and in other man-made structures such as bridges, highways, electric power lines, railroads, and underground pipes. Spatially heterogeneous subsidence is also the cause of changes in the drainage patterns, with canals that no longer carry their original design flows. For all these reasons, there is an increasing demand for accurate monitoring of land subsidence in order to improve the understanding of the phenomena and to give recommendations for a sustainable use of the underground resources.
Abstract
Introduction
T. Strozzi and Urs Wegmiiller are with Gamma Remote Sensing, Thunstrasse 130, 3074 Muri BE, Switzerland (strozzi8 gamma-rs.ch). L. Tosi is with ISDGM, CNR, S. Polo 1364, 30124 Venezia, Italy. G. Bitelli is with DISTART, University of Bologna, Viale Risorgimento 2, 40136 Bologna, Italy. V. Spreckels is with the Institute for Photogrammetry and Engineering Surveys, University of Hannover, c/o Deutsche Steinkohle AG (DSK),DG, Gleiwitzer Platz 3, 46236 Bottrop, Germany.PHOTOGRAMMmRIC ENGINEERING & REMOTE SENSING
The problems of land

subsidence are a known environmental aspect in many cites and were among those included in Internathe list of research projects recommended-~~o's tional Hydrological Decade, which began in 1965, and the International Hydrological Program, a continuing UNESCOprogram beginning in 1975 (Barends et al., 1995). Traditionally, land subsidence is determined with precision leveling surveys. Wide networks of benchmarks were set up in most of the subsiding cities, and measurements are periodically repeated (FigueroaVega, 1976;Poland, 1984;Balestra and Villani, 1991; Brighenti, 1991;Capra et al., 1991; Bertoni et al., 1995; Carbognin et al., 1995; Gottardi et al., 1995).Leveling allows the monitoring of land subsidence in selected locations with high precision but is, for large areas, time consuming and expensive. Recently, GPS stations were implemented in few significant points of large subsiding areas in order to have a fast, cheap, and global measurement of subsidence (Bitelli et al., 2000). In the context of traditional surveying techniques, differential SAR interferometry has the potential to provide very important subsidence information over urban areas because of its two-dimensional spatial coverage, its high vertical accuracy, the SAR data availability, and its competitive cost. The potential of differential synthetic aperture radar (SAR) interferometry to map coherent displacements at centimeter to millimeter resolutions resulted in recent years in spectacular new results for the geophysical sciences. Earthquake displacement (Massonetet al., 1993),volcano deformation (Massonetet al., 1995),glacier dynamics (Goldstein et al., 1993;Joughin, 1996, Mohr et al., 1998),and land subsidence (Strozzi and Wegmiiller, 1999; Strozzi et al., 1999;Fruneau et al., 1999; Ferretti et al., 1999; Strozzi et al., 2000b) were mapped. In the SARinterferometric approach, two SAR images are combined to exploit the phase difference of the signals. The interferometric phase is sensitive to both surface topography and coherent displacement along the look vector occurring between the acquisitions of the interferometric image pair. The basic idea of differential SAR interferometry is to subtract the topography related phase from the interferogram to derive a displacement map. For a specific case, different approaches can be applied depending on the availability of a digital elevation model (DEM), on the characteristics of the SAR data with respect to spatial baseline, on acquisition time difference and coherence, on the displacement rates and shapes, on the land cover, and on the topography. In this contribution we review the principles of the differential interferometric SAR technique for subsidence monitoring, discuss how the approach is adapted to a specific case, and present significant results for sites in Germany, Mexico, and
Photogrammetric Engineering & Remote Sensing Vol. 67, No. 11, November 2001, pp. 1261-1270.0099-1112/01/6711-1261$3.00/0O 2001 American Society

for Photogrammetry
and Remote SensingNovember ZOO/1261
Italy. Because the results are based on SAR data from the European Remote Sensing Satellites ERS-1and ERS-2, emphasis is given to the characteristics of these sensors. Taking into account the limitations of SAR interferometry in non-urban environments, we also propose a strategy for the integration of leveling, GPS, and SAR data in order to achieve an accurate, rational, and cost-effective monitoring.
Differential SAR Interferometry for Land Subsidence MonltorlngSAR Interferometry
Here, R is the distance from the SAR to the scatterer, 0is the incidence angle, and B, is the component of the baseline perpendicular to the line-of-sight direction. The baseline is small be assumed to be the same for both senenough that R and @can sor positions. For the ERS SAR configurations, with a wavelength of 5.66 cm, a nominal incidence angle of 23", and a nominal slant range of 853 krn,Equation 4 reduces to
In the interferometric approach, two SAR images acquired from slightly different orbit configurations and/or at different times are combined to exploit the phase difference of the signals (Zebker et al., 1994;Bamler and Hartl, 1998;Wegmiiller et al., 1998;Mohr and Madsen, 1999).By assuming that the scattering phase is the same in both images, the interferometric phase 4 is a very sensitive measure of the range difference R2 - R1: i.e.,
Here, and are the phases of the first and second SAR images, respectively; R1is the distance from the SARto the scatterer by the first acquisition; R, is the distance by the second acquisition; and A is the wavelength. The interferometric phase is determined as the argument of the normalized interferogram y, defined as the normalized complex correlation coefficient of the complex electromagnetic fields sl and s, backscattered by the illuminated elements at positions R1 and R2:i.e.,
For a perpendicular baseline component of 50 m, the topographic height change associated with a 2 n-cycle is 188 m. On the other hand, an error of 10 m in the estimation of the topographic height (e.g., as a consequence of an inaccurate DEM) results, for a perpendicular baseline component of 50 m, in a phase error of 0.11n-. For interferograms with small perpendicular baseline components, the effect of the topographic related phase is very small but not negligible for extremely slow subsidence rates. Differentiation of the first term of Equation 3 with respect to the slant-range distance y yields the fringe frequency in slant-range for a flat surface (Bamler and Hartl, 1998):i.e.,
Differentiation of the first term of Equation 3 with respect to the azimuth distance x yields the fringe frequency in azimuth for a flat surface
Here, the brackets (s) stand for the ensemble averaging of s, and s* is the complex conjugate of s. Coherent averaging over statistically independent samples reduces the variance of the interferometric phase. The coherence ) A, a measure of the phase noise, is the

magnitude of the normalized interferogram. The coherence is limited to values in the interval from 0 to 1. The interferometric phase is sensitive to both surface topography and coherent movement of the scatterer in the lineof-sight direction between the two observations, with inhomogeneous propagation delay (so-called atmospheric artifacts) and phase noise introducing the main error sources. For the purpose of this analysis, the interferometric phase is linearized as (Mohr and Madsen, 1999)
the change of BII along the line of flight Here, aBlllaxrepresents as a result of non-parallel orbits. The phase term related to the coherent displacement of the scattering center along the radar look vector, rdi,,,may be transformed to a vertical displacement, r,,,, if deformation in vertical direction can be assumed: i.e.,
where BII denotes the component of the baseline, i.e., of the separation of the two antennas in space, parallel with the line-ofsight direction, rdisqis the displacement of the scatterer in the line-of-sight direction, and ram,represents path length changes caused by different atmospheric conditions during the acquisitions of the two SAR images. The interferometric phase 4 is still ambiguous to within integer multiples of 2 n-. The problem of overcoming this ambiguity through phase unwrapping is addressed in the next section. Differentiation of the first term of Equation 3 with respect to the surface height yields the relationship between a change in the topographic height z and the corresponding change in the interferometric phase 4 (Bamler and Hartl, 1998):i.e.,1262 November 2001
Here, coherent means that the same displacement of adjacent scatters is observed. In the case of ERS SAR data, a phase change of 2n-corresponds to a vertical displacement of 3.07 cm. In general, however, surface deformations may also occur in directions other than the vertical. Strictly speaking, a single SAR interferometric observation does not allow the full determination of the magnitude and direction of a surface deformation, but allows only an observation of the deformation component along the radar look direction. Changes in the effective path length between the SAR and the surface elements as a result of changing permittivity of the atmosphere lead to a non-zero atmospheric delay, r,,,. The relevant changes in the atmospheric permittivity conditions are mainly associated with water vapor content. Heterogeneity of the water vapor content in the atmosphere can cause heterogeneity of the propagation delay. Currently, it does not seem to be possible to completely avoid these limitations. We address the minimization of errors introduced by atmospheric effects with the interferogram stacking technique (in a later section of this paper).PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING
Finally, the standard deviation of the phase noise gn0i,,, reached asymptotically for a large number of looks N, is a function of the coherence (Rodriguez et al., 1992):i.e.,
(

4
The coherence can be pided into three dominant contributions (Bamler and Hartl): the influence of the signal-to-noise ratio ythemal,spatial decorrelation y,,,,l related to the imaging geometry (interferometricbaseline, local incidence angle), and temporal decorrelation yiempo,l. The thermal decorrelation is usually not a problem for the E R S SAR. Geometric decorrelation is related to the different look angles of the two SAR acquisitions and leads to a critical baseline over which the interferometric phase is pure noise. In the ERS case, the critical baseline for horizontal terrain is about 1150m. The geometrical decorrelation effect related to surface scattering can be avoided using processing filters tuned to different center frequencies, a procedure known as spectral shift filtering, at the cost of a reduced range resolution (Gatelli et al., 1994).However, effects of the terrain slope and of volume scattering are usually not considered in the spectral shift filtering and reduce the range of baselines useful for subsidence studies to values smaller than a few hundred of meters. Temporal decorrelation is the result of the changed geometrical configuration of the scatterers between the acquisitions of the interferometric pair (Strozzi et al., 2000a). For example, water and forest show very low coherence even for acquisition time intervals of one day. Agricultural fields show intermediate coherence, still useful for phase interpretation after one or two months in temperate regions during the winter. Urban and arid areas, finally, show high coherence also for acquisition time intervals larger than one year. The main problems of high phase noise are the lack of spatial information for certain land-use classes and the difficulties in the phase unwrapping. The basic idea of differential SAR interferometry is to separate the topography and displacement related phase terms to allow the mapping of displacement (Wegmiiller et al., 1998).In the socalled two-pass differential interferometry approach, the topography-related phase term is calculated from a DEM. Alternatively, in the so-called three-pass and four-pass approaches, the topographic-relatedphase term is estimated from an independent interferometric pair with a short acquisition time interval, such as an ERS-112 Tandem pair (one-day acquisition time interval). It is referred to as the three-pass approach if one SAR image is common to both interferometric pairs. In many cases, the use of a DEM turns out to be more robust and operationally practicable. The phase unwrapping required in the multi-pass approach is often difficult to resolve and far from operational for low coherence areas, especially in rugged terrain. In addition, gaps in the unwrapped topographic phase for areas of too low coherence may be present, depending on the phase unwrapping method used. Previous to subtraction from the complex interferogram, the topography-related phase term has to be scaled according to the perpendicular

baseline component (Equation 4). The exact knowledge of the baseline is therefore very important. At present, we use various estimation methods based on the orbit data, the image-to-imageregistration offsets, and the fringe rate of the interferograms. For two-pass differential SAR interferometry, but also in order to combine the information retrieved by the SARsystem with other information in map coordinates, geocoding (i.e., the transformation between the range-Doppler coordinates of the SAR and orthonormal map coordinates) is required. For automatic terrain corrected SAR geocoding with a DEM in orthonormal coordinates, we use the approach presented by Wegmiiller (1999).DlfferentlalSAR Interferometry
Once the differential interferogram is computed, the phase has to be unwrapped. From the interferogram, the interferometric ~ h a s is e onlv known to Modulo 2 nand the correct multiple of &has to be added. Phase unwrapping is a problematic step due to discontinuities and inconsistencies as a result of high phase noise. In a first step, filtering and multi-looking are used to reduce the phase noise (Goldstein and Werner, 1997). Then, phase unwrapping is done by applying a region-growing algorithm to the filtered interferogram (Rosen et al., 1994).Critical areas with very low coherence and inconsistencies are identified and avoided in the phase unwrapping. In a final step, the unwrapped phase is converted to a displacement value. A geocoded map of displacement values per pixel is the final product of differential SARinterferometry. Phase distortions caused by spatially heterogeneous atmospheric (in particular, water vapor) and ionospheric conditions limit the accuracy of differential SAR interferometry. Large atmospheric distortions can often be identified by their specific shape, by cross-comparison of interferograms, or possibly based on meteorological data. It is not clear, though, how to best integrate such tests in an operational processing chain. The goal of our analysis is to achieve an operational and robust technique. For this reason, we take a conservative assumption on the atmospheric distortions in single interferograms. As a typical atmospheric distortion, we consider a relative high phase error of n12, which corresponds for the ERS SAR configurations to an error of approximately 0.7 cm in the estimation of the displacement. In order to obtain reliable displacement estimates, the displacement signal should dominate over the error term. In order to keep the expected error on the order of 5 percent, the displacement phase term has to be 20 times the assumed atmospheric distortion, i.e., l o r or 5 fringes, corresponding to around 14 cm of displacement for the ERS SAR configurations. Because the displacement increases with increasing acquisition time interval, interferometric pairs with long acquisition time intervals are preferred to reduce the effect of the atmospheric distortions. For instance, we preferably select an interferogram wi

th three or more years acquisition time interval to map the displacement of an area with velocities up to 5 cmlyear. In order to map faster displacement velocities, shorter intervals are preferred to study the temporal dynamic of the displacement. For slower displacement velocities, even longer intervals would be required. The time interval cannot be extended much beyond a few years because of data availability and because the use of very long time intervals introduces excessive temporal decorrelation that precludes interpretation of the data. A welcome approach to improve the ratio between the displacement signal and the atmospheric phase errors is the stacking of interferograms. Under the assumption of a stationary process (i.e., assuming that the displacement velocity is constant in time), the displacement terms add up linearly, whereas the atmospheric errors (for which we can assume statistical independence for independent interferograms]increase only with the square root of the number of pairs considered. The mathematical framework of the stacking technique is simple. Let us consider n independent interferograms with different acquisition time intervals ti and unwrapped phase cPp The sum of all tjresults in a total acquisition time interval tc, and the sum of all q5jresultsin a cumulative unwrapped phase &,. The average displacement velocity along the look direction v disp is computed aslnterferogram StacMng Technique
and the displacement velocity estimation error asNovember 2001
I
PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING
1263
Here, E is the assumed phase error of a single interferograms (e.g., 7r/2) and A is the wavelength. Equations 10 and 11can be used to determine the cumulative time t,,, required to map a certain displacement velocity vdjspwith a predefined expected relative estimation error. The potential of the stacking technique is demonstrated by the fact that the stacking of more than ten independent interferograms with atmospheric distortions up to ?r/2 allows a cumulative time interval of more than 20 years to be reached and results in an expected displacement velocity estimation error of around 1mmlyear. Such interferogram stacking is made possible by the immense ERS SAR data archive. In a real case, not all interferograms may be statistically independent because a sAR image may be used for various pairs. In this case, the expected relative error is larger than the one computed with Equation 11.In addition, due to different ground conditions and spatial baselines, the coherence of the interferogramsmay be different, resulting in a different coverage with unwrapped-phase information. In our methodology, a displacement value is only determined for points with more than a certain number of valid values. Finally, because interferogram stacking is based on the assumption of stationary displacement velocities, information on the temporal dynamic of the displacement is lost. The successful application of differential SAR in

terferometry for surface displacement mapping depends on factors related to the phenomenon itself (displacement velocitv field, temuoral and spatial gradients), scene parameters (topography,l&d use), sensor ~arameters [data availabilitv. acauisition time interval, basiline), the availability and accuracy of a DEM, and data processing related factors (baseline estimation, phase unwr-apping, geocoding, stacking of interferograms). The assum~tion of a known disulacement direction is usually appropriate in the case of subsLdence caused by ground water extraction where the displacement is nearly vertical. In the case of mining-induced ground movement, on the other hand, the deformation geometry is more complicated and even a combination of SARdata acquired in ascending and descending modes does not allow the resolution of the complete threedimensional displacement field without the use of additional information. Information on the displacement velocity field and, in particular, on its temporal and spatial gradients is, if available, very important for the selection of the most appropriate SAR data with respect to the acquisition time interval. In order to transform the displacement phase term to an absolute displacement velocity, a reference point is required. Often the scene content does permit the specification of stable points. In other cases, leveling surveys or GPS measurements may provide the reference ~ o i n t . ~ c e ntopography i and land use also have an influence on the accuracy and robustness of the interferometric technique. The estimation of the displacement phase term requires the prior estimation of the topographic phase term. Inaccuracies in the DEM directly translate into errors in the displacement measurement (see Equation 4). In addition, in areas of very rugged topography, the successful application of SAR data is strongly reduced by the particular imaging geometry of the sensor. The land use affects the coherence. For acquisition time intervals longer than one year, high coherence is mainly found in urban and suburban areas. Sensor parameters such as the wavelength, the spatial resolution, and the incidence angle have a strong influence on the feasibility of the SAR interferometric displacement mapping.4.L
Error Esthtlons for a S u c c e s s f u I Application of the Technique
Such factors are not discussed here, because the examples shown are based only on SAR data from the ERS satellites. Data availability is another factor that affects the possibility of generating subsidence maps. While for most of the European sites a large number of ERS acquisitions are available, allowing the optimization of the data selection with respect to acquisition dates and interferometric baselines, for some regions in developing countries only a small amount of ERS SAR data was acquired. Data processing related aspects, which influence the robustness and operationality of the application, include baseline estimation, phase unwrapping, geocodin

g, and the stacking of interferograms. Not so much the optimization of the processing for a specific example, but a robust processing scheme that works under most circumstances is the major objective. In our processing chain, we identified phase unwrapping and baseline estimation as the most critical points. The geocoding scheme, on the other hand, turns out to be operational and robust when an external DEM is used (Wegmiiller, 1999). Phase unwrapping is a difficult task because the coherence in long time-delay interferograms tends to be very low for large parts of the image. In some cases, patches of coherence sufficiently high to reliably interpret the interferometric phase are disconnected in space and separated by areas with high phase noise. The unwrapping of sparse data may be successfully achieved if accurate subsidence models are available or if the subsidence geometry is simple and the displacement rate slow (Costantini and Rosen, 1999). Errors in the estimation of the baseline affect the scaling of the topography-related phase and the flattening of the interferogram (i.e., the removal of the fringe frequency corresponding to a flat surface) in range and azimuth directions. We first consider the relationship between a change in the topographic height az and the corresponding change in the interferometric phase a4 (Equation 4). An error SB, in the estimation of the perpendicular baseline affects the scaling of the topographic height, resulting in a phase difference error S(a4). By transEquation 5 for forming S(a4) to a displacement error 8(rdisp), the ERS configuration can be rewritten as
A typical estimation error of the perpendicular baseline of 0.1 m for a height change of 100 m results in a displacement error of 0.03 mm. For relatively flat terrain, the scaling of the interferogram is therefore not a significant problem. Further, we consider baseline estimation errors affecting the flattening of the interferogram in the slant-range direction y (Equation 6). An error SB, in the estimation of the perpendicular baseline results in an error S(a4) in the phase difference in the range direction. In the case of the ERS configuration, this of corresponds to a displacement error 6(rdjsp)
A typical estimation error for the perpendicular baseline of 0.1 m results in a displacement error of 2.8 mrn in a 10-krn distance. Therefore, flattening is a critical step in subsidence measurements of large areas, and the required baseline-estimation accuracy is high. Stacking of interferograms may reduce the errors caused by flattening because the baseline estimation error can be assumed to be statistically independent for the different interferograms. Finally, we consider baseline estimation errors affecting the flattening of the interferogram in the azimuth direction x in the estimation of the change (Equation 7). An error S(aBllIax) of BII along the line of flight results in an error a(&$) in the phasePHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING
TAB

LE 1 . ESTIMATED SUBSIDENCE VELOCITIES FOR THE SELECTED SITES. INALL CASES ERS SAR DATA AREUSED.THEEXPECTED ACCURACY IS COMPUTED ASSUMING AN ATMOSPHERIC DISTORTION OF r / 2 . NOOTHER ERRORS ARECONSIDERED Site Ruhrgebiet (Germany) Mexico City (Mexico) Bologna (Italy) Euganean Geothermal Basin (Italy) Velocities [cmlyear] 0-200 0-40 0-6 Monitoring interval1month 3 months
Number of interferograrns
Cumulative time interval 1 month 6 months 4 years 20 years
Expected accuracy [cmlyear]
1 26 10
0-0.4
1 year 5 years
difference in the azimuth direction. Transforming the phasedifference error into a displacement error S(rdisp) results, for the ERS configuration, in
A typical estimation error for the parallel baseline change of 5 cm per 100 km results in a displacement error of 5 rnm for a 10km distance in azimuth. The flattening of the interferogram in the azimuth direction is again a critical step for the monitoring of large areas. The ERS satellites carry SAR instruments suitable for interferometric data analysis. Their %-day repeat-orbit cycle is appropriate for subsidence monitoring in most of the cases. The large archive of ERS SAR data available since 1991 allows subsidence studies with SAR interferometry in many regions around the world. For this analysis, we selected four sites characterized by
Results of Land Subsidence Maps
different land use, topography, and temporal and spatial dependence of the movement (see Table 1):the Ruhrgebiet (Germany), Mexico City (Mexico),the Euganean Geothermal Basin (Italy),and Bologna (Italy). While a single successful example is sufficient to demonstrate the potential of a technique, the analysis of several cases allows us to assess its maturity and operational readiness. A more detailed description of these studies can be found in Wegrniiller et al. (2000),Strozzi and Wegmiiller (1999), Strozzi et a]. (1999), and Strozzi et al. (2000b). In all cases, SAR processing was performed to the same Doppler centroid in order to maximize the fringe visibility (Werner eta].,2000). All the complex SAR images were then registered with sub-pixel accuracy to the same reference geometry. The interferograms were calculated with five looks in azimuth and one look in range (ten looks in azimuth and two looks in range for Mexico City) using spectral shift filtering. Based on the orbit parameters and the average interferogram fringe frequency, the baseline was estimated. For the two Italian
1
(a)
(b)
(c)
Figure 1. Ruhrgebiet: Deformation maps for an active coal excavation site. The image width i s 4 km. The isocinetics indicate the deformation along the SAR observation direction with intervals of 0.5 cm in 35 days. For the image brightness a backscattering intensity image is used. The time period is (a) 20 September 1995-30 November 1995, (b) 30 November 1995-04 January 1996, (c) 05 September 1996-10 October 1996, (d) 02 February 1997-13 April 1997.
I
PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING
November 2001
1265

本文来源：https://www.bwwdw.com/article/41hl.html