Open access peer-reviewed chapter

Satellite SAR Interferometry for Earth’s Crust Deformation Monitoring and Geological Phenomena Analysis

By Giuseppe Solaro, Pasquale Imperatore and Antonio Pepe

Submitted: November 16th 2015Reviewed: May 16th 2016Published: September 8th 2016

DOI: 10.5772/64250

Downloaded: 465

Abstract

Synthetic aperture radar interferometry (InSAR) and the related processing techniques provide a unique tool for the quantitative measurement of the Earth’s surface deformation associated with certain geophysical processes (such as volcanic eruptions, landslides and earthquakes), thus making possible long-term monitoring of surface deformation and analysis of relevant geodynamic phenomena. This chapter provides an application-oriented perspective on the spaceborne InSAR technology with emphasis on subsequent geophysical investigations. First, the fundamentals of radar interferometry and differential interferometry, as well as error sources, are briefly introduced. Emphasis is then placed on the realistic simulation of the underlying geophysics processes, thus offering an unfolded perspective on both analytical and numerical approaches for modeling deformation sources. Finally, various experimental investigations conducted by acquiring SAR multitemporal observations on areas subject to deformation processes of particular geological interest are presented and discussed.

Keywords

  • deformation modeling
  • geodesy
  • SAR interferometry

1. Introduction

Synthetic aperture radar interferometry (InSAR) is a consolidated technique that can be used to measure crustal deformation (associated with volcanic and seismic activity) by exploiting the phase of coherent electromagnetic signal. Specifically, theoretical foundation of the spaceborne (across-track) SAR interferometry and multitemporal advanced processing (e.g., persistent-scatters and small-baseline based) methods are introduced. First, we methodologically address the InSAR methods allowing the detection, mapping and monitoring of the Earth’s crust dynamic processes (surface displacements) over large temporal and spatial scales with centimeter to millimeter accuracy. Then, emphasis is placed on the geological processes taking place within the Earth's crust, such as the movement of a seismogenic fault, the accumulation of magma, variation of pressure in the magmatic reservoirs, subsidence. All these phenomena can cause deformations of the Earth's surface and can then be investigated by suitably exploiting satellite observations. For such a purpose, different approaches are possible; most of them are based on the inversion of a suitable model describing the underlying geophysical phenomenon. Specifically, in order to model the deformation sources both analytical and numerical approaches have been adopted. Within the analytical framework, we first address the most commonly adopted models, which can reproduce the observed deformations in a sufficiently realistic way by using simple functions characterized by a limited number of parameters. Although these analytical models neglect several aspects (e.g., the properties of magma inside the source, including its compressibility, the asperities along the fault plane, the crustal heterogeneity), they still constitute a valuable tool for a preliminary evaluation on the localization and geometric characteristics of the sources. Numerical modeling, which is a powerful tool allowing a realistic simulation of geophysical processes, using heterogeneous information and efficient computational methods, is also discussed. Specifically, various numerical modeling techniques exist; one of the most used in the Earth Sciences community is the finite element method (FEM) technique. In fact, both the increase in knowledge about geophysical systems and technological development of numerical techniques have enabled the implementation of complex modeling approaches, which are able to represent the spatiotemporal variability of the geophysical parameters of interest. In this context, the use of FEM multiphysics tools represent a new frontier for the understanding of the spatial and temporal evolution of different geodynamic settings, such as volcanic and seismic areas and those with a hydrogeological instability. Therefore, a comprehensive and updated perspective is offered in this chapter, encompassing advanced remote sensing and geophysical methodologies addressed to the analysis of several natural phenomena resulting in the deformation of the Earth’s crust. Furthermore, a wide range of case studies is shown, which have systematically been investigated by considering data acquired by different SAR sensors (e.g., ENVISAT, RADARSAT-2) on diverse hazardous geologically zones of interest (e.g., areas interested by seismic and volcanic activity).

2. SAR interferometry principles

Synthetic aperture radar (SAR) [13] is a coherent active microwave remote sensing system widely used for the Earth remote sensing. SAR instruments can be mounted on-board aircraft or satellite platforms; they work by transmitting microwave pulses toward the Earth surface and by measuring the microwave echoes scattered back to the sensor platform. SAR is an imaging system with all-weather, day and night sensing capability that nowadays plays a key role for the remote sensing of the environment, and in particular it is extensively used for the monitoring and analysis of several geophysical phenomena. A SAR image can be represented as a two-dimensional (2D) complex signal in the (range, azimuth) plane, whose amplitude gives information about the backscattering coefficient of the ground and the phase includes information about the distance traveled by the emitted electromagnetic pulses from the transmitting to the receiving antennas (i.e., twice the sensor-to-target distance). Range (or cross-track) direction is associated with the “line-of-sight” distance from the radar to the target, whereas azimuth (along-track) direction is parallel to the flight track.

One of the major applications of the SAR technology is represented by the SAR interferometry (InSAR) technique [48], which is based on the measurements of the phase pattern difference between two complex-valued SAR images acquired from two different orbital positions, and allows the measurements of geomorphological characteristics of the ground, such as the topography height and its modifications over time (e.g., the surface deformation) due to earthquakes, volcano eruptions, or other geophysical phenomena. Historically, the main application of InSAR was the retrieval of the terrain topography [46]. Depending on the time when SAR acquisitions are collected and the orbital position of the SAR platform, different InSAR configurations can be distinguished. Cross-track interferometry is a basic SAR interferometric configuration in which two antennas are arranged across the track of the platform, as sketched in Figure 1.

Figure 1.

SAR interferometric configuration. The black lines show radar signal paths for an interferogram pair formed by the antennas M and S.

Within this context, two different acquisition modes can be distinguished: single-pass mode is characterized by two distinct antennas on the same platform (in the standard form, the former (master) operating in a receive/transmit mode and the latter (slave) in the receive mode only), the repeat-pass mode concerns two separate passes of a single SAR mission over the same area [8]. In addition to the standard cross-track interferometric configuration, we also mention the along-track interferometry (ATI), which is a single-pass configuration with two antennas displaced with a baseline parallel to the direction of motion: airborne ATI has been mainly used for measurement of ocean currents.

Let us consider again the imaging geometry depicted in Figure 1, where the first SAR image (i.e., the master image) is taken from the orbital position labeled to as M, and the second one (i.e., the slave image) is captured from the orbital position labeled to as S, at a distance b (typically referred to as baseline) from M. Taking into account simple geometrical considerations relevant to the considered geometry, it is possible to uniquely locate each imaged targets on the ground and get an estimate of their heights (namely, z) above the reference plane. As evident by inspection of Figure 1, if a same target (namely, T) is observed from two orbital positions (master and slave), the difference between the path lengths to the target can be correctly measured and the target height z above the assumed zero-altitude plane can be unambiguously determined. This is obtained by taking into account the following two equations (see Figure 1):

(r+δr)2=r2+b22rbsin(ϑα)E1
z=HrcosϑE2

where δr and δ + δr represent the radar ranges from the corresponding antennas to the target point being observed, ϑ is the radar look angle, α represents the angle of the baseline relative to the horizontal, z denotes the scatterer height above the flat-earth reference, H is the height of the sensor above the reference surface, and b is the physical separation of the antennas that is referred to as the baseline of the interferometer. Notice that (1) derives from the application of the cosine rule to the MST triangle and (2) is a simple geometric relationship linking the target topography (z), the sensor height (H), and the radar side-looking angle (ϑ). The ability in successfully reconstructing the unknown topography (z) is strictly dependent on the capability to precisely measure the slant-range difference δr, which represents one of the known terms of the system of Eqs. (1) and (2).

Historically, a first methodology to get an estimate of δr was represented by the radar stereometry [8]. In such a method, the master/slave sensor-to-target slant-range difference δr is measured by searching for the position of the same target in the two coregistered SAR images (being the coregistration the operation needed to spatially aligned one SAR image to another) [9, 10]. As a matter of fact, the attainable accuracy in estimating δr is on the order of the system slant-range resolution. However, it can be proved that the errors in the estimation of δr is magnified by a factor on the order of the ratio (rb)when they are transferred to height measurements [3], thus leading to an inaccurate measurement of the target height (z). For instance, we consider ENVISAT platform parameters (rb=800km100m)and suppose being able to discriminate reasonable range displacements of 1/16th of the pixel spacing through use of correlation digital processing (i.e., the accuracy in measurement of δr is equal to 0.5 m). Accordingly, the achievable height accuracy turns out to be on the order of kilometers, and it is evidently unacceptable. This is the main reason of InSAR success with respect to radar-stereometry. Indeed, the intrinsic limitation of radar stereometry due to the low attainable accuracy of topography is fully overcome by SAR interferometry, which allows estimates of the master/slave slant-range difference δr with centimeter accuracy over region of hundreds of kilometers in size at a resolution of a few meters.

In the following, we primarily refer to the repeat-pass cross-track SAR interferometry configuration. Let us consider again the imaging geometry depicted in Figure 1 and assume the radar system has infinite bandwidth and hence with point-wise image pixels [4]; under this condition the master and slave complex-valued SAR images (pixel-by-pixel) can be mathematically represented as follows:

γ^1=γ1exp[j4πλr]E3
γ^2=γ2exp[j4πλ(r+δr)]E4

where γ1 and γ2 are the complex reflectivity functions of the master and slave scene, respectively, and λ denotes the operative radar wavelength. It is worth mentioning that the phase of each single-channel radar signal is composed of two parts: the first represents the propagation phase that depends on the radar-scene distance, the second depends on the inherent electromagnetic scattering process. The interferometric phase map (so called interferogram) is formed on a pixel-by-pixel basis starting from two coregistered (complex) SAR images as follows. For each pixel, the phase difference between the two SAR images is extracted by simply multiplying the first image (master) by the complex conjugate of the second image (slave) and then by extracting its phase term.

From (3), we get the radar observable (interferometric phase):

ψ˜=arg[γ^1γ^2*]=arg[γ1γ2*exp(j4πλδr)]E5

where the asterisk denotes the complex conjugate operation, and the symbol arg[·] refers to the phase extraction operation (i.e., the operator that extracts the phase of a complex number restricted to the ]− π, π] interval). Assuming that the scattering mechanism on the ground has not significantly changed (arg[γ1] = arg[γ2]) between the two passages of the sensor over the illuminated area (mutually coherent observations), the measured interferometric phase ψ˜depends upon purely geometric information on the path difference δr only:

ψ˜=arg[exp(j4πλδr)]E6

The observed interferometic phase ψ˜is 2π-ambiguous, and the obtained image is called an interferogram; the pattern formed by the iso-phase contours is commonly referred to as fringe pattern. Since the ambiguity of the phase measured modulo 2π, the information on range difference δr is then retrieved from the interferogram by applying the phase unwrapping operation [11, 12], thus estimating the inherent absolute interferometric phase ψ, which is given by:

ψ=4πλδrE7

Note also that: ψ˜=W(ψ), where W is the so called wrapping operator [13].

The difference in range from the scatterer to the two aperture phase centers is well approximated (since br, the commonly referred to as parallel-ray assumption is reasonable) as δr = −b sin(ϑα), where b|| = −b sin(ϑα) is just the projection of the baseline along the line of sight (LOS) (Figure 1). Thus, the interferometric phase is given by:

ψ=4πλbsin(ϑα)E8

It is worth highlighting the height sensitivity of ψ, through the dependence of the actual look angle ϑ, on the altitude z = Hr cos ϑ, where H is the height of the sensor above the reference surface. By considering the standard interferometric configuration depicted in Figure 1, it is possible to relate the computed interferometric phase to the (unknown) height topography [4]. At first order, we obtain:

ψψ0+ψzz=4πλbsin(ϑ0α)4πλbrsinϑ0zE9

where z is the topography height above the flat earth reference, ϑ0 is the look angle to the point target assuming zero local height, b = b cos(ϑ0α) represents the projection of the baseline normal to the line of sight from the radar to the target and it is an important parameter referred to as orthogonal baseline. The first term in (9), ψ0=4πλbsin(ϑ0α), accounts for phase contribution generated by an ideally flat-earth (z = 0); this term is present even in the absence of any height elevation above the reference surface. Indeed, across the image swath there will be an equivalent flat-earth variation in phase resulting from the corresponding change of incidence angle from near to far swath edge. In order to avoid that the result be biased with position across the swath, the flat earth variation needs to be removed from the recorded phase, thus removing (interferogram flattening) the high-frequency modulation induced by the “flat earth” phase variations to facilitate further processing. The second term in (9), Δψ=ψzz, is the resulting “flattened” phase difference, with the height sensitivity of the interferometer given by ψz=4πλbrsinϑ0. From (9), it is clear that the sensitivity of the interferometer could be improved by increasing the baseline. The perpendicular baseline, however, cannot exceed the limiting case (critical baseline) for which the variation in the interferometric phase difference across a single ground range resolution element is 2π. Indeed, the arising decorrelation phenomena lead to significant noise disturbances in the computed interferogram [14], hence fraction of the critical baseline are typically used in practice. As a result, a compromise is needed for the selection of the optimal baseline: on the one hand, large interferometric baselines would guarantee more accurate estimates of height topography, on the other hand, large baseline interferograms are more affected by decorrelation noise.

2.1. Detecting surface deformation

Figure 2.

Differential SAR interferometry geometry. Note that r2−r1=(r2−r˜)+(r˜−r1)≅ΔdLOS+δr, where δr=r˜−r1 is the path difference in the absence of any ground displacement, and the LOS displacement, ΔdLOS, is given by ΔdLOS = Δd sin(ϑ – αD), with Δd representing the amplitude of the displacement from P1 to P2.

In this section, we shortly review the basic principles of differential SAR interferometry. Indeed, satellite SAR interferometry nowadays is mostly used for the detection/monitoring of surface changes occurring between the two passages of the radar sensor over the same scene. In such a case, as a slightly change across the two SAR acquisition times occurs in the imaged scene (due to, for instance, subsidence, landslide, or earthquake phenomena), an additive term associated with the radar line of sight (LOS) component of the surface displacement arises in the interferometric phase, in addition to the phase dependence on topography. By the inspection of the imaging geometry depicted in Figure 2, at the first-order we get:

Δψ=ψzΔz+ψdLOSΔdLOS=4πλbrsinϑΔz+4πλΔdLOSE10

where ΔdLOS represents the projection of the surface-displacement vector onto LOS (range) direction, ϑ is the look angle to the point target with respect to the nominal local height, and Δz denotes the residual topographic variation. Note that it is reasonable to assume that the radar echoes remain correlated since the surface displacements are assumed small with respect to a resolution cell. It is also important to note that a much more sensitive dependence of phase (10) results from surface displacement ΔdLOS than from residual topographic variation Δz, insofar as the distance r typically is very much greater than the orthogonal baseline distance b. Notice that, in order to isolate (measure) the interferometric phase term associated with the displacement, it is necessary to remove the interferometric phase contribution pertinent to the underlying topography in Eq. (10). Specifically, the so-called differential SAR interferometry (DInSAR) basically consists in the synthesis of a simulated topographic phase screen from an available digital elevation model (DEM) of the area (using the so so-called back-geocoding process) and to subtract on pixel basis these synthetic fringes leaving only the terms associated with the displacement (see Eq. (10)) [4].

In this ideal configuration, the DInSAR technique gets an unambiguously measurement of the LOS displacement of the order of fractions of wavelength: note that a differential phase change of 2π is converted to a LOS displacement of λ/2. As an example, since the error on the estimate is of a fraction of π and the wavelength is of the order of centimeters (e.g., for the ERS-1/2 case λ = 5.6 cm), we could measure LOS displacement down to millimeter accuracy, provided that coherence of the differential interferograms is sufficiently high. Computed differential SAR interferograms however contain, in addition to the deformation component, some (unwanted) phase terms arising from unavoidably inaccuracies in the knowledge of the actual topographic pattern and/or of the orbital parameters. In particular, the variation of the interferometric phase can be expressed more in general in the form:

Δψ=Δψdisp+Δψtopo+Δψorb+Δψprop+ΔψnoiseE11

where:

  • Δψdisp=4πλΔdLOsaccounts for a possible displacement of the scatterer between observations, where ΔdLOS denotes the projection of the relevant displacement vector on the line of sight;

  • Δψtopo=4πλbrsinϑΔzrepresents the residual-topography induced phase due to a nonperfect knowledge of the actual height profile (i.e., the DEM errors Δz);

  • Δψorb accounts for residual fringes due the use of inaccurate orbital information in the synthesis of the topographic phase;

  • Δψprop denotes the phase components due to the variation of propagation conditions (pertinent to the change in the atmospheric and ionospheric dielectric constant) between the two master/slave acquisitions;

  • Δψnoise accounts for decorrelation phenomena: spatial baseline decorrelation, Doppler centroid decorrelation, thermal decorrelation, and temporal decorrelation (including any change in scattering behavior) [14].

As a final remark, we observe that another source of misinterpretation upon the measured deformation is intrinsic to the InSAR technique itself, and it is due to phase unwrapping errors. Evidently, phase unwrapping errors are integer multiples of 2π but they can propagate within the inversion process, thus significantly affecting the deformation measurements [3].

Figure 3.

Geometric scheme for the deformation components.

Availability of InSAR results computed from SAR data obtained from ascending and descending orbits allows the retrieval of the east-west (E-W) and the up-down (U-D) components of the detected deformation [15, 16]. Let us assume the target “observed” from both the ascending and the descending satellite passes, and assume the displacement components along the ascending and descending radar LOS directions have been estimated. For the sake of simplicity, the following assumptions are made: (i) ascending and descending radar LOS directions (dLOS(asc)anddLOS(desc), respectively) lay on the plane identified by east and –z directions, and (ii) the sensor look angle ϑ is approximately the same for both the ascending and descending observations. In particular, for all the pixels that are common to both radar geometries, the sum and the difference of LOS-projected deformations computed (over approximately the same time period) for the ascending and the descending orbits can be calculated. Based on simple geometric considerations, the E-W and up-down components of the measured surface deformation can be estimated as follows:

dLOS(East)=dLOS(desc)dLOS(asc)2sinϑE12
dLOS(Up)=dLOS(desc)+dLOS(asc)2cosϑE13

Notice that, because of the namely polar sensor orbit direction, the north-south (N-S) component of the deformation cannot be reliably singled out. Geometric scheme to interpret the deformation component is portrayed in Figure 3.

Finally, we emphasize that a fundamental advantage of InSAR technology, with respect to global positioning system (GPS) networks, resides in its dense spatial sampling of the displacement field.

3. Multichannel SAR interferometry

Differential SAR interferometry methodology has first been applied to investigate single deformation events. At the present days, however, it is chiefly applied for the computation of displacement time-series through the so-called multitemporal (or multichannel) interferometric SAR approaches [1725]. These advanced methods are based on the processing of sequences of multitemporal interferograms relevant to an area of interest and are aimed at recovering the expected LOS-projected time-series of deformation. A short overview of the main algorithms proposed up to now is here reported. Generally speaking, multichannel interferometric techniques can be categorized into two broad families, those focused on analyzing persistent scatterers, that is to say point-like targets that are not significantly affected by decorrelation effects [1719], and the small baseline (SB) [2026] methodologies, relying on the investigation of deformation signals related to distributed scatterers (DS) on the ground, which can be however severely corrupted by decorrelation effects. In this latter case, an a priori selection of the exploited SAR data pairs with small baseline values is required to reduce the noise level in the generated interferograms. Despite of their intrinsic differences, both the PS and SB algorithms have successfully been used to detect and monitor deformation phenomena, due to several natural and anthropic hazards, such as volcanic events, earthquakes, landslides, damages to man-made infrastructures in urbanized areas caused by underground, and tunneling excavations and/or gas and water exploitation [2737]. Very recently, a plethora of different PS- and SB-oriented approaches have been implemented and public InSAR toolboxes [1726] are available to users. Recently, some innovative approaches based on the joint exploitation of spatial and temporal relationships among sequences of interferograms and of the statistical characteristics of SAR images involved in their formation have been proposed for the analysis of deformations affecting DS targets [3842]. In particular, the method proposed in [40], which is known in literature to as SqueeSAR, is aimed at retrieving the displacement time-series of DS that are identified by preliminarily studying the statistical homogeneity of adjacent pixels in long sequences of amplitude SAR images, and then by averaging the interferometric phase only on the set of statistically homogeneous (SH) pixels [40, 41]. In addition, the average interferometric phases (associated with couples of images) are jointly employed (for each pixel of the SAR scene) to obtain estimates of the phase relevant to SAR acquisitions [40], thus finally retrieving (for each SH target) a time-series of deformation. This method allows increasing the number of detectable DS targets, but at the expenses of ad-hoc processing for the generation of average (multilook) InSAR interferograms. At variance with the SqueeSAR and other recently proposed multitemporal filtering techniques (e.g., [41, 42]), the method proposed in [41] (and also detailed in [42]) used conventional multilook interferograms, which can be generated by using any of existing InSAR toolboxes and without any preselection of SH targets. This leads to the nonapplicability of statistical framework adopted in [40], which is based on the distributed scattering hypothesis under which the probability density function (pdf) of the complex-valued SAR image may be regarded as being a zero-mean multivariate circular normal distribution. This issue is not considered a very limiting factor in [43], where “conventional” multilook interferograms (also potentially prefiltered using other space-based noise filtering techniques [44]) are filtered in time with the aim to isolate and discard the noise components that are not conservative in time from generated time-series of deformation. The mathematical framework of this new improved SBAS-oriented processing chain is illustrated in [43, 45] where the method is fully detailed. In the following, we focus on the small baseline subset (SBAS) algorithm, originally proposed in [20], by analyzing the underlying basic principles. Let us consider a set of Q single-look-complex (SLC) SAR data acquired over a certain area of interest. One of them is assumed as the reference (master) image, with respect to which the images are properly coregistered. This set is characterized by the corresponding acquisition times {t1,…, tQ} and perpendicular baselines {b⊥1,…, bQ} evaluated with respect to the reference image. Application of the standard SBAS technique starts with the generation of a number, namely M, of small baseline multilook (differential) interferograms. The multichannel phase unwrapping (MCh-PhU) problem consists in the jointly retrieval of the original (unwrapped) phase signals from the modulo-2π measured (wrapped) phases relevant to the considered stack of interferograms. The MCh-PhU operation can be straightforwardly implemented through various 2D [4648] and 3D approaches [44, 4951] (and/or hybrid ones [13, 52]). The variation of the interferometric phase pertinent to the kth SAR data pair can be expressed as (see also (10)):

Δψk=4πλΔdLOSk4πλbkrsinϑkΔz+Δψorbk+Δψpropk+ΔψnoisekE14

where k ∈ {1,…, M} specifies the considered interferometric pair (master/slave) of the multiple baseline configuration used for the generation of the relevant interferogram. Readers are referred to [20] for further details. Once the phase associated to each SAR acquisition, as well as the residual topography, are estimated, the phases are converted to deformation and the atmospheric phase screen (APS) is computed and filtered out from the obtained deformation time-series. APS removal is achieved by exploiting the assumption that APS is spatially correlated and uncorrelated in time, thus processing atmospheric corrupted time-series is performed with a spatial low-pass (LP) filter and a time high-pass (HP) filter [17, 20]. The quality of retrieved LOS time-series is finally evaluated pixel-by-pixel by calculating the values of the temporal coherence factor, defined in [52]. Residual orbital fringes are also estimated and filtered out in the conventional SBAS processing chain by searching for (in each interferogram) any possible phase ramp, which can be directly related to errors in the knowledge of sensor position along its orbit. Such residual phase ramps (see also [38, 39]) are jointly analyzed to correct orbits state vectors. Finally, for pixels with high temporal coherence the map of LOS mean deformation rate over the analyzed time-periods is computed. Note that, whenever ascending/descending SAR data-tracks are available, SBAS processing can be applied for the two complementary orbits. Thus, the ascending/descending rates of deformation can be composed, as described in the previous section, to retrieve the east-west and up-down displacement rates over the time-period span by the available SAR scenes. Finally, we highlight that a parallel computational model for SBAS algorithm is discussed in [13, 26].

4. Geological models and applications

In this section, we describe the technical aspects related on how to retrieve the characteristics of a deformation source from satellite InSAR data, focusing the attention on the seismic, volcanic, and landslide activities. We present the state-of-the-art of the techniques concerning this problem, describing the most commonly used analytical and numerical models, and also providing appropriate geological examples for each kind of modeling approach.

4.1. Analytic modeling

The increasingly widespread use of space geodesy has resulted in numerous, high-quality surface deformation data sets. For example, a dense array of more than 1000 continuous GPS (global positioning system) stations are distributed throughout Japan [53] and more than 700 GPS stations are operating in the California area [http://earthquake.usgs.gov/monitoring/deformation/]. Many geologically active areas, such as Hawaii, Mt. Etna, Campi Flegrei, and Long Valley caldera, have regional GPS networks as well [5558]. At the same time, DInSAR is a well-established technique for studying and analyzing tectonically active areas characterized by wide spatial extent of the inherent deformation [5]. These geodetic data can provide important constraints on fault geometry, its slip distribution and also type and position of a magmatic source. For this reason, over last years, many researchers have developed robust and semiautomatic methods for inverting suitable models to infer the source type and geometry from surface deformation [54]. Most of these methods use elasticity theory and a trial-and-error approach to find geologically plausible deformation models that fit the major features of the observed deformation field [55]. Other authors have systematically searched through a large set of feasible models, comparing the model predictions to the data, and choosing the model that minimizes the misfit [56].

The knowledge of source geometry from geodetic data primarily requires a forward model describing how the crust responds to various kinds of deformation sources. The most commonly used crustal model is the homogeneous, isotropic, linear and elastic half-space [57]. In spite of its limitations, the elastic half-space model is widely used to simulate surface deformation, primarily due to its mathematical simplicity. Sources models commonly used in many papers are [58, 59]:

Figure 4.

Four types of buried point dislocation sources: tensile, dilatational, strike slip, and dip slip (from [58]).

  • The elastic dislocation of a finite rectangular source (Okada model): it is one of the most used model to simulate the surface displacement due to an earthquake, represented as a shear dislocation over a finite rectangular fault [60]. Moreover, the Okada model can also be used to describe magma intrusion like sills or dykes, interseismic and postseismic displacement. Pertinent source parameters are east and north position, depth, length, width, strike angle, dip angle, dislocation (or slip), dislocation angle (rake), opening (for magmatic intrusion) (Figure 4).

  • The point pressure source (Mogi model): it is one of the simplest and effective sources used in volcanology, as its description requires only four parameters: depth, east and north position, volume/pressure variation [61] (Figure 5).

  • The finite spherical pressure source relies on the assumption that the radius of the source cannot be separated from the pressure change. This means that we can only obtain estimates of the depth, location, and power of the source [58]. In [62], the mathematical expressions to approximate the deformation due to a pressurized finite spherical cavity were derived by applying higher order corrections for stresses reflected back on the source by its image.

  • The closed pipe: a model for a plugged conduit or a cigar-shaped magma chamber. It includes a conduit to transport magma from the chamber to the surface. During quiescence period, the magma tends to cool and forms a plug, and the pressure in the magmatic system can increase [58]. The distribution of surface deformation from inflation of a closed pipe is quite different from that previously described for a sphere and this is mainly related to two main aspects: (1) most conduits are quite small relative to magma chambers, and (2) the near-field deformation from an elongate embedded source depends on the value of Poisson’s ratio [63].

  • The open pipe: a composite model for the filling of an open conduit. In [64], it is presented a dislocation model for surface deformation from magma rising in a conduit that is open at the top. A constant cylindrical dislocation is used to model portions of the conduit subject to a uniform pressure change [58].

Figure 5.

Surface deformation from an embedded point pressure source (Mogi model) (from [58]).

In spite of its limitations, the elastic half-space models are widely used to model surface deformation caused by uniform rectangular dislocations [60] and point sources [61]. Moreover, until recently, most geodetic data were not of sufficiently high quality to justify more complex crustal models.

For almost all the listed models, the geometric parameters (position, depth, dimension, orientation, etc.) are nonlinearly related to the surface displacement. On the contrary, other parameters, as the dislocation for the Okada model or the pressure change for a Mogi model, have a linear dependency with the surface displacement [59]. The estimation of nonlinear and linear parameters from geodetic data follows different inversion strategies, which are explained in the next section.

4.1.1. Inversion strategies for source parameters estimation

The relationship between the measured deformation field (which for instance can be inferred through InSAR technique, as discussed in Sections 2 and 3) and the source geometry can be expressed by the following equation:

d=G(m)+εE15

where d is the deformation data vector, m is the source parameter vector (e.g., for a fault, length, width, depth, dip, strike, location, slip are the source parameters to be estimated), and G describes the specific functional form. The ε term is a vector of observation errors. For the source geometry estimation problem the data, in general, are nonlinearly related to the source parameters. For this reason, source estimation reduces to nonlinear optimization [54]. Therefore, we systematically search the finite dimensional parameter space for m, using G to predict the deformation field for a given m. The geodetic signal contains unmodeled deformation such as those arising from elastic heterogeneity or anisotropy, which may contribute to the misfit, thus our best estimated source model is always conditional on the assumptions intrinsic to the forward model.

Derivative-based algorithms, Levenberg-Marquardt or the method of conjugate gradients, offer straightforward and efficient strategies for solving the mentioned optimization problem [54]. These algorithms depend on the gradient and higher-order derivatives to guide them through misfit space; however, due to the nonlinear nature of the G functional form, they can get trapped in the first local minimum that they encounter and never find or even approach the global minimum. Consequently, these algorithms work well only when the initial guess is near the global minimum. A priori information can often provide a good initial guess. Clearly, whether a derivative-based method reaches the global minimum depends on where it starts. Moreover, in [54], it was found that particularly in the case of low measured displacement, the misfit space often contains numerous local minima and lacks a deep, well-defined global minimum. Therefore, derivative-based methods offer a practical approach for retrieving the solution to the geodetic inversion problem only in cases characterized by high measured displacement and good geologic insights, such as the type and location of the deformation source [54].

In spite of their inefficiency, exhaustive and random searches do not have the limitation to remain trapped in a local minimum. In the past, mathematicians have sought algorithms that combined the efficiency of a derivative-based method with the robustness of a random search. The result was the Monte Carlo class of algorithms. The common feature that all algorithms of this class share is an element of randomness that permits an occasional uphill move, that is, the algorithms will not always move from a candidate model with higher misfit to a model with lower misfit [54]. The most common Monte Carlo algorithms are the simulated annealing [65] and the random cost algorithm [66]. Another class of Monte Carlo algorithm includes the genetic algorithms [67].

Simulated annealing. In such a kind of algorithm, the possibility to choose a higher misfit model compared to a lower one mainly depends on the state of the annealing process at the time of the choice [54]. The algorithm gives an estimate of this state dependence in terms of a global time-varying parameter called temperature. At high temperatures, all source models have roughly equal chances of getting picked, whereas at low temperatures the algorithm favors low misfit models. The specific annealing algorithm adopted here follows from the work by Yu and Rundle [65] and Berg [68]. It is called the “heat bath” algorithm and proceeds as follows. The initialization procedure consists of two steps: (1) set bounds on the values for all the model parameters (these bounds can come from geologic constraints or physical limitations) and (2) randomly pick an initial starting model. Cycle through the individual model parameters. The most significant complication to the simulated annealing algorithm is the cooling schedule, i.e., how the temperature changes as the annealing progresses. This plays a crucial role in the successor failure of the optimization. In [69], a critical temperature at which the bulk of the annealing should, for maximum efficiency, occur was defined. In brief, at the critical temperature the system remains cool enough to favor low misfits but still high enough to escape local minima.

Random cost. This algorithm is an alternative Monte Carlo approach for nonlinear optimization problems characterized by many local minima in a broad misfit space [66]. It considers a stochastic process to enforce a random walk in misfit space, which allows it to overcome the increase of misfit and to find the global minimum. In [54], the authors indicate that this algorithm is significantly less efficient than simulated annealing, but it is much easier to implement because it does not require a specific cooling schedule. The random cost approach begins by generating a set of trial models that span a region about an arbitrary a priori model [54].

4.1.2. Geological applications

In this section, we present two examples of deformation sources in volcanic (Lazufre, Chile) and seismic (2012, Emilia earthquake, Italy) environment, by applying the analytic modeling. In the first case, the simulated-annealing-based approach is adopted, while in the second case we apply the Levemberg-Marquardt algorithm (see Section 4.1.1).

4.1.3. Sill and finite spherical sources: the case of Lazufre (Chile) volcano

The Lazufre volcanic area is located on the Chilean-Argentinean border at ~300 km east of the subduction trench (Figure 6). The area contains several morphologically distinct volcanic centers [71, 72]. Only one of these, the Lastarria volcano (~5700 m asl), shows strong and persistent fumarolic activity localized on the recent crater borders and on the western flank (Figure 6).

Figure 6.

Deformation at the Lazufre volcanic area: (a) location of Lazufre; (b) InSAR observation for the period June 1995‒December 2006 acquired by ERS; (c) InSAR observation for the period April 2003‒January 2008 acquired by ENVISAT; (d) details of Lastarria volcano; (e) NNW-SSE profiles across the deformation areas for the ERS dataset (black) and for the ENVISAT dataset (gray); (f) photograph of the Lastarria volcano from the northwest, 10 km distant from the summit [70].

Through InSAR observations, a large-scale elliptical deformation pattern was detected during the period from 1995 to 2008, with a deformation rate ranging from 1.8 to 3.2 cm/year [70]. The observed displacement rate at LAS reaches up to 2 cm/year from 2003 to 2008, with a part of this signal being related to the large-scale deformation field. To retrieve the mean deformation velocity maps of the area the SBAS algorithm (see Section 3) was applied to two SAR datasets acquired by the European Satellite missions ERS-1/2 and the ASAR sensor onboard the ENVISAT satellite, operating both in descending orbits.

Figure 7.

Inversion results of the Lazufre deformation data from 2003 to 2008: (1) observation data, (2) analytic models, and (3) residuals. Lastarria displacement result by simulating a finite spherical source showing (4) the observation data, (5) the analytic model, and (6) residuals highlighting three fumarolic areas (black circles). Dashed lines indicate flank movements (FM), on the western flank of the Lastarria volcano [70].

In order to quantify the sources that are responsible for the observed two-scale deformations [70], the considered analytical models were inverted by applying the simulated annealing method. To isolate the displacement pattern the authors followed two main steps: (1) a linear Pearson correlation coefficient [73] and a search of pixels falling within 95% of a similar trend to the maximum displacement observation point (see CEN in Figure 6) were applied; pixels that are not affected by the deformation were automatically excluded; (2) a subsampling of the cross-correlated dataset using a regularly spaced grid (1 km), thus reducing significantly the computational time without affecting the parameter estimation performance, was applied. Because the observed main deformation pattern is very extended in space and its source is likely laterally extended, in [70] an expanded dislocation plane acting as a sill source model [65] has been assumed, and then its parameters has been estimated. For the sake of simplicity, the models were performed in an elastic half-space medium with a Poisson’s ratio v = 0.25 and a Young’s modulus of E = 50 GPa. Residuals are generally less than 0.2 cm/year with the exception of a near radial-symmetric deformation signal with uplift rates larger than 1 cm/year centered on the Lastarria volcano affecting an area of about 50 km2 (Figure 7).

To further investigate this residual deformation a spherical source model approximation is applied [62]. The residuals are again generally less than 0.2 cm/year, with the exception of the area where the three main fumarolic fields are located, which still shows a discrepancy (i.e., the difference between the satellite observation and retrieved model) up to 0.5 cm/year (Figure 7). The best fitting model suggests a shallow spherical source located between 0.6 and 0.9 km below the Lastarria summit. The source radius is ~0.3 km (between 230 and 360 m) and subject to a volume change of ~13,000 m3/year [70].

4.1.4. Okada fault model: the case of the Emilia (Italy) earthquake

Figure 8.

(a) April 30‒June 17, 2012 RADARSAT-2 InSAR interferogram; b⊥ =447m (perpendicular baseline), ϑ = 30° (look angle); the black square represents the InSAR reference point. Note that the red and blue colors correspond to a sensor-target range decrease and increase, respectively. (b) Analytic modeling of the RSAT-2 displacement map. The blue stars are the locations of the two main shock events and the black rectangles represent the surface projection of the best-fit Okada plane solutions. (c) Residuals map; the blue triangles indicate the locations of the Ml ≥ 5.0 aftershocks occurred after May 20. Table 1 reports the retrieved fault parameters for IFT and MFA (modified from [74]).

On May 20, 2012, a local magnitude (Ml) of 5.9 earthquake occurred near the town of Finale Emilia, in the Central Po alluvial Plain, Italy. The seismic sequence evolved with some decreasing magnitude aftershock events (Ml ≤ 5.1), until May 29, when a Ml = 5.8 seismic event occurred around the Mirandola village, about 10 km SW of the May 20 main shock epicenter (Figure 8). The focal mechanisms for these two seismic events show both a WNW-ESE and E-W oriented nodal planes, respectively, and a ~N-S compressional kinematics [74]. The large amount of data available for the considered area, acquired through InSAR analyses, geophysical and deep borehole geological investigations, allows extensively studying the relationship between the ground deformation fields and the activated fault segments associated with the Ml 5.9 and Ml 5.8 main shocks.

To this aim, an analytic modeling was performed [74] by investigating a RADARSAT-2 (RSAT-2) interferogram (see Section 2) that, encompassing the two main earthquakes, allowed quickly identifying the upper crust regions affected by the faulting processes. In particular, in [74], the authors searched for the faults parameters, by using a nonlinear inversion based on the Levenberg-Marquardt Least-Squares approach [75]; the DInSAR data were subsampled through a QuadTree algorithm [76] over a mesh of about 4600 points. The best fit solution consists of two distinct reverse fault planes, corresponding to the south dipping Inner Ferrara Thrust (IFT) and Mirandola Ferrara Anticline (MFA) for the May 20 and the May 29 events, respectively (Figure 8b and Table 1) (more details are provided in [74]). The model shows a good fit with the measured InSAR data, as clearly highlighted by the residual map in Figure 8c, where values smaller than 2 cm are generally found. However, small areas with higher residuals are also noted; they appear at the locations corresponding to the few aftershocks with Ml ≥ 5.0 (not considered in the inversion procedure) occurred in the same time period covered by the RSAT-2 interferogram.

4.2. Numerical modeling: finite element method

Most of the analytical formulations are based on the assumption of a geologic source (seismic or magmatic) embedded in a homogeneous elastic half-space medium [77]. Analytical elastic models are attractive because of their straightforward formulation. However, active geological areas are usually characterized by severe heterogeneities, nonelastic rheologies and complex topography, which are responsible for significant shallow and surface effects. To meet this need, different numerical procedures can be applied in ground deformation studies to estimate how heterogeneity and topography can affect the deformation field solution.

To make numerical simulations practical, it is necessary to reduce the number of degrees of freedom of the object under study to a finite number. The reduction is called discretization. The product of the discretization process is the discrete model. The most popular numerical techniques in structural mechanics are finite element method and boundary element method (BEM). FEM is the most widely used. The basic concept in the physical FEM is the subdivision of the model into disjoint (nonoverlapping) components of simple geometry called finite elements. The response of each element is expressed in terms of a finite number of degrees of freedom characterized as the value of an unknown function, or functions, at a set of nodal points. The response of the model is then considered to be approximately that obtained by connecting or assembling the collection of all elements. A detailed discussion, which is however beyond the scope of this chapter, can be found in [78, 79].

4.2.1. Geological applications

Two examples of deformation sources in landslide (Ivanchich, Italy) and seismic (2012, Emilia earthquake, Italy) environment obtained by applying the numerical modeling are shown in the next sections.

4.2.2. A steady-state creep model: the case of the Ivancich (Italy) landslide

The Ivancich landslide is located in the southeast part of the historical town of Assisi municipality (Italy) and is affected by an active slow motion. Recurrent damages to buildings and infrastructures caused by the slow landslide evolution led local authorities to carry out geological and geotechnical investigations aimed at implementing effective remedial works and mitigation strategies. The kinematical evolution of the Ivancich unstable mass has been simulated by performing a two-dimensional time-dependent FEM of the active ground deformation [80]. We briefly report here the main results achieved in [80].

Figure 9.

(A) The landslide inventory map of Assisi area; the location of four considered inclinometers is also reported. (B) ERS-ENVISAT mean deformation velocity map with location of the six considered SAR pixels. The thick black line shows the longitudinal cross section A-A’ used for modeling, along which the sectors subdivision is reported. (C) A-A’ 2D section reporting the model geometry of the landslide area with geological units, superimposed on the triangular FE mesh. For further details, see [80].

The longitudinal section along the A-A’ line (Figure 9) has been reconstructed by using the available borehole information, the geomorphological evidences and the inclinometer readings.

In [80], the authors subdivided the slope modeling domain into four geomechanical units: (i) the landslide deposit (unsorted debris), (ii) the upper part of the slope is constituted by the limestone bedrock, (iii) the central part is the pelitic-sandstone bedrock, and (iv) the shear zone, with a thickness lower than 2 m, at a depth ranging between 20 and 60 m. In addition, the analysis of geomorphological evidences and InSAR displacement measurements allowed us to identify four areas showing similar kinematical behavior. InSAR data cover almost 20 years of ERS-1/2 and ENVISAT SAR images acquired between April 1992 and November 2010 and processed through the SBAS technique (see Section 3). Four different subsectors along the landslide shear band, characterized by different creep rate parameters, have been assumed in the mesh domain. The authors chose a deviatoric creep model characterized by a creep rate, depending on the stress state deviatoric component to simulate the behavior of the soil in the shear band [81]. In particular, they proposed that the creep strain rate of the soil in the shear band is the unknown parameter, which can be obtained through an optimization procedure with field data. In Figure 10, a comparison between the time series of six selected SAR pixels and those calculated with the in LOS-projected model is shown. According to the authors, the modeling results highlights that a quasi-linear trend in LOS projection can reasonably describe the variation of the slope displacement over time. Higher displacement rates are calculated for the central portions of the landslide, whereas significantly lower rates are predicted in the upper and lower portions of the slope. Moreover, for the same creep model, they showed the comparison between the time series of the displacement at the top of four inclinometers located along the examined longitudinal section and the model results, and they found a good agreement between field data and model results for all considered inclinometers [80].

Figure 10.

Comparison between the time series of six SAR pixels and the calculated secondary creep model in LOS [from 80].

4.2.3. Discretization of faults model: the case of Emilia (Italy) earthquake

Figure 11.

(a and b) 2D numerical model along the BB’ and CC’ profiles of (a) with the indication of the used boundaries and subdomain settings. The parameters rho, E, and n represent the density, Young’s modulus, and Poisson’s ratio, respectively (see [74] for more details). (c and d) Comparison of RSAT-2 (blue triangles), analytical model (green triangles), and FEM model (red triangles) data evaluated along the BB0 and CC0 profiles, respectively. (e and f) Sections of the displacement maps of the Ml 5.9 and Ml 5.8 seismic events, respectively. The arrows indicate the mean displacement directions. (g and h) Locations of the Okada (green lines) and FEM (red lines) fault solutions superimposed on the numerical model mesh. W1, W2, and W3 as well as Fx and Fy are the widths and active loads along the optimized faults, respectively (see [74] for more details).

A numerical modeling for the retrieved ground deformation of the two Emilia earthquakes, already described in Section 4.1.2, was performed in [74] by using FEM. This modeling approach permits us to take into account geological (rock types) and geophysical information available for the considered area. The two seismic events were analyzed in a structural mechanical context under the plane strain approximation mode, in order to solve for the retrieved displacements [82]. Figure 11a and b reports the geological and structural conditions on which the subdomain setting of the FEM model is based. In [74], a 2D structural geometric domains of the region at depth along the AA’ line (Figure 8a) was derived. A 2D optimization was performed: the two BB’ and CC’ profiles, shown in Figure 8a, cross the areas of maximum deformation associated with the Ml 5.9 and Ml 5.8 seismic events, respectively. The model was made to evolve through two stages: during the first stage (preseismic), the model compacted under the weight of the rock successions (gravity loading) until it reached a stable equilibrium. At the second stage (coseismic), where the stresses were released through a nonuniform slip along the faults, an iterative optimization procedure based on a trial and error approach [82] was used, allowing us to follow the evolution of the faulting processes within the best fit solution retrieval. In [74], the authors applied the following boundary conditions (Figure 11a and b): the upper boundary representing the Earth’s surface was not constrained; the bottom boundary was a fixed constraint; a symmetry condition was assumed for the SSW and NNE areas to make the edge effects as negligible. Moreover, they considered three different boundary settings to simulate the sedimentary and tectonic contacts between different rock successions (Figure 11a and b): (i) free mechanical constrains where the faults are kept locked; (ii) roller constraints, which allow the faults to freely slip under the applied stress field, thus the mechanical discontinuities are considered as active; (iii) boundary loads along which the forces are concentrated and transferred to the boundary subdomains. In Figure 11c and d, a comparison between the best fit solutions for the RSAT-2 data with the analytic and the heterogeneous FEM models along the BB’ and CC’ lines, respectively, is shown. From this analysis, a good fit between the FEM models developed along these profiles and the observed ground deformation pattern is evident, in terms of shape and amplitude of the signal, for both seismic events.

5. Conclusion

This chapter offers an updated and applications-oriented perspective on the satellite InSAR technology, with emphasis on subsequent geophysical investigations. Various phenomena occurring in hazardous geologically zones of interest (e.g., areas interested by earthquake, volcanic activity, or landslide), for which the inherent Earth’s crust deformation pattern can be obtained by suitably processing data acquired by SAR sensors (e.g., ENVISAT, RADARSAT-2), have been investigated. Moreover, the adoption of appropriate geophysical models for the considered scenarios has also permitted to consistently explain the resulting deformation patterns. Finally, the obtained information can be suitably stored in geographic information system (GIS) for the geospatial data management, with important implications in terms of the assessment of geological risks (such as volcanic and seismic), damage assessment, and the proper prevention/planning of human activities.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Giuseppe Solaro, Pasquale Imperatore and Antonio Pepe (September 8th 2016). Satellite SAR Interferometry for Earth’s Crust Deformation Monitoring and Geological Phenomena Analysis, Geospatial Technology, Pasquale Imperatore and Antonio Pepe, IntechOpen, DOI: 10.5772/64250. Available from:

Embed this chapter on your site Copy to clipboard

<iframe src="http://www.intechopen.com/embed/geospatial-technology-environmental-and-social-applications/satellite-sar-interferometry-for-earth-s-crust-deformation-monitoring-and-geological-phenomena-analy" />

Embed this code snippet in the HTML of your website to show this chapter

chapter statistics

465total chapter downloads

2Crossref citations

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Collaborative Uses of Geospatial Technology to Support Climate Change Adaptation in Indigenous Communities of the Circumpolar North

By Megan Sheremata, Leonard J.S. Tsuji and William A. Gough

Related Book

First chapter

Ground Based SAR Interferometry: a Novel Tool for Geoscience

By Guido Luzi

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More about us