The Multiple Aperture SAR Interferometry (MAI) Technique for the Detection of Large Ground Displacement Dynamics: An Overview

Abstract

This work presents an overview of the multiple aperture synthetic aperture radar interferometric (MAI) technique, which is primarily used to measure the along-track components of the Earth’s surface deformation, by investigating its capabilities and potential applications. Such a method is widely used to monitor the time evolution of ground surface changes in areas with large deformations (e.g., due to glaciers movements or seismic episodes), permitting one to discriminate the three-dimensional (up–down, east–west, north–south) components of the Earth’s surface displacements. The MAI technique relies on the spectral diversity (SD) method, which consists of splitting the azimuth (range) Synthetic Aperture RADAR (SAR) signal spectrum into separate sub-bands to get an estimate of the surface displacement along the azimuth (sensor line-of-sight (LOS)) direction. Moreover, the SD techniques are also used to correct the atmospheric phase screen (APS) artefacts (e.g., the ionospheric and water vapor phase distortion effects) that corrupt surface displacement time-series obtained by currently available multi-temporal InSAR (MT-InSAR) tools. More recently, the SD methods have also been exploited for the fine co-registration of SAR data acquired with the Terrain Observation with Progressive Scans (TOPS) mode. This work is primarily devoted to illustrating the underlying rationale and effectiveness of the MAI and SD techniques as well as their applications. In addition, we present an innovative method to combine complementary information of the ground deformation collected from multi-orbit/multi-track satellite observations. In particular, the presented technique complements the recently developed Minimum Acceleration combination (MinA) method with MAI-driven azimuthal ground deformation measurements to obtain the time-series of the 3-D components of the deformation in areas affected by large deformation episodes. Experimental results encompass several case studies. The validity and relevance of the presented approaches are clearly demonstrated in the context of geospatial analyses.

1. Introduction

The exploitation of Earth observation (EO) methods, based on the use of instruments operating in the microwave region of the electromagnetic spectrum, represents a common practice nowadays in the scientific community [1,2,3,4,5,6]. In this context, the development of new remote sensing techniques and the consolidation of well-established ones is fostered by increasing amounts of EO data collected by several sensors (on-board space or aerial vectors) that have emerged in recent years. A significant role is performed by the technologies based on their use of synthetic aperture radar (SAR) data that allow long-lasting and extensive monitoring of the Earth’s surface. SAR sensors are active instruments [7] that have been extensively employed for studying and monitoring the Earth’s surface [8,9,10], as well as the surfaces of other celestial bodies such as the solar system planets [11,12,13]. SAR instruments exploit the interaction of packets of electromagnetic waves emitted by RADAR sensors that are backscattered from objects within an observed area. The signals used by radar sensors usually lie in the microwave band (i.e., in the range of wavelengths between 1 cm and 1 m). This means that radar instruments can be used efficiently under any meteorological condition (i.e., in the presence of a significant cloud cover) with a full day-and-night operational capability. In this framework, several techniques have been developed to process and analyze SAR data. Among these, in recent decades, SAR interferometry techniques (InSAR) [14,15,16] have gained a role of particular relevance in the study of the Earth’s surface displacement phenomena. Ground deformations can be induced by numerous natural phenomena such as broad fractures due to big earthquakes [17,18,19,20,21,22], slow-moving landslides [23,24,25], eruptions of active volcanoes [26,27,28,29,30] and glacier movements [31,32]. Human beings also severely contribute to changes in the environment. In this framework, InSAR can monitor changes in highly urbanized zones due (for instance) to ground-water extraction as well as the deterioration of buildings and private and public facilities [33,34,35,36,37,38,39]. InSAR can monitor these modifications over time with high accuracy, allowing the retrieval of useful information for world governments to mitigate the risks for people and the environment. However, one of the main limitations of the InSAR methodology is that only sensor-to-target line-of-sight (LOS) projections of ground deformation can be detected and analyzed. Nevertheless, the combination of complementary information obtained by multi-satellite/multi-track InSAR observations can allow the discrimination of overall ground deformations in 3-D profiles (up–down, east–west, north–south). However, almost all modern spaceborne satellites fly along near-polar orbits (i.e., parallel to the north–south direction). This condition poses a problem in discriminating the north–south components of ground deformations measured from LOS-projected deformation signals, which have a low sensitivity to north–south displacements. Such a problem has been overcome in recent years through the application of pixel-offset (PO) tracking [40] and multiple aperture SAR interferometry (MAI) [41] approaches, which both permit computing the dynamics of large along-track ground surface changes. PO techniques rely on the measurement of local sub-pixel offsets between a couple of amplitude SAR images. As a consequence, the accuracy of estimating the azimuth/range offsets with amplitude matching can only reach a fraction (i.e., 1/30th) of the pixel spacing, which means the smaller the pixel size, the higher the accuracy. In a practical sense, considering the first-generation SAR satellites (e.g., the Advanced Synthetic Aperture Radar (ASAR), mounted onboard the Envisat platform /ASAR), PO can reach accuracy values that are at best greater than 10 cm. PO methods initially were mostly applied to studying single deformation events [32,40,42,43,44,45,46], but recently some efforts have also been made to extend these techniques to investigating SAR Doppler anomalies in focused SAR data [47] and for the generation of 3-D ground displacement time-series [48]. The latter is obtained by combining the information derived from both sequences of multi-temporal DInSAR interferograms and amplitude maps (e.g., [48]). A substantial improvement in measuring along-track displacement has been made with the MAI technique [41], which exploits spectral diversity (SD) processing [49], consisting of splitting the azimuth (range) SAR signal spectrum into separate sub-bands, to get an estimate of the surface displacement along the azimuth (sensor line-of-sight) direction. In particular, the larger the spectral separation between sub-bands, the higher the accuracy of the attainable measurements. Splitting band approaches were initially exploited by the oceanography community to infer direct measurements of ocean surface velocity, analyses and interpreting the Doppler shift between two RADAR SAR images acquired over the ocean [50,51]. The accuracy of the PO and MAI methods has been assessed and cross-examined in several studies [41,52,53]. The accuracy of MAI measurements has had a range between 3 and 7 cm, depending on the SAR sensor used, processing parameters and coherence of the MAI. Considering the parameters of an L-band system (e.g., the Phased Array type L-band Synthetic Aperture Radar (PALSAR 1-2), mounted onboard the Advanced Land Observing Satellite (ALOS) 1-2 platforms ALOS-PALSAR), the MAI technique can reach an accuracy of about 3–4 cm. MAI/PO techniques can complement results obtained using advanced multi-temporal InSAR approaches for the generation of ground deformation time-series [54,55,56,57,58,59,60,61,62,63,64]. A few solutions for acquired 3-D deformation time-series by combining multi-satellite/multi-observation (ascending/descending passes) data tracks have been proposed in the literature [65,66,67,68,69,70,71,72,73,74,75]. Among these, one method called Minimum Acceleration (MinA) combination [75,76] allows the generation of 2-D (3-D) displacement time-series of the Earth’s surface. This is achieved by combining LOS-projected deformation time-series obtained from different SAR platforms by imposing that the reconstructed 2-D (3-D) components of the ground deformation have a minimum acceleration. Noteworthy, the SD approaches are also currently used extensively to improve the performance of co-registration of the SAR images acquired with the novel Terrain Observation with Progressive Scans (TOPS) mode [77,78,79,80,81], which is the acquisition mode of the interferometric wide (IW) Sentinel-1 data [82]. Moreover, SD methods are used to estimate and correct atmospheric phase artefacts [83,84,85,86,87] that corrupt InSAR-driven ground deformation measurements.

In this work, we first shortly summarize the fundamentals of InSAR technology. Next, we introduce the principles of the SD and MAI techniques, discussing their applications in real scenarios and addressing their theoretical and empirical accuracies. A significant concern is mainly on the estimation of the along-track (north–south) components of the ground deformations. Furthermore, we propose a new method for combining multi-track (multi-satellite) along-track ground displacement time-series, as obtained through advanced, multi-temporal MAI approaches (e.g., [88]), and recovering the combined time-series of the ground displacement along the north–south direction. To this end, we have adequately adapted the technique initially proposed in [75] to overcome the challenge of combing along-track MAI-driven measurements. Henceforth, by applying both the proposed method and original MinA technique, we have retrieved the north–south and the two-dimensional (up–down, east–west) ground displacement time-series of the investigated areas, thus recovering the (whole) 3-D dynamics of the observed deformation phenomena. Experimental results are presented by processing sequences of SAR data collected by the Envisat/ASAR, COSMO-SkyMed and Copernicus Sentinel-1 radar sensors relevant to the case studies of Afar Rift depression, Ethiopia and the 2019 Ridgecrest earthquakes, California. The results of some investigations proposed by other authors in recent years are also referenced throughout this manuscript.

The paper is organized as follows. Section 2 illustrates the InSAR rationale, emphasizing key strengths and potential limitations. Section 3 provides an overview of the SD and multi-aperture InSAR techniques for generating ground deformation maps related to single episodes as well as for the production of time-series of deformation by multi-temporal SAR acquisitions. In Section 4, we present a method for combining MAI-driven deformation measurements recovered from multiple data tracks to get enhanced N-S ground deformation time-series. Discussion and conclusions are finally addressed in Section 5.

2. Fundamentals of InSAR Technology and Applications

SAR Interferometry (InSAR) technology [14,15,16] plays a significant role in the context of microwave remote sensing. The InSAR methodology is applied to retrieve information on the state of the Earth’s surface, its topography and its modifications. It relies on the extraction of the phase difference between two complex-valued SAR images related to the same imaged scene, usually referred to as master and slave, respectively. The multiplication between the master image and the complex conjugate of the slave image leads to the generation of a complex image (i.e., the complex interferogram) whose phase value is the phase difference between the two interfering SAR images. Such a complex interferogram takes into account the amplitude and phase returns of the observed target, and also contains a significant phase contribution due to the optical path from the target to the sensor at the collected times. Such a geometric phase contribution catches useful information on the geometry of the observed scene, such as the topography of the imaged scene or modifications (i.e., the ground displacement) that have occurred along the sensor-to-target line-of-sight illumination. Depending on the relative positions of the sensing transmitting/receiving antennas, different InSAR investigations configurations arise (e.g., [7,89,90]):

(1)

Across-track interferometry. In this case, the radar carrier has one or two sensors mounted on-board that are spaced along the across-track direction. Depending on the number of antennas on-board (one or two), the inferring SAR data pair can be acquired from different positions at different times (repeat pass interferometry), or from different locations at the same time (single-pass interferometry) (see Figure 1A,B). Across-track InSAR configuration is mainly used to map the Earth’s crust and detect and monitor its changes [21,91,92,93].

(2)

Along-track interferometry (ATI). In this case, two sensing antennas that are spaced in the along-track direction are mounted on-board the same carrier. Accordingly, two interfering SAR images can be acquired at different times from the same position (see Figure 1A,C). Along-track interferometry (ATI) configuration is usually used to estimate the radial motion of a surface point. In such a context, the ATI was firstly used for mapping the tidal ocean surface currents [94,95,96,97]. ATI is also principally applied for traffic monitoring with spaceborne SAR data [98,99,100].

We review the basic principles of across track interferometry in the following subsections. Interested readers can find a more comprehensive overview of the InSAR techniques and their applications in the literature (e.g., [101,102]).

InSAR for Topography Estimation

The SAR interferometry technique exploits the difference of viewing angles in two (or more) acquisitions (e.g., across-track configuration) to estimate the topography of an observed scene (e.g., [9]). Considering repeat-pass across-track configuration (see Figure 1B), two interfering SAR images can be expressed in complex notation as:

$$I_1\left(x,r\right)={\gamma }_1\left(x,r\right)e^{j\frac{4\pi }{\lambda }R}$$

(1a)

$$I_2\left(x,r\right)={\gamma }_2\left(x,r\right)e^{j\left(\frac{4\pi }{\lambda }R+{\delta }_R\right)}$$

(1b)

where ${\gamma }_{1,2}\left(x,r\right)$ represents the complex-valued reflectivity function related to the target of azimuth and range coordinates $\left(x,r\right)$ and $j=\sqrt{-1}$ is the imaginary unit. Note that R and R + ${\delta }_R$ are the range distances between the sensor and the target in the two images, with ${\delta }_R$ being the additive range distance term due to slightly different viewing geometries. The phase difference of the two SAR images is obtained by extracting the phase of the complex interferogram [9], which is given by the complex conjugate product operation between the two SAR images, as follows:

$$I\left(x,y\right)=I_1\left(x,y\right)·I_2^{*}\left(x,y\right)=\left|I_1\left(x,y\right)·I_2^{*}\left(x,y\right)\right|exp\left[j\frac{4\pi }{\lambda }{\delta }_R\right]$$

(2)

By taking into account the acquisition geometry shown in Figure 2 and applying the law of cosines, after trivial mathematical manipulations [103], it is straightforward to demonstrate that:

$$\Delta \varphi =\frac{4\pi }{\lambda }{\delta }_R\cong \frac{4\pi }{\lambda }bsin\left(\vartheta -\alpha \right)$$

(3)

where $\Delta \varphi$ is the extracted interferometric phase. Note that $\lambda$ is the operational wavelength, b is the interferometric baseline (i.e., the distance between orbital positions from which the two SAR images are taken), $\vartheta$ is the sensor side-looking angle and $\alpha$ is the angle the baseline makes with respect to a horizontal reference plane. Equation (3) contains useful information on the target height $h$. Indeed, by expanding Equation (3) around the angular position $\vartheta ={\vartheta }_0$, representing the illumination angle for flat terrain (i.e., h=0), we have:

$$\Delta \varphi =\frac{4\pi }{\lambda }·b\left[sin\left({\vartheta }_0-\alpha \right)+\cos \left({\vartheta }_0-\alpha \right)\left(\vartheta -{\vartheta }_0\right)\right]\cong \Delta {\varphi }_{flat}-\frac{4\pi }{\lambda }·\frac{b_{\perp }}{Rsin\left(\vartheta \right)}h=\Delta {\varphi }_{flat}-\Delta {\varphi }_{topo}$$

(4)

where $\Delta {\varphi }_{flat}$ is the so-called flat-Earth phase contribution, and $b_{\perp }$ is the component of the baseline that is perpendicular to the radar-to-sensor line-of-sight direction. The first term on the right-hand side of Equation (4), representing the flat-Earth phase term, can be derived (e.g., [93]) by computing the local phase frequency in the range direction as:

$$\Delta {\varphi }_{flat}=-\frac{4\pi }{\lambda }·\frac{b_{\perp }}{Rtan\left({\vartheta }_0-\mathsf{\Omega }\right)}\rho$$

(5)

where $\mathsf{\Omega }$ is the local slope of the terrain and $\rho =r$ is the pixel range coordinate of the imaged target. Conversely, the second phase term is the topographic phase signature $\Delta {\varphi }_{topo}$. More general relations arise if the curvature of the Earth is also taken into account (e.g., [7]). If the range-dependent phase term is compensated for in Equation (4), the topography of the imaged scene h is finally computed as [9]:

$$h=-\frac{\lambda Rsin\theta }{b_{\perp }}·\frac{\Delta {\phi }_{topo}}{4\pi }$$

(6)

From Equation (6), the height $h$ can be estimated with an accuracy that is baseline-dependent; this means that the larger is the perpendicular baseline, the more accurate the estimation of topography will be [7]. However, very large baseline values are responsible for severe decorrelation noise artefacts [104]. They corrupt the measured interferometric phase $\Delta \varphi$ and, subsequently, lead to incorrect height profile estimations. It is also worth highlighting that the interferometric phase is restricted to the [-π,π[ interval. The phase unwrapping (PhU) operation is a crucial step in any InSAR processing tool, and several approaches have been developed [105,106,107,108]. PhU involves the searching of (unknown) 2π-integer multiples that must be added to wrapped phases to compute unwrapped (full) phases.

As an evolution of the InSAR approach, a processing methodology known to as the differential synthetic aperture interferometric (DInSAR) technique [9,15,89] has been developed. It represents a common practice nowadays in the remote sensing scientific community for the detection and monitoring of the Earth’s surface displacement phenomena. The key factor of the DInSAR technique with respect to other conventional approaches (e.g., GPS and levelling measurement campaigns) is that it allows the continuous monitoring of displacement phenomena with a dense grid of measurement points. To convey the Earth’s surface displacements, the DInSAR technique relies on the generation of a so-called differential SAR interferogram [89]. Considering the acquisition geometry depicted in Figure 3, it is assumed that two complex-valued SAR images of the imaged scene were acquired at different times (see Figure 1B) and in different orbital positions. In a case where the ground surface is displaced by the 3-D vector d between the two flight passages of the sensor along the scene, the computed interferometric phase difference between the two SAR images—namely, $\Delta \varphi$—is made up of two main phase components:

$$\Delta \varphi =\Delta {\varphi }^{\prime }+\Delta {\varphi }_{defo}=\Delta {\varphi }^{\prime }+\frac{4\pi }{\lambda }{d}_{LOS}$$

(7)

where $\Delta {\varphi }^{\prime }$ is the phase contribution in the absence of deformation, as expressed in Equation (4), and $\Delta {\varphi }_{defo}$ is the additional phase contribution related to the ground displacement. Therefore, the phase signal associated with the deformation can be recovered by synthetically reconstructing the phase component contribution $\Delta {\varphi }^{\prime }$ and subtracting it, modulo $2\pi$, to the measured phase difference. In particular, $d_{LOS}$ is the projection of the displacement vector d along the sensor-to-target LOS direction. The synthetic phase term $\Delta {\varphi }^{\prime }$ is simulated from the knowledge of external information on the acquisition geometry such as a digital elevation model (DEM) of the observed area, the orbit state-vectors and the operational parameters of the radar instrument [7,9]. In a more general case, the differential interferometric phase is expressed as:

$$\Delta \varphi =\Delta {\varphi }_{defo}+\Delta {\varphi }_{topo}+\Delta {\varphi }_{orb}+\Delta {\varphi }_{atmo}+\Delta {\varphi }_{noise}$$

(8)

where $\Delta {\varphi }_{topo}=\frac{4\pi }{\lambda }·\frac{b_{\perp }}{rsin\theta }\Delta z$ represents the residual topography phase component induced by the DEM errors $\Delta z$; $\Delta {\varphi }_{orb}$ represents the residual phase term related to inaccurate orbital parameter information, which is used to estimate the aforementioned synthetic phase term; $\Delta {\varphi }_{atmo}$ denotes the phase components relative to the propagation variation of the RADAR signal caused by the absorption effects of Earth’s atmosphere; and, finally, $\Delta {\varphi }_{noise}$ is representative of additive noise contributions (i.e., spatial and temporal decorrelation, incorrect focusing of SAR raw data, etc.), that corrupt the interferometric phase [104,109]. In particular, the spatial decorrelation is less pronounced for short perpendicular baseline configurations.

From Equation (7), the accuracy of the displacement measurements depends on the operational wavelength $\lambda$. In this case, the maximum detectable unambiguous displacement that corresponds to a full phase cycle of $2\pi$ is equal to $\frac{\lambda }{2}$, denoting measurement accuracy of a fraction of the employed wavelength that is dependent on the amount of noise present in the measured phase.

Interested readers can find an overview of DInSAR and its main applications in [93,110].

3. Spectral Diversity and Multiple Aperture Interferometry

In this section we present spectral diversity approaches that can be exploited by the multiple aperture interferometry methods and used for the computation of large deformation phenomena of the solid Earth, as well as for correction of the co-registration errors of TOPS mode Sentinel-1 SAR data.

3.1. Spectral Diversity

Spectral diversity (SD) techniques [111,112] represent a class of methods based on the calculations of the spectral properties of pairs of complex-valued SAR images. The basic method consists of calculating the spectral separation of sub-bands of SAR images obtained by cutting two different frequency slices along the azimuth (Doppler band) of range directions, or by splitting the azimuth (range) spectrum of a full-band SAR image using a proper band-pass filter. Let us describe the SD process by initially referring to the azimuth case, where we assume the availability of two full-band Single Look Complex (SLC) images, referred to as the master and slave images $I_1\left(x,r\right)$ and $I_2\left(x,r\right)$, respectively. If the two SLC images are focused (squinted geometry) with the same Doppler centroid (DC) frequency—namely, $f_{DC}$—the two SAR images can be expressed after the co-registration step [85] as:

$$I_M\left(x,r\right)={\gamma }_M\left(x,r\right)·exp\left[j\frac{4\pi }{c}{f}_c·r\right]·exp\left[j2\pi \frac{f_{DC}}{v}·x\right]$$

(9a)

$$I_S\left(x,r\right)={\gamma }_S\left(x,r\right)·exp\left[j\frac{4\pi }{c}{f}_c·\left(r+\delta r\right)\right]·exp\left[j2\pi \frac{f_{DC}}{v}·\left(x+\delta x\right)\right]$$

(9b)

where $f_c$ is the radar carrier frequency, $v$ is the platform velocity and $\delta x$ represents the azimuthal misregistration of the slave image with respect to the reference master image. Accordingly, an additional interferometric phase term is present. Indeed,

$$\Delta \varphi =\frac{4\pi }{\lambda }{d}_{LOS}+2\pi \frac{f_{DC}}{v}\delta x$$

(10)

Consequently, if two azimuth sub-bands are centered on low and high DC frequencies ($f_{DC}^{-}$ and $f_{DC}^{+}$, respectively (see Figure 4), four SAR images can be formed:

$$I_M^{-}\left(x,r\right)={\gamma }_M\left(x,r\right)·exp\left[j\frac{4\pi }{c}{f}_c·r\right]·exp\left[j2\pi \frac{f_{DC}^{-}}{v}·x\right]$$

(11a)

$$I_S^{-}\left(x,r\right)={\gamma }_S\left(x,r\right)·exp\left[j\frac{4\pi }{c}{f}_c·\left(r+\delta r\right)\right]·exp\left[j2\pi \frac{f_{DC}^{-}}{v}·\left(x+\delta x\right)\right]$$

(11b)

$$I_M^{+}\left(x,r\right)={\gamma }_M\left(x,r\right)·exp\left[j\frac{4\pi }{c}{f}_c·r\right]·exp\left[j2\pi \frac{f_{DC}^{+}}{v}·x\right]$$

(11c)

$$I_S^{+}\left(x,r\right)={\gamma }_S\left(x,r\right)·exp\left[j\frac{4\pi }{c}{f}_c·\left(r+\delta r\right)\right]·exp\left[j2\pi \frac{f_{DC}^{+}}{v}·\left(x+\delta x\right)\right]$$

(11d)

This means that, for each image of the SAR data pair, two images are determined that are related to the low and high components of SAR data pair spectra. Subsequently, from these two data pairs, two interferograms are generated. After the split-spectrum operation, the interferometric phase terms related to the lower and higher azimuthal sub-bands are expressed as follows:

$$\Delta {\varphi }^{-}=\frac{4\pi }{\lambda }{d}_{LOS}+2\pi \frac{f_{DC}^{-}}{v}\delta x$$

(12a)

$$\Delta {\varphi }^{+}=\frac{4\pi }{\lambda }{d}_{LOS}+2\pi \frac{f_{DC}^{+}}{v}\delta x$$

(12b)

Therefore, the wrapped phase difference between the lower and higher spectral bands is

$$\Delta \tilde{\varphi }=Wr\left(\Delta {\varphi }^{-}-\Delta {\varphi }^{+}\right)=2\pi \frac{f_{DC}^{-}-f_{DC}^{+}}{v}\delta x$$

(13)

where $Wr(·)$ is the wrapping operator that wraps out the phase in the $\left[-\pi ,\pi \right[$ interval. As a consequence, the azimuth misalignment $\delta x$ can be computed from Equation (13). Assuming that such a misalignment is small enough to avoid the measured phase becoming ambiguous, the term $\delta x$ can easily be derived as

$$\delta x=\frac{v\Delta \tilde{\varphi }}{2\pi \left(f_{DC}^{-}-f_{DC}^{+}\right)}$$

(14)

From Equation (14), it is evident that the larger the spectral separation between the central frequencies of the two sub-bands, the more accurate the estimation of $\delta x$ is. We would like to stress that in non-stationary scenarios, the imaged target at azimuth location x is subject to displacement, —namely, $\delta x_{disp}$—and, in this case, the measured $\delta x$ term contains both a contribution due to the azimuthal misregistration $\delta x_{misr}$ and the azimuthal $\delta x_{disp}$ displacement. Of course, misalignments of SAR data pairs can be thoroughly controlled in the imaging process to mitigate phase decorrelation effects. Perfectly aligned SAR data are then used to infer geophysical displacement measurements, as they have been via the multiple aperture SAR interferometry (MAI) technique in cases of large ruptures in the ground [41,113].

We would like to remark that SD and splitting-band methods can also be applied along the range direction by extracting the two sub-bands in the range. In particular, we begin with two range sub-bands centered on the two frequencies $f_c^{-}$ and $f_c^{+}$ (see Figure 5), which can be expressed, similarly to Equations (11a-d), as

$$I_M^{-}\left(x,r\right)={\gamma }_M\left(x,r\right)·exp\left[j\frac{4\pi }{c}{f}_c^{-}·r\right]·exp\left[j2\pi \frac{f_{DC}}{v}·x\right]$$

(15a)

$$I_M^{-}\left(x,r\right)={\gamma }_M\left(x,r\right)·exp\left[j\frac{4\pi }{c}{f}_c^{-}·r\right]·exp\left[j2\pi \frac{f_{DC}}{v}·x\right]$$

(15b)

$$I_M^{+}\left(x,r\right)={\gamma }_M\left(x,r\right)·exp\left[j\frac{4\pi }{c}{f}_c^{+}·r\right]·exp\left[j2\pi \frac{f_{DC}}{v}·x\right]$$

(15c)

$$I_S^{+}\left(x,r\right)={\gamma }_S\left(x,r\right)·exp\left[j\frac{4\pi }{c}{f}_c^{+}·\left(r+\delta r\right)\right]·exp\left[j2\pi \frac{f_{DC}}{v}·\left(x+\delta x\right)\right]$$

(15d)

In this case,

$$\Delta {\varphi }^{-}=\frac{4\pi }{c}{f}_c^{-}·d_{LOS}+2\pi \frac{f_{DC}}{v}\delta x+\Delta {\varphi }_n^{-}$$

(16a)

$$\Delta {\varphi }^{+}=\frac{4\pi }{c}{f}_c^{+}·d_{LOS}+2\pi \frac{f_{DC}}{v}\delta x+\Delta {\varphi }_n^{+}$$

(16b)

where $\Delta {\varphi }_n^{-}$ and $\Delta {\varphi }_n^{+}$ are the negative and positive sub-band phase noise contributions, respectively. Accordingly, the wrapped phase difference of the higher- and lower-band phase terms gets an alternative estimate of the LOS displacement as

$$d_{LOS}=\frac{c\Delta \tilde{\varphi }}{4\pi \left(f_c^{-}-f_c^{+}\right)}$$

(17)

where, similarly to Equation (13), $\Delta \tilde{\varphi }$ represents the wrapped phase difference $Wr\left(\Delta {\varphi }^{-}-\Delta {\varphi }^{+}\right)$ between the two phase looks (see Equation (16b)). Note that $d_{LOS}$ is assumed to be small enough to avoid phase measurements ambiguities. In this framework, a procedure that is known as spectral shift filtering [114] implements proper filters that are tuned on specifically chosen different central frequencies to compensate for the geometrical decorrelation effects that corrupt the interferometric phase. This strategy has been applied correctly to process coupled SAR images collected by the ERS and the Envisat systems, which were characterized by a spectral shift of 31 MHz [115,116].

As final remark, we would like to stress that only a fraction $\beta$ of the total bandwidth of the two images is used after the application of the split-spectrum operation. Accordingly, being the spatial resolution inversely proportional to the band of the signal [89], the spatial resolution of the two sub-band interferograms is reduced by the same fraction.

3.2. Multiple Aperture Interferometry Principles

As said in Section 2, InSAR techniques have successfully been used to study and analyze many different phenomena that characterize the Earth’s surface. However, these techniques are limited to measuring only the 1-D components of ground deformation along the RADAR LOS direction. Conversely, enabling 3-D displacement measurements is crucial for better analyses and to study the deformation phenomena that characterize areas of interest. In the last two decades, many studies have been carried out enabling the measurement of the 3-D (up–down (U-D), east–west (E-W) and north–south (N-S)) ground displacement [66,67,117]. Multi-platform, multi-track InSAR techniques have been developed to discriminate the U-D and E-W components of surface displacement, but estimation of N-S ground displacement data has been more challenging. In this respect, SAR amplitude tracking methods—namely, pixel offset (PO) [32,42,48,65,68,118] and multiple aperture interferometry (MAI) [41,49,52,53]—have been applied.

It has been demonstrated that the measurement accuracy of N-S displacement that is attainable with PO methods is low. In particular, it reaches only a fraction of pixel spacing (1/30th), as the azimuthal pixel spacing of first-generation and present-day SAR instruments is in the order of 3–5 cm [32,44,45,46,48]. For example, taking into account that the azimuth pixel spacing of the TerraSAR-X sensor is of about 3 m, the expected accuracy of PO measurements is about 10 cm. Of course, fine and ultra-fine SAR data with enhanced spatial resolutions are envisaged to apply the PO method. It is worth noting that a substantial improvement in the discrimination of the N-S ground displacement component from InSAR data is a key factor in performing extended geophysical investigations.

3.3. Multiple Aperture Interferometry for the Along-Track Measurement

Monitoring the Earth’s surface phenomena requires more comprehensive analyses, such as estimating the 3-D field of ground surface displacement that affects a specific area. Since SAR satellite platforms travel in almost near-polar orbits, the main problem in the discrimination of the north–south component of deformation (approximately along-track displacement) is that the projection of the ground displacement along the north–south direction is not accurate enough [41].

In recent years, a big improvement in the discrimination of along-track displacement has been through development of the MAI technique. This technique uses a split beam process called “spectral diversity” (see Section 3.1) to determine backward- and forward-looking interferograms [111] from a SAR data pair that is related to the same scene but were acquired at different times (see Figure 6). Considering Equation (13), the backward- and forward-looking interferograms can be expressed as follows:

$${\varphi }_{fw}=\frac{2\pi }{v}·f_{DC}^{+}x$$

(18a)

$${\varphi }_{bw}=\frac{2\pi }{v}·f_{DC}^{-}x$$

(18b)

Specifically,

$$f_{DC}^{+}=f_{DC,c}+n·\frac{\Delta f_D}{2}$$

(19a)

$$f_{DC}^{-}=f_{DC,c}-n·\frac{\Delta f_D}{2}$$

(19b)

where $f_{DC,c}=\frac{f_{DC,m}+f_{DC,s}}{2}$ and $\Delta f_D=\frac{2v}{l}$ are the average Doppler centroid (DC) frequencies of the two full-bandwidth SAR images, with $f_{DC,m}$ and $f_{DC,s}$ being the DCs of the master and slave SAR images, respectively. Note that $\Delta f_D$ is the effective Doppler bandwidth, $l$ indicates the azimuthal antenna length (see Figure 7A) and $n$ is the fraction of the azimuth bandwidth (e.g., $\Delta f_D$; see Figure 7A) used in the split-beam process, also indicated as the normalized squint angle of the sub-aperture process. Consequently, starting from the measured phase difference between the forward- and backward-looking interferometric phases, the MAI phase is determined as follows:

$${\varphi }_{MAI}=Wr({\varphi }_{fw}-{\varphi }_{bw})=\frac{2\pi }{v}n·\Delta f_D·\Delta x=\frac{2\pi }{v}n\left(\frac{2v}{l}\right)\Delta x=\frac{4\pi }{l}n\Delta x+\eta$$

(20)

In Equation (20), it is clear that the MAI phase is proportional to the along-track deformation $\Delta x$, where, for the sake of simplicity, we have imposed $f_{DC,c}=0$, and $\eta$ is the interferometric noise. Naturally, Equation (20) represents the simplest case. In particular, during the MAI interferogram determination, the SD process generates the forward- and backward-looking interferometric phases, respectively, splitting the full azimuthal spectral bandwidth of the master and slave images (see Figure 7B).

If the master and slave acquisition orbits are not perfectly parallel when considering the backward- and forward-looking imaging geometries (see Figure 6), different flat-Earth and topographic phase contributions arise in the forward- and backward-looking interferograms. Considering Equations (4) and (5), Jung et al. [113] addressed this problem by considering the expression of the additional spurious phase terms related to the flat-Earth and the topographic contributions, which are related to the baseline difference between the forward- and backward-looking imaging geometries, namely, $\Delta B_{\perp }$:

$${\varphi }_{flat}=\frac{4\pi \Delta B_{\perp }}{\lambda R\tan \left(\theta \right)}\rho$$

(21a)

$${\varphi }_{topo}=\frac{4\pi \Delta B_{\perp }}{\lambda R\sin \left(\theta \right)}h$$

(21b)

where $\rho$ and $h$ are the slant range and topographic height, respectively. The atmospheric phase screen (APS) may be considered negligible in the final MAI phase, as APS is influential in the same way both the backward- and forward-looking interferograms are.

In this context, Jung et al. [113] presented a MAI improvement process on which they estimated and corrected for the mentioned flat-Earth and topographic spurious phase terms. Afterwards, we will present two case studies related to the application of the MAI technique to the investigation of a single deformation episode.

The first case study concerns the Afar depression. Afar is a region situated in Ethiopia, Africa, that was affected in 2005 by a big earthquake associated with the Ethiopian Great Rift fault mechanism [119,120].

Specifically, from September to October 2005, a seismic sequence consisting of 131 events affected the Afar region as a result of the Dabbahu volcano eruption (see Figure 8A). Since the initial main events of the drifting episode, a further sequence of 13 discrete dyke events was detected from 2005 to 2010 along the entire Dabbahu rift segment [119].

We performed an experiment by considering one pair of Envisat/ASAR data that captured the effects of the Dabbahu rifting episode. SAR data were acquired along ascending orbits (Track 200) on 19 December 2005 (master) and 25 August 2008 (slave), as shown in Figure 9A,B. Starting from this pair of SAR data, a MAI interferogram was generated (see Figure 10A). The forward- and backward-looking interferograms of the Dabbahu region area are shown in Figure 10B,C. In particular, we have highlighted an area in the MAI interferogram with a red circle that showed sensitive ground displacement and was in correspondence with a trace of the 2005 activated rift segment.

The second case study was performed by processing a couple of SAR data acquired over the north-east Ridgecrest town area, situated in California, USA, (see Figure 8B), which was struck by a big (Mw) 7.1 earthquake on 6 July 2019, the strongest in the region in at least two decades. The phenomena included other previous main shock events, and in the days and weeks that followed, thousands of aftershocks rumbled beneath Southern California [121].

In particular, starting from a single pair of COSMO-SkyMed acquisitions—the master collected on 4 July 2019, and the slave on 20 July 2019 (see Figure 11A,B)—a MAI interferogram (Figure 12A) was generated from a pair of forward- and backward-looking interferograms (Figure 12B,C). Specifically, in the interferogram depicted in false colors, we can clearly see the big fault related to the main seismic event circled in red (despite the presence of the interferometric noise) at the center.

Considering the parameters of the X-band COSMO-SkyMed RADAR instrument of the Italian Space Agency (ASI) and the Envisat/ASAR of European Space Agency (ESA)—with azimuthal antenna lengths of 5.7 and 11.1 m, respectively—we evaluated the along-track measurement enabled by the MAI technique. In particular, taking into account Equation (20) and neglecting the noise term (for the sake of simplicity), we observe that an ambiguous maximum along-track deformation $\Delta x$ of about $5.7$ and 11.1 m for every in $2\pi$ phase cycle can be measured. This result implicitly confirms the useful utility of the MAI technique for estimating and analyzing large displacements of the Earth’s surface. A sounder analysis of MAI measurement performance is presented in Section 3.4. The processing scheme of MAI operations is shown in Figure 13.

3.4. Multiple Aperture Interferometry Accuracy and Noise propagation

Based on Equation (20), the accuracy of the MAI technique can be theoretically estimated by the following Equation (22) that relates the standard deviation of the MAI phase measurement ${\sigma }_{{\phi }_{MAI}}$, with the standard deviation of the along-track deformation measurements ${\sigma }_x$ [113]:

$${\sigma }_x=\frac{l}{4\pi ·n}{\sigma }_{{\phi }_{MAI}}$$

(22)

From Equation (22), it is clear that the MAI measurement accuracy strictly depends on the standard deviation of the interferometric phase noise [112], which is given by

$${\sigma }_{{\phi }_{MAI}}=\sqrt{{\sigma }_{fw}^2+{\sigma }_{bw}^2-2{\sigma }_{fw,bw}}$$

(23)

where ${\sigma }_{fw}^2,{\sigma }_{bw}^2,{\sigma }_{fw,bw}$ are the variances of forward- and backward-looking interferograms and their co-variance, respectively. Thus, based on the assumption that the probability density function (pdf) of the single/multi-look interferometric phases can be expressed as a function of the correlation coefficient $\gamma$ [122], and assuming the forward- and backward-looking interferograms are two statistically independent random variables (e.g., ${\sigma }_{\phi ,fb}$ is zero) in the limit of Cramér–Rao [112] (which is valid for high values of $\gamma$), ${\sigma }_{\phi ,f}$ and ${\sigma }_{\phi ,b}$ can be inferred as follows:

$${\sigma }_{fw}={\sigma }_{bw}\approx \frac{1}{\sqrt{2NL}}·\frac{\sqrt{1-{\gamma }^2}}{\gamma }$$

(24)

where NL is the effective number of looks. Finally, using Equations (23) and (24), the standard deviation of the interferometric MAI phase is reduced to

$${\sigma }_{{\phi }_{MAI}}\approx \frac{1}{\sqrt{NL}}·\frac{\sqrt{1-{\gamma }^2}}{\gamma }$$

(25)

where we assume that the forward- and backward-looking interferograms are uncorrelated. From Equations (22) and (25), it is obvious that improving the coherence of the MAI interferogram and having a larger spectral separation between the sub-band both lead to better accuracy in the estimation of along-track deformation measurements. For the Cramér–Rao boundary, it was demonstrated in [77,112] that the optimal separation between the sub-bands is equal to two-thirds of the full bandwidth, and that:

$${\sigma }_{{\phi }_{MAI,CR}}\approx \sqrt{\frac{3}{NL}}·\frac{\sqrt{1-{\gamma }^2}}{\gamma }$$

(26)

By considering Equation (26), we can plot the standard deviation of the azimuthal displacement vs the spatial coherence of the MAI interferogram (see Figure 14A,B) for different values of the effective look numbers (NL), which comprise the operational parameters of the Envisat/ASAR, COSMO-SkyMed and Sentinel-1 platforms listed in

Table 1. Envisat/ASAR stripmap mode acquisition platform parameters.

Parameter	Value	Unit
Wavelength (Centre Frequency)	0,05624624	m
Band type	C
Pulse repetition frequency	1.652.4157	Hz
Polarisation options	Single VV, HH or Dual VV+HH, VV+VH, HH+HV
Azimuth antenna size	11.1	m
Incidence angle range	15° - 45°	deg.
Swath width	100	km
Resolution	6(az)x9(rg)	m
Azimuth pixel spacing	4.31	m
Range pixel spacing	7.8	m

,

Table 2. COSMO-SkyMed stripmap (ping-pong) mode acquisition platform parameters.

Parameter	Value	Unit
Wavelength (Centre Frequency)	0,031228381	m
Band type	X
Pulse repetition frequency	3.554.5024	Hz
Polarisation options	Single VV, HH or Dual VV+HH, VV+VH, HH+HV
Azimuth antenna size	5.6	m
Incidence angle range	25° - 50°	deg.
Swath width	30 km	km
Resolution	15(az)x15(rg)	m
Azimuth pixel spacing	2.53	m
Range pixel spacing	1,56	m

and

Table 3. Sentinel-1 operational parameters in interferometric wide (IW) acquisition mode.

Parameter	Value	Unit
Wavelength (Centre Frequency)	0.055465763	m
Band type	C
Pulse repetition frequency	486.4863	Hz
Polarisation options	Dual HH+VV, VV+HH or Single HH,VV
Azimuth antenna size	12.3	m
Incidence angle range	29.1° - 46.0°	deg.
Azimuth steering angle	±0.6°	deg.
Swath width	250	km
Number of sub-swaths	3
Resolution	5(rg)x20(az)	m
Azimuth pixel spacing	14.1	m
Range pixel spacing	2,3	m

, respectively. In particular, as outlined in Section 3.6, the Doppler bandwidth of the TOPS-mode SAR data is reduced for a single point target with respect to the conventional stripmap case by a factor of about three (see Equation (20)).

Interested readers can find a more comprehensive analysis by Jung et al. in [123].

3.5. Multiple Aperture Interferometry for the Generation of Along-Track Deformation Time-Series

SAR interferometric techniques (such as DInSAR and MAI) allow many geophysical events that characterize Earth’s surface to be studied in depth through the generation of single deformation maps. However, ground displacement phenomena are processes-characterized by gradual or sudden changes in the Earth’s surface elevation over time. In this context, to better characterize the displacement phenomena in an area of interest, an analysis of temporal modifications to ground displacement is necessary. Over the years, several methodologies have been developed that are useful for extending the use of the differential interferometric SAR techniques to time-monitor displacement phenomena [54,55,56,57,58,59,60,61,62,63]. These methods are based on the inversion of properly selected sequences of DInSAR interferograms, allowing the computation of LOS-projected deformation time-series. Two main categories of these techniques are present in literature. The class of the persistent scatterer (PS) methods [54,58] is based on the identification of coherent point-wise targets exhibiting high phase stability over a sequence of DInSAR interferograms that have been generated from a set of SAR data related to the same scene but acquired at different times. Conversely, the class of small baseline (SB) techniques [55,59,63,124] is devoted to the analysis of distributed targets on the ground by processing sequences of small baseline interferograms. More recently, a new advanced method for the characterization of DS targets, called SqueeSAR [61], has also been developed. SqueeSAR uses a statistical approach that exploits phase and amplitude information of a sequence’s SAR data to identify very highly coherent DS targets by the use of a maximum-likelihood optimization.

Very recently, temporal analysis methods have been used to extend the MAI technique and retrieve the temporal evolution of the along-track ground displacement. Specifically, a first attempt was made to estimate mean deformation along-track velocity by stacking a sequence M MAI interferograms $\varphi =\left[{\varphi }_{MAI}^0,{\varphi }_{MAI}^1,\dots ,{\varphi }_{MAI}^{M-1}\right]$ generated from a stack of N SAR images. This was done by writing out the following system of equations:

$$\left[\begin{array}{c}\Delta t_0\\ ⋮\\ \Delta t_{N-1}\end{array}\right]\left[v\right]=\left[\begin{array}{c}{\varphi }_{MAI}^0\\ ⋮\\ {\varphi }_{MAI}^{N-1}\end{array}\right]$$

(27)

where ${\varphi }_{MAI}^i$ and $\Delta t_i$ are the ith MAI interferogram of the sequence and its temporal baseline, respectively, and $v$ is the unknown deformation velocity to estimate. Finally, by resolving the system of Equation (27) in the least-squares (LS) sense, the deformation velocity is retrieved. A novel processing chain to stack a sequence of MAI interferograms was presented in Jo et al. [125] and used to measure the slow-moving azimuthal displacements of the Kilauea volcano on the Big Island, Hawaii. Taking into account that forward- and backward-looking interferograms are affected by noise, they suggested that—instead of estimating the MAI phase from the interferograms and applying the stacking procedure—the multi-temporal forward- and backward-looking residual interferograms could be individually stacked, in order to generate a MAI velocity map as:

$$v_{MAI}=\frac{l}{4\pi n}·\frac{\left[\left\{{{\displaystyle \sum }}_{i=1 }^N{\phi }_{fw,res}^i\right\}-\left\{{{\displaystyle \sum }}_{i=1 }^N{\phi }_{bw,res}^i\right\}\right]}{{{\displaystyle \sum }}_{i=1}^N{\Delta \mathrm{t}}_{\mathrm{i}}}$$

(28)

where ${\phi }_{fw.res}^i$ and ${\phi }_{bw.res}^i$ are the interferometric phases of ith residual forward- and backward-looking interferograms, respectively, which are generated by removing a full aperture multi-look differential interferogram from the forward- and backward-looking one; ${\Delta \mathrm{t}}_{\mathrm{i}}$ is the ith time duration; and N is the total number of forward- and backward-looking interferograms.

In particular, the theoretical variance of the estimated deformation velocity map measurements for the MAI stacking was evaluated to be:

$${\sigma }_{v,MAI\_stacking}=\frac{l}{4\pi n}·\frac{{\sigma }_{\varphi ,MAI}}{\sqrt{M}·\overline{\Delta t}}$$

(29)

where $\overline{\Delta t}$ is the average temporal baseline of the whole InSAR distribution. The three-dimensional (3-D) displacement velocity map of the Kilauea Volcano, obtained by combining multi-stacked DInSAR and MAI ground displacement rate maps, is shown in Figure 12 of [125].

Another important temporal study of the Earth’s surface displacement via the MAI technique was carried out by Gourmelen et al. [126]. In their work, the authors combined displacement velocity maps that were obtained by applying the DInSAR and MAI techniques to map the ice surface velocity of the Langjokull and Hofsjokull ice caps in Iceland in 1994. In particular, their approach improved the accuracy of the ice flow measurement by a factor of two for E-W, U-D measurements, and up to a factor of 10 for N-S measurements, as compared with a velocity solution based on InSAR and pixel offset (PO).

Another class of studies based on the use of MAI interferograms is the multi-temporal extension of the MAI approach, which is used to generate along-track ground deformation time-series [88].

In this instance, a network of M InSAR data pairs is identified from a set of N SAR images. If $t_M$ and $t_S$ represent the acquisition times of the master and slave images of a given interferogram, respectively, then the MAI phase is expressed via Equation (20) as:

$${\varphi }_{MAI}=\frac{4\pi }{l}n\left[\Delta x\left(t_S\right)-\Delta x\left(t_M\right)\right]$$

(30)

where $\Delta x\left(t_S\right)$ and $\Delta x\left(t_M\right)$ are the unknown along-track deformations at the master and slave image acquisition times, respectively. This leads to the solution of a system of linear equations that can be expressed, using matrix formalism, as:

$$\mathit{A}·\mathit{x}={\mathbf{\varphi }}_{\mathit{M}\mathit{A}\mathit{I}}$$

(31)

where $\mathit{x}$ represents the vector of the unknown time-series deformation, ${\mathbf{\varphi }}_{\mathit{M}\mathit{A}\mathit{I}}$ is the vector of the phase components relative to the MAI interferogram and $\mathit{A}$ is the incidence-like matrix of the InSAR network graph, whose jth row is defined as follows:

$$\{\begin{array}{c}\mathit{A}\left(j,I_{M,j}\right)=-\frac{4\pi }{l}n\\ \mathit{A}\left(j,I_{S,j}\right)= \frac{4\pi }{l}n\\ \begin{array}{cc}otherwise& 0\end{array}\end{array}$$

(32)

where in ${\mathit{I}}_{\mathit{M}}$ and ${\mathit{I}}_{\mathit{S}}$ are the M-length vectors of the indices of the master and the slave time acquisitions, respectively. The solution of the system of Equation (31) can be obtained in the LS sense as:

$$\tilde{\mathit{x}}={\mathit{A}}^{+}·{\mathbf{\varphi }}_{\mathit{M}\mathit{A}\mathit{I}}=\left({\left({\mathit{A}}^{\mathit{T}}·\mathit{A}\right)}^{-\mathbf{1}}·{\mathit{A}}^{\mathit{T}}\right)·{\mathbf{\varphi }}_{\mathit{M}\mathit{A}\mathit{I}}$$

(33)

where ${\mathit{A}}^{+}$ indicates the left pseudo-inverse of the matrix $\mathit{A}$. Note that $\mathit{A}$ is a full-rank matrix if all the SAR acquisitions form one single connected set of data. In a more general case where several different subsets of data are present, however, this can lead to a rank-deficient matrix $\mathit{A}$. In this situation, the same strategy adopted in the small baseline subset (SBAS) technique can be applied, which consists of a proper manipulation of the system of Equation (31) and the application of the singular value decomposition (SVD) method (see [57] for additional details). We note that, for typical values of the expected along-track displacement, the unwrapping operation is not required in our case. Furthermore, the phase residuals:

$$\mathit{e}=\mathit{A}·\tilde{\mathit{x}}-{\mathbf{\varphi }}_{\mathit{M}\mathit{A}\mathit{I}}$$

(34)

can be used to identify sets of well-processed SAR pixels. In particular, in this paper we propose exploiting the following temporal coherence factor to evaluate the quality performance of the inversion strategy:

$$\mathsf{{\rm Y}}=\frac{1}{N}\left|{\displaystyle \sum }_{i=0}^{N-1}exp\left[j{e}_i\right]\right|$$

(35)

where $e_i$ is the ith phase residual value, and $\mathsf{{\rm Y}}$ has the same mathematical expression of the temporal coherence factor initially proposed in [127]. However, in this case, the phase residuals are not a result of time-inconsistent phase unwrapping mistakes, but rather phase noise of the MAI interferograms.

The deformation time-series have been obtained by processing a sequence of 14 MAI interferograms generated from a sequence of seven Envisat/ASAR acquisitions of the Afar depression spanning the period between 9 December 2005 and 10 August 2009. The used data are listed in

Table 4. Envisat/ASAR acquisition data of the Afar depression zone. Date format is day/month/year.

Acquisition n.	Date
1	19/12/2005
2	27/02/2006
3	04/12/2006
4	10/09/2007
5	28/01/2008
6	25/08/2008
7	10/08/2009

, and the distribution of the 14 MAI interferograms in the time/perpendicular baseline plane (indicating a Delaunay graph network [128]), is depicted in Figure 15.

Figure 16 shows the mean deformation velocity of the Afar area (see Figure 16B) superimposed on a SAR amplitude image of the investigated area, and the relevant time-series of deformation (see Figure 16C–F). Only the deformation values of the coherent, well-processed SAR pixels are shown in the velocity map. They are identified by computing the map of the temporal coherence through Equation (35), as shown in Figure 16A.

The area in Figure 16B with a sensitive deformation (blue color) is located across the fault, and we highlighted four zones of interest that concerned the time-series analysis (points (B), (C), (D) and (E)). These major fault areas correspond to: (i) a major fault zone whose deformation time-series (Figure 16C) showed an almost constant along-track deformation rate annually, with a leap of about 50 cm in the deformation during the 2005–2006 time-period, when the main seismic event occurred; (ii) a zone of medium-magnitude deformation (Figure 16D–E), where we can perceive opposite behaviors, as if the deformations that characterize the AFAR depression should balance out somehow; and (iii) a zone of low-magnitude deformation (Figure 16F) where the deformation time-series remained constant.

A more quantitative estimate of the accuracy of the achieved along-track deformation time-series can be obtained by applying the basic principles of error noise propagation [129]. In particular, given a full covariance matrix for the vector of MAI phases ${\mathit{C}}_{MAI}={\sigma }_{MAI}^2{\mathit{I}}_{M\times M}$, the covariance matrix of the along-track deformation time-series is given by ${\mathit{C}}_X={\left({\mathit{A}}^T·{\mathit{C}}_{MAI}^{-1}·\mathit{A}\right)}^{-1}$, where we have assumed the simplified hypothesis that the MAI phases are independent and uncorrelated, and that the standard deviation of the MAI phase is given by Equation (22). A more extensive treatment is required in the case that matrix $\mathit{A}$ is rank deficient, and a given correlation is assumed between the MAI interferograms. This is a matter for further investigation, which would require a deep understanding of the statistical relationships among the interferograms, as well as a characterization of the studied PS and DS targets on the ground.

More recently, other advanced multi-temporal deformation investigations have been carried out using the MAI technique. Specifically, He et al. [88] performed a temporal analysis of big landslides that had occurred in the Fushun west open-pit mine (FWOPM), China. They mapped the two-dimensional deformation field time-series of landslide phenomena that occurred between 2007 and 2011 by coupling DInSAR-SBAS with MAI-SBAS. They validated the generated time-series (across-track and along-track directions) with GPS displacement measurements over Fushun mine. The time-series obtained by He et al. are presented in Figure 10 of their article [88]. In particular, the results in the two graphics demonstrate that the time-monitoring analysis of the DInSAR and MAI techniques returned information that was highly consistent with the landslide kinematic pattern measured by the GPS sensor in that area.

All of the analyses presented for the generation of along-track ground deformation time-series demonstrate that the key factor is the reduction of decorrelation noise errors. In Liu et al. [130], a very accurate investigation of the effects of these spurious phase components is presented that perturbs the final measurements. Specifically, Liu et al. classified these errors into three categories: (i) random error, (ii) systematic error and (iii) gross error. These categories provide a more comprehensive understanding of the errors in inverting interferometric DInSAR and MAI phases when computing time-series and 3-D maps of deformation of an observed scene.

3.6. Overview of the Enhanced Spectral Diversity for TOPS Mode SLC Image co-registration

An enhanced spectral diversity (ESD) approach has also been developed and extensively used for the fine co-registration of pairs of single-look-complex SAR images acquired through the TOPS mode (e.g., IW products). Indeed, the growing demands in having broader swath coverage led to the development of a new advanced imaging mode. In particular, the capability to steer the radar antenna beam along the azimuth direction led to the design of the Terrain Observation by Progressive Scans (TOPS) [131] mode, which is used to achieve wide-swath coverage. As with the conventional ScanSAR mode [132], the TOPS imaging mode achieves wide-swath coverage by switching the antenna beam among range directions from swath to swath (often referred to as sub-swaths) [132].

However, in the TOPS mode (in contrast with the ScanSAR), the antenna beam has an angular velocity $v_{rot}$ that rotates from backward to forward throughout the acquisition with respect to a virtual rotation center (see Figure 17). This allows all targets to be illuminated during the data acquisition duration (burst) within a large portion of the azimuth antenna pattern, thus mitigating the image quality degradation effects such as ambiguity to signal ratio or scalloping effects. However, these benefits are offset by a worsened azimuthal resolution of SLC-focused images [131], which are now reduced by a factor of $A=\frac{vg}{vs}$, where $vg$ and $vs$ are the ground and the satellite velocities, respectively.

Let us first introduce the problem by considering two SLC IW images, which are usually acquired over several bursts in three overlapped swaths [79]. The fine co-registration of the two images represents one of the main issues of the Sentinel-1 (S-1) SAR data processing: to achieve this task, it is widely recognized that the azimuth co-registration of the SAR images has to be as accurate as possible to achieve a maximum error not exceeding 0,001 azimuth samples [78,80,81,133,134]. Usually, a conventional geometric co-registration step is performed preceding every SLC burst image.

To introduce the problem, let us consider the two contiguous ith and (i+1)th bursts of the master and slave images—namely $M_i$ and $M_{i+1}$ and $S_i$ and $S_{i+1}$, respectively. In the TOPS case, the Doppler frequency varies linearly for every burst with azimuth time $t_a$—namely, $f_{dop}\left(t_a\right)=K_T·t_a$, where $K_T$ is the Doppler frequency rate. Some considerations of TOPS acquisition geometry have led to the conclusion that the Doppler frequency rate within each burst is as follows [81]:

$$K_T=\frac{K_a·K_{rot}}{K_a-K_{rot}}$$

(36)

where $K_a=-2v_{eff}^2/\lambda r$ is the Doppler rate of the target, and $K_{rot}=-2v_{eff}^2/\lambda r_{rot}=2v_s^2{k}_{\psi }/\lambda$. Finally, $k_{\psi }$ is the angular rotation rate. Note that $v_{eff}=v_g{v}_s$ is the effective velocity target, where $v_g$ and $v_s$ are the velocities of the ground and the satellite platform, respectively; $r$ and $r_{rot}$ are the sensor-to-target range distance and the distance between the target on the ground and the virtual rotation center of the acquisition, respectively; and $\lambda$ is the operational wavelength. The (i+1)th and the ith complex burst interferograms—namely, ${\mathsf{\Gamma }}_{i+1}$ and ${\mathsf{\Gamma }}_i$—can be obtained as follows (see Figure 18):

$${\mathsf{\Gamma }}_{i+1}=M_{i+1}·S_{i+1}^{*}$$

(37a)

$${\mathsf{\Gamma }}_i=M_i·S_i^{*}$$

(37b)

Accordingly, in the overlapped region between the two consecutive ith and (i+1)th bursts, a differential interferogram is formed as

$${\varphi }_{ovl}=arg\left\{{\mathsf{\Gamma }}_{i+1}·{\mathsf{\Gamma }}_i^{*}\right\}=Wr({\varphi }_{i+1}-{\varphi }_i)$$

(38)

where

$${\varphi }_{ovl}={\varphi }_{azDispl}+{\varphi }_{geomErr}+{\varphi }_{Misreg}$$

(39)

The first term on the right side of Equation (39), ${\varphi }_{azDispl}$, represents the contribution of the Earth’s surface ground displacement along the azimuth direction; the second term, ${\varphi }_{geomErr}$, is the residual geometrical phase due to problems of local co-registration during the geometrical coregistration phase; and the last term, ${\varphi }_{Misreg}$, accounts for the constant (time) azimuth misregistration $\Delta t$ after the co-registration processes. Given a stationary scene and assuming negligible residual geometrical phases, the differential phase over the overlapped region (see Figure 18) is

$${\varphi }_{ovl}=2\pi ·\Delta f_{ovl}·\Delta t$$

(40)

where $\Delta f_{ovl}$ is the Doppler frequency difference between the two bursts at the given SAR pixel in the overlapped region between the two bursts (see Figure 19), which can be computed as follows:

$$\Delta f_{ovl}=K_T·\left(t_{a,i+1}-t_{a,i}\right)$$

(41)

where $t_{a,i}$ and $t_{a,i+1}$ are the relative times of the ith and i+1th bursts with respect to the starting times of the bursts. The time misregistration is converted in the azimuth misregistration, in terms of the fraction of a pixel, as $\Delta x=\Delta t/\tau$, where $\tau =1/PRF$ is the time duration of a single azimuth SAR pixel and PRF is the sensor pulse repetition frequency (PRF).

Therefore, a constant azimuth misregistration leads to a spurious phase ramp along the azimuth direction, which turns out into the presence of phase jumps between consecutive burst interferograms if not adequately compensated for. Note that the phase is measured in the $\left[-\pi ,\pi \right[$ interval. For typical values of the TOPS mode as implemented in Sentinel-1 system (see the parameters listed in Table 3), the phase in Equation (40) will become ambiguous when the azimuth misalignment (measured in samples) exceeds

$$\Delta x_{max}=\pm \frac{1}{2\tau ·\Delta f_{ovl}}=\pm \frac{PRF}{2\Delta f_{ovl}}$$

(42)

which corresponds to $\Delta x_{max}=\pm 0.0119$ when considering the parameters listed in Table 3. The (unknown) azimuth misregistration can usually be estimated by processing some contiguous bursts from the wrapped phase histogram using a simple periodogram [81] as

$$\widehat{\Delta }x=\underset{\Delta x}{\mathrm{argmin}}\left\{Re\left({\displaystyle \sum }_{p}exp\left[j\left({\varphi }_{ovl,p}-2\pi ·\frac{\Delta f_{ovl,p}}{PRF}·\Delta x\right)\right]\right)\right\}$$

(43)

where $Re(·)$ is the operator that extracts the real part of a complex number. However, if we assume that the preliminary geometrical coregistration is accurate, the measured phase can be considered not ambiguous [134,135], and simply be estimated as

$$\widehat{\Delta }x=\frac{PRF·{{\displaystyle \sum }}_p{\varphi }_{ovl,p}}{2\pi ·{{\displaystyle \sum }}_p\Delta f_{ovl,p}}$$

(44)

To improve the coherence of the burst overlap interferograms, adaptive azimuth band-filtering is applied to compensate for the different mean Doppler centroids of the bursts and the orbit crossing angles of the two bursts [78].

To show the potential of the ESD approach, let us consider a set of Sentinel-1 SAR data collected over the area of Ridgecrest, California, which was struck by an earthquake in July 2019. The two-burst interferograms elated to the InSAR data pair between 11 and 23 May 2019 (e.g., in the stationary case) are shown in Figure 20A, whereas Figure 20B shows the phase difference between the two contiguous bursts (before application of the ESD-based co-registration procedure) in the overlapped region (highlighted with a red rectangle). In this case, as is evident, the phase difference only accounts for the residual co-registration errors, as there are no significant along-track deformations present.

The performance of ESD has been evaluated using Terra-SARX [78] and Sentinel-1 data [79,136]. The attainable coregistration accuracy is given by

$${\sigma }_{ESD}=\frac{PRF}{2\pi \Delta f_{ovl}}·\sqrt{\frac{3}{NL}}·\frac{\sqrt{1-{\gamma }^2}}{\gamma }$$

(45)

where $\gamma$ is the coherence. Similarly, if the MAI technique is applied to TOPS data, the expected accuracy in the Cramér–Rao limit is [112]

$${\sigma }_{MAI}=\frac{PRF}{2\pi \Delta f_{MAI}}·\sqrt{\frac{3}{NL}}·\frac{\sqrt{1-{\gamma }^2}}{\gamma }$$

(46)

The cross-comparison between the expected accuracies in both cases is shown in the plot of Figure 21. The results suggest that only the targets with high coherence have to be selected to estimate the co-registration mistakes. A method to select and exploit the coherent targets has been proposed in [135].

Within the framework of generating surface displacement maps with multi-temporal InSAR methods [55,56,57,58,59,60,61,62,63], it is crucial to implement proper multi-temporal co-registration procedures that might be efficient for the overall set of InSAR data pairs to be processed. ESD approaches can be successfully applied to estimate the constant misregistration error along the azimuth direction of single SAR data pairs.

If we consider a network composed by M multi-temporal InSAR data-pairs, we can compute the M-length vector of the azimuthal misregistration—namely, $\Delta \mathit{X}=\left[\Delta x_1,\Delta x_2,\dots ,\Delta x_M\right]$—and infer the misregistration errors related to every single SAR image—namely, $\mathit{X}$—by solving the following system of linear equations:

$$\Delta \mathit{X}=\mathit{A}·\mathit{X}$$

(47)

in a least-squares or a weighted-least-squares way, where $\mathit{A}$ is the incidence-like matrix of the network graph. In [81,134,135,137], different families of InSAR networks are shown. Similarly, as demonstrated for the problem in Section 3.5, the accuracy of the estimated misregistration values for every image is finally obtained as

$${\mathit{C}}_X={\left({\mathit{A}}^T·{\mathit{C}}_{\Delta X}^{-1}·\mathit{A}\right)}^{-1}\cong {\sigma }_{\Delta X}^2{\left({\mathit{A}}^T·\mathit{A}\right)}^{-1}$$

(48)

where ${\sigma }_{\Delta x}\approx 0.001$ [86]. Nevertheless, several challenges are still open that require further effort to perform the coregistration of TOPS SAR data pairs when the scenes are not stationary—that is to say, when a significant azimuth displacement occurs in the analyzed case-study area. In that case, the intra-burst phase difference of the overlapped regions also depends on the azimuthal displacement, thus making estimation of the constant misregistration very challenging and inducing co-registration mistakes in the processed SAR data. To get an idea, Figure 22B shows one bursting interferogram related to the SAR data pairs of 11 May 2019 to 10 July 2019 following the 7.1 Mw earthquake that hit the area.

The burst overlap phase difference estimates the azimuthal displacement [81,123]. In these scenarios, constant misregistration values are usually calculated by excluding the zones with large displacements or by assuming knowledge of a model for the azimuthal ground deformation. Another source of non-stationarity azimuth misregistration is the inhomogeneity in ionospheric propagation delay, as discussed in Section 3.7. Figure 20A and Figure 22A show the pre-seismic and the co-seismic interferograms.

3.7. SD Atmospheric Artefacts Retrieval

The inhomogeneity in ionospheric propagation delay causes variable azimuth misregistration that cannot be easily corrected. An important class of phenomena influencing SAR remote sensing measurements concerns the Earth’s atmosphere [138], which is the medium through which the electromagnetic radiation propagates. From the InSAR point of view, the uncompensated atmospheric phase screen (APS) represents a factor of great disturbance, because it introduces delays in the measured InSAR phase that result in additional errors in the obtained DInSAR products. Atmospheric disturbances can be referred typically as the result of one of two sources: (i) tropospheric, due to water vapor absorption of electromagnetic radiation, and (ii) ionospheric, due to interactions of electromagnetic radiation with the total electron content (TEC) in the ionosphere. To this aim, several techniques have been developed in recent decades for the mitigation of atmospheric disturbance effects relative to tropospheric and ionospheric effects [54,138,139]. In particular, the absorption effect of water vapor generates an apparent phase signal that masks deformation signals, especially in areas of high relief where topographical phase components are correlated to the atmospheric phase [140,141,142]. In this context, most of the techniques developed for the correction of the water vapor effects exploit the assumption that atmospheric signals have rapid variations over time (between multiple and temporally different acquisitions of the same scene) and slow variations in space (within a single acquisition) [54,143]. In addition, hybrid methodologies using signal-processing techniques have been developed over the last two decades to attain information about the physical characteristic of water vapor, which are estimated from numerical models of weather forecasting such as those developed by the European Centre for Medium Range Weather Forecast (ECMWF) [140,141,144]. In the case of ionospheric influence, the interactions of electromagnetic waves with electrons contained in the ionosphere can cause path delays on radar signals, particularly with L- and P-band SAR systems, because the ionospheric phase distortion effects are inversely proportional to the frequency of the used signal. Generally, high TEC and low RADAR frequencies cause significant ionospheric artefacts. In this framework, new MAI-based ionospheric correction techniques have very recently been proposed [83,84,85,86] to correct the InSAR phase. These techniques exploit the fact that the azimuth derivative of the ionospheric phase distortion on RADAR interferograms is linearly proportional to the apparent pixel azimuth displacement $\Delta x$:

$$\frac{d{\phi }_{ion}}{dx}=\alpha \frac{4\pi }{\lambda }\Delta x$$

(49)

where $\alpha$ is a system- and geometry-dependent factor. Note that ${\phi }_{ion}$ can be expressed as

$${\phi }_{ion}=-\frac{4\pi K}{cf}·\frac{1}{cos\theta }\Delta \mathrm{TEC}$$

(50)

where K is equal to 40.28 $\left(m^3/s^2\right)$, c is the speed of light, f is the radar operational frequency and $\Delta \mathrm{TEC}$ is the variation of TEC. By using the MAI technique, the estimation of $\Delta x$ can be carried out accurately, as in the study made by Jung et al. [83]. In [83] the authors exploited the potential of the MAI technique to get an accurate estimate and compensation of the ionospheric phase components with a high spatial resolution. Specifically, they applied the method on a pair of ALOS PALSAR (L-band) radar sensor acquisitions affected by strong ionospheric effects. After the generation of the MAI interferograms, several steps were performed to calculate the ionospheric phase components for every azimuth x and range y, as follows:

$${\phi }_{ion}\left(x,y\right)={\displaystyle \sum }_{u=1}^x\left[\left(\alpha ·{\overline{\phi }}_{MAI}\left(u,y\right)+\beta \right)·{\Delta }_{az}\right]+\mathrm{C}\left(\mathrm{y}\right)$$

(51)

where $\alpha$ and $\beta$ are two parameters determined by a polynomial fitting; ${\overline{\phi }}_{MAI}=-\frac{l}{n\lambda }{\phi }_{MAI}$ is the scaled MAI phase; ${\Delta }_{az}$ is the multi-looked azimuth pixel spacing; and C(y) is an integral constant that varies along with the range position. Finally, subtracting the estimated ionospheric phase term from the InSAR algorithm, the corrected InSAR phase ${\widehat{\phi }}_{InSAR}$ is obtained as

$${\widehat{\phi }}_{InSAR}\left(x,y\right)=Wr\left({\phi }_{InSAR}\left(x,y\right)-{\phi }_{ion}\left(x,y\right)\right)$$

(52)

In Figure 7a of [83], Jung et al. show that the performance obtained by their method in correcting InSAR interferograms was affected by ionospheric effects. Furthermore, in Figure 9a,b of [83], Jung et al. show a 2-D Fourier analysis, relative to the power spectra of the uncorrected and corrected interferograms. Very recently, the same authors of [83] have presented a more-efficient ionospheric correction method that improves the performance of the aforementioned work by using a multivariate regression method [84] to correct both ionospheric and orbital artefacts. Finally, Fatthai et al. [86] very recently proposed an algorithm based on the use of the MAI technique for the time-series estimation of ionospheric phase delay from a stack of SAR acquisitions.

4. Generation of Multi-Track 3-D Ground Displacement Time-Series

In this Section, we summarize recent developments in the generation of 3-D ground displacement velocity maps, as well as in the retrieval of 3-D ground deformation time-series. Interested readers can find a more comprehensive overview on this topic in [76]. In our work, the focus is principally on the generation of combined MAI/DInSAR-driven 3-D maps. Moreover, we propose a new algorithm for the generation of 3-D ground displacement time-series, complementing the recently developed minimum acceleration combination (MinA) method [75] (applied for the generation of up–down and east–west deformation components) with the north–south ground displacement time-series estimated by applying the methods shown in Section 3.2.

4.1. Overview of the Techniques for the Generation of 3-D Displacement Time-Series

One of the limits of the conventional DInSAR methods is that only the projection of the deformation along the LOS direction can be estimated. Nevertheless, the availability of DInSAR data products (i.e., mean deformation rate and corresponding displacement time-series) computed from SAR data collected from ascending and descending data tracks allows a simple combination of the mean displacement velocity measurements as estimated from (at least) two complementary viewing orbital geometries. By considering the geometry portrayed in Figure 23, it can be demonstrated that the east–west and up–down deformation rate components can be obtained as [68,76]

$$d_{LOS}^{\left(East\right)}\approx \frac{d_{LOS}^{\left(Desc\right)}-d_{LOS}^{\left(Asc\right)}}{2sin\left(\vartheta \right)}$$

(53a)

$$d_{LOS}^{\left(Up\right)}\approx \frac{d_{LOS}^{\left(Desc\right)}+d_{LOS}^{\left(Asc\right)}}{2cos\left(\vartheta \right)}$$

(53b)

Equation (53a,b) is obtained by applying simple geometrical relations and assuming that ascending/descending measurements are taken by the same sensor, illuminating the imaged scene with approximately the same (absolute) side-looking angle. Clearly, the sum of the ascending/ descending LOS-projected displacement measurements is related to the vertical components of the ground deformation, whereas the difference of the ascending/descending components gives an estimate of the E-W components of the deformation. A more general relation arises when DInSAR measurements are obtained from complementary viewing angles but from different radar platforms, observing the scene with distinctive side-looking angles, namely, ${\theta }_{asc}$ and ${\theta }_{desc}$ (see Figure 16 in [76]):

$$\left(\begin{array}{c}{d}_H\\ d_v\end{array}\right)={\left(\begin{array}{cc}\sin \left({\vartheta }_1\right)& \cos \left({\vartheta }_1\right)\\ -\sin \left({\vartheta }_2\right)& \cos \left({\vartheta }_2\right)\end{array}\right)}^{-1}·\left(\begin{array}{c}{d}_{LO{S}_1}\\ d_{LO{S}_2}\end{array}\right)=\frac{\left(\begin{array}{c}{d}_{LO{S}_1}cos\left({\vartheta }_2\right)-d_{LO{S}_2}cos\left({\vartheta }_1\right)\\ d_{LO{S}_1}sin\left({\vartheta }_2\right)-d_{LO{S}_2}sin\left({\vartheta }_1\right)\end{array}\right)}{sin\left({\vartheta }_1+{\vartheta }_2\right)}$$

(54)

Nonetheless, in the derivation of Equation (54), there is no valuable information related to the Earth’s ground deformation along the north–south direction. Indeed, because modern spaceborne radar systems are mounted on satellites that fly nearly polar orbits [41], the north–south (N-S) components of the deformation cannot be reliably measured using conventional DInSAR methods. Several algorithms have been developed to generate 3-D average displacement maps as well as 3-D displacement time-series by combining multi-satellite/multi-orbit data tracks [65,66,68,72,74]. In our review paper, we mainly concentrate on the MinA technique [75], whose basic rationale is exploited to extend its applicability at the generation of north–south displacement time-series, as presented in Section 4.3. The MinA technique assumes the availability of K independent SAR datasets collected by multiple SAR platforms at the ordered acquisition times $T^{\left(j\right)}=\left[t_0^{\left(j\right)},t_1^{\left(j\right)},\dots ,t_{Qj-1}^{\left(j\right)}\right]$, where $Q_j$ is the number of SAR images of the jth SAR dataset. Every single data-track is separately processed, and the LOS-projected displacement time-series are evaluated in correspondence with the location of coherent targets, which are well detected during each data processing step. After geocoding, the DInSAR products are efficiently combined over the common coherent points of the geocoded grid. The combination relies on the application of the mathematical relationships existing between the measured LOS-projection deformations and their corresponding 3-D components (see Figure 23):

$$d_{Los}=\mathit{d}·{\widehat{\mathit{i}}}_{LOS}=sin\left(\vartheta \right)cos\left(\phi \right)d_{East-West}-cos\left(\vartheta \right)d_{Up-Down}+sin\left(\vartheta \right)sin\left(\phi \right)d_{North-South}$$

(55)

where $\phi$ is the radar heading angle (see Figure 23B). In the simplified case, the dependence of the LOS-projected displacement from the north–south components can be neglected. Accordingly, Equation (55) simplifies as follows:

$$d_{Los}=\mathit{d}·{\widehat{\mathit{i}}}_{LOS}\cong sin\left(\vartheta \right)d_{East-West}-cos\left(\vartheta \right)d_{Up-Down}$$

(56)

where we have assumed that $\phi \cong 0$. The MinA technique relies on the application of Equations (55) and (56) for the generation of 2-D (and 3-D) ground displacement time-series. To this aim, the whole set of chronologically ordered acquisition times—namely, $T={{\displaystyle \cup }}_{j=1}^K{T}^{\left(j\right)}$—is considered, and, for every coherent point, the vectors of the velocities among consecutive time acquisitions (see [55]) for the east–west and up–down components—namely, ${\mathit{V}}_{\mathit{E}\mathit{a}\mathit{s}\mathit{t}–\mathit{W}\mathit{e}\mathit{s}\mathit{t}}$ and ${\mathit{V}}_{\mathit{U}\mathit{p}–\mathit{D}\mathit{o}\mathit{w}\mathit{n}}$—are considered as the unknowns of the problem at hand. Accordingly, the following system of linear equations is built (see [75] for the details):

$$\mathit{B}·\left[\begin{array}{c}{\mathit{V}}_{\mathit{E}\mathit{a}\mathit{s}\mathit{t}–\mathit{W}\mathit{e}\mathit{s}\mathit{t}}\\ {\mathit{V}}_{\mathit{U}\mathit{p}–\mathit{D}\mathit{o}\mathit{w}\mathit{n}}\end{array}\right]=\left[\begin{array}{c}{\mathit{d}}_{\mathit{L}\mathit{O}\mathit{S}}^{\left(\mathbf{1}\right)}\\ \begin{array}{c}{\mathit{d}}_{\mathit{L}\mathit{O}\mathit{S}}^{\left(\mathbf{2}\right)}\\ ⋮\\ {\mathit{d}}_{\mathit{L}\mathit{O}\mathit{S}}^{\left(\mathit{K}\right)}\end{array}\end{array}\right]$$

(57)

where ${\mathit{V}}_{\mathit{E}\mathit{a}\mathit{s}\mathit{t}–\mathit{W}\mathit{e}\mathit{s}\mathit{t}}$ and ${\mathit{V}}_{\mathit{U}\mathit{p}–\mathit{D}\mathit{o}\mathit{w}\mathit{n}}$ are the vectors of the LOS-projected deformation time-series from the K different, complementary viewing geometries. The matrix $\mathit{B}$, which takes into account the temporal relationships of combined displacement time-series, is expressed as:

$$\mathit{B}=\left[\begin{array}{cc}{\mathit{B}}^{\left(1\right)}sin{\vartheta }^{\left(1\right)}cos{\phi }^{\left(1\right)}& -{\mathit{B}}^{\left(1\right)}cos{\vartheta }^{\left(1\right)}\\ ⋮& ⋮\\ {\mathit{B}}^{\left(K\right)}sin{\vartheta }^{\left(K\right)}cos{\phi }^{\left(k\right)}& -{\mathit{B}}^{\left(K\right)}cos{\vartheta }^{\left(K\right)}\end{array}\right]$$

(58)

where ${\mathit{B}}^{\left(j\right)}j=1,2,\dots ,K$ is the jth incidence-like matrix (see [55] for further details) of the linear transformation relating the LOS displacement time-series with velocity deformation rates ${\mathit{V}}_{\mathit{E}}$, ${\mathit{V}}_{\mathit{U}}$.

To solve the problem, Equation (57) is regularized by imposing the condition that the (unknown) 2-D (E-W, U-D) deformation time-series at a given SAR pixel location has a minimum acceleration [75]. This condition is obtained by adding the following equations to Equation (57):

$$\mathit{C}=\{\begin{array}{c}\mathsf{\delta }\left(v_{{\mathrm{E}}_{\mathrm{i}+1}}-v_{{\mathrm{E}}_{\mathrm{i}}}\right)=0\hspace{1em}\mathrm{i}=1,\dots ,\mathrm{Q}-2\\ \mathsf{\delta }\left(v_{{\mathrm{U}}_{\mathrm{i}+1}}-v_{{\mathrm{U}}_{\mathrm{i}}}\right)=0\hspace{1em}\mathrm{i}=1,\dots ,\mathrm{Q}-2\end{array}$$

(59)

where $\mathsf{\delta }$ is a regularization factor. The regularized system of linear equations can be re-written as

$$\left[\begin{array}{c}\mathit{B}\\ \mathit{C}\end{array}\right]·\left[\begin{array}{c}{\mathit{V}}_{\mathit{E}\mathit{a}\mathit{s}\mathit{t}–\mathit{W}\mathit{e}\mathit{s}\mathit{t}}\\ {\mathit{V}}_{\mathit{U}\mathit{p}–\mathit{D}\mathit{o}\mathit{w}\mathit{n}}\end{array}\right]=\left[\begin{array}{c}{\mathit{d}}_{\mathit{L}\mathit{O}\mathit{S}}^{\left(\mathbf{1}\right)}\\ \begin{array}{c}{\mathit{d}}_{\mathit{L}\mathit{O}\mathit{S}}^{\left(\mathbf{2}\right)}\\ ⋮\\ \begin{array}{c}{\mathit{d}}_{\mathit{L}\mathit{O}\mathit{S}}^{\left(\mathit{K}\right)}\\ \mathbf{0}\end{array}\end{array}\end{array}\right]$$

(60)

This system of linear equations is solved in the least-squares (LS) sense as

$$\mathit{V}={\rm Y}·\mathit{D}$$

(61)

where $\mathit{V}$ is the (unknown) vector of the LOS-projected E-W and U-P deformation time-series, ${\rm Y}$ is the generalized inverse of the matrix $\left[\begin{array}{c}\mathit{B}\\ \mathit{C}\end{array}\right]$ and $\mathit{D}$ is the vector of geocoded LOS-displacement time-series relative to the K datasets.

Finally, the east–west and up–down deformation time-series are computed by time integration (pixel by pixel) of the obtained 2-D deformation velocity components. An error budget analysis of the MinA technique has previously been discussed in [145].

A regularized problem that is similar to MinA was proposed within the Multidimensional- SBAS (MSBAS) algorithm [146]. Differently from [75], the combination problem discussed in MSBAS was directly applied to the sequences of unwrapped multiple-track differential SAR interferograms. As opposed to MinA, the MSBAS method relies on searching for a minimum-velocity-norm (MN) solution. Therefore, MSBAS requires the simultaneous inversion of several (a few hundred or more) unwrapped interferograms for the retrieval of 3-D components of deformation.

4.2. Experimental Results

Here, we present some results achieved by combining the LOS-projected ground displacement signals collected after applying the multi-temporal SBAS technique to two sets of ascending/descending SAR data. The data were collected by Sentinel-1 radar instruments of the EU Copernicus Sentinels constellation based on the case-study of the area of Ridgecrest in California. The lists of available SAR data are indicated in in

Table 5. Sentinel-1 ascending SAR dataset. Date format is day/month/year.

Acquisition n.	Date
1	11/05/2019
2	23/05/2019
3	16/06/2019
4	28/06/2019
5	10/07/2019
6	22/07/2019
7	03/08/2019
8	15/08/2019
9	27/08/2019
10	08/09/2019
11	02/10/2019
12	14/10/2019
13	26/10/2019

and

Table 6. Sentinel-1 descending SAR dataset. Date format is day/month/year.

Acquisition n.	Date
1	05/05/2019
2	17/05/2019
3	29/05/2019
4	10/06/2019
5	22/06/2019
6	04/07/2019
7	16/07/2019
8	28/07/2019
9	09/08/2019
10	21/08/2019
11	02/09/2019
12	14/09/2019
13	08/10/2019
14	20/10/2019

for the ascending and descending data tracks, respectively.

Figure 24A,B portrays geocoded LOS-projected mean displacement maps of the area, superimposed on an amplitude image of the investigated area. The geocoding operations have been performed using precise orbital information and a 3-arc-second Shuttle RADAR Topography Mission (SRTM) [147] digital elevation model (DEM).

Subsequently, the MinA combination technique was applied to the two sets of geocoded data, to decompose the LOS-projected deformation measurements along the east–west and up–down directions, as shown in Figure 25A,B, respectively

4.3. Generation of North–South Displacement Maps with MAI Approaches

In this subsection, we finally propose an adaptation of the MinA multi-temporal technique [75] for the generation of north–south displacement time-series using a combination of K independent sets of multiple-satellite/multiple-viewing angles. Similar to previous subsections above, we assume here that SAR data were acquired at ordered times $T^{\left(j\right)}=\left[t_0^{\left(j\right)},t_1^{\left(j\right)},\dots ,t_{Qj-1}^{\left(j\right)}\right]$, where $Q_j$ is the number of the available SAR scene in the given data track. $T={{\displaystyle \cup }}_{j=1}^K{T}^{\left(j\right)}$ is the whole set of ordered acquisition times. Let us also assume that the K sets of data are geocoded over a common grid of points, and $\left\{d_{az}^{\left(1\right)},d_{az}^{\left(2\right)},\dots ,d_{az}^{\left(K\right)}\right\}$ are the K vectors of the azimuthal displacement time-series obtained by applying the MAI approaches described in the previous sections.

By applying simple geometric considerations (see Figure 23B), we can straightforwardly derive that

$$d_{az}=\mathit{d}·{\widehat{\mathit{i}}}_{az}=cos\left(\chi \right)·d_{N-s}-sin\left(\chi \right)·d_{E-W}$$

(62)

where $\chi$ is the inclination angle of the orbit with respect to the N-S direction, which is approximately equal to 6° for the Sentinel-1 case. For each coherent point that is well processed in all of the K independent processing stages, a system of linear equations can be built. In particular, if we consider the velocity deformations as unknowns across contiguous time acquisitions for the east–west and north–south directions, the following system of equations can be written:

$$\mathit{B}·\left[\begin{array}{c}{\mathit{V}}_{\mathit{N}\mathit{o}\mathit{r}\mathit{t}\mathit{h}–\mathit{S}\mathit{o}\mathit{u}\mathit{t}\mathit{h}}\\ {\mathit{V}}_{\mathit{E}\mathit{a}\mathit{s}\mathit{t}–\mathit{W}\mathit{e}\mathit{s}\mathit{t}}\end{array}\right]=\left[\begin{array}{c}{\mathit{d}}_{\mathit{a}\mathit{z}}^{\left(\mathbf{1}\right)}\\ ⋮\\ {\mathit{d}}_{\mathit{a}\mathit{z}}^{\left(\mathit{K}\right)}\end{array}\right]$$

(63)

where ${\mathit{V}}_{\mathit{N}\mathit{o}\mathit{r}\mathit{t}\mathit{h}–\mathit{S}\mathit{o}\mathit{u}\mathit{t}\mathit{h}}$ and ${\mathit{V}}_{\mathit{E}\mathit{e}\mathit{a}\mathit{s}\mathit{t}–\mathit{W}\mathit{e}\mathit{s}\mathit{t}}$ are the vectors of the azimuthal projected deformation time-series from the K different, complementary viewing geometries. The matrix B, similar to the case accounted for in Section 4.1, takes into account the temporal relationships of combined displacement time-series and, in this case (see the imaging geometry of Figure 23B), is expressed as follows:

$$\mathit{B}=\left[\begin{array}{cc}{\mathit{B}}^{\left(\mathbf{1}\right)}cos{\chi }^{\left(1\right)}& -{\mathit{B}}^{\left(\mathbf{1}\right)}sin{\chi }^{\left(1\right)}\\ ⋮& ⋮\\ {\mathit{B}}^{\left(\mathit{K}\right)}cos{\chi }^{\left(K\right)}& -{\mathit{B}}^{\left(\mathit{K}\right)}sin{\chi }^{\left(K\right)}\end{array}\right]$$

(64)

Similar to what has already been discussed for the combination of LOS measurements [75], the system of Equation (63) is regularized as

$$\left[\begin{array}{c}\mathit{B}\\ \mathit{C}\end{array}\right]·\left[\begin{array}{c}{\mathit{V}}_{\mathit{N}\mathit{o}\mathit{r}\mathit{t}\mathit{h}–\mathit{S}\mathit{o}\mathit{u}\mathit{t}\mathit{h}}\\ {\mathit{V}}_{\mathit{E}\mathit{a}\mathit{s}\mathit{t}–\mathit{W}\mathit{e}\mathit{s}\mathit{t}}\end{array}\right]=\left[\begin{array}{c}{\mathit{d}}_{\mathit{a}\mathit{z}}^{\left(\mathbf{1}\right)}\\ \begin{array}{c}{\mathit{d}}_{\mathit{a}\mathit{z}}^{\left(\mathbf{2}\right)}\\ ⋮\\ \begin{array}{c}{\mathit{d}}_{\mathit{a}\mathit{z}}^{\left(\mathit{K}\right)}\\ \mathbf{0}\end{array}\end{array}\end{array}\right]$$

(65)

The estimated velocity vector for the north–south components are time-integrated and the N-S ground displacement is computed. Because the projection along the azimuth direction has a very limited accuracy with respect to w.r.t. to the E-W direction, the computed E-W ground displacement time-series (which may be obtained by combining the MAI measurements) definitely has a worse accuracy with respect to w.r.t. obtained using the LOS-projected displacement measurements. The error budget of the N-S displacement time-series can be derived by extending the analysis provided in [145]. The comprehensive theoretical and quantitative analysis of the error budget for the retrieved 3-D ground displacement time-series would require the processing of several independent SAR datasets, as well as a comparison with external ground deformation measurements. This is clearly outside the scope of the current investigation, but remains a matter for future analysis. Figure 26 shows the results of this novel combination method applied to the sets of ascending/descending SAR data relevant to the area of the Ridgecrest earthquake in California. As a result, we have obtained a map of the north–south mean deformation rate of the area. Finally, we selected two points across the fault, labeled (I) and (II) in Figure 25A,B and Figure 26, and Figure 27 shows the plots of the up–down, east–west and north–south displacement time-series.

Additional remarks on the achieved results are now in order. Concerning the Ridgecrest case-study area, the fringe rate in the near proximity of the fault line that was activated by the July 2019 Mw 7.1 earthquake was very high, as is evident from Figure 22A. This led to loss of coherence and unavoidable phase unwrapping errors. The application of conventional InSAR methods, properly complemented with processing strategies to combine multiple-orbit information, was effective for retrieving the up–down and the east–west ground deformation components. The results shown in Figure 25 prove that most of the deformations were horizontal, with a maximum rate of about 250 cm/year corresponding to an abrupt lateral terrain displacement of about 1 meter. The application of the MAI approach to both ascending and descending Sentinel-1 data tracks also allowed us to obtain an estimate of the north–south displacements according to the generation of relevant deformation time-series (see Figure 26). As shown in Figure 12, we also used one single pair of COSMO-SkyMed SAR images to produce a map of the north–south displacement across the fault line. It should be noted that the magnitude of the retrieved relative N-S deformation across the main seismic event of July 2019 was in general agreement with the deformation time-series of N-S components shown in Figure 26 (with northward and southward deformations of up to 1 meter across the fault). Of course, the accuracy of COSMO-SkyMed and Sentinel-1 MAI-driven measurements were not comparable (see plots shown in Figure 14), as the expected accuracy of MAI measurements with Sentinel-1 data is in the order of 2 cm, with an average coherence value of 0.8. Finally, we would like to remark that most of the experiments were carried out using Sentinel-1 data. The aim was to enhance the strength of spectral diversity approaches for the TOPS mode Sentinel-1 data co-registration, as well as to investigate the potential of the MAI technique applied to TOPS mode SAR data.

5. Conclusions

A review of the existing methods used to study large deformation episodes was carried out based on the application of spectral diversity methods for the generation of combined 3-D ground displacement time-series and mean deformation velocity maps. Specifically, we focused on the description of spectral diversity/multiple aperture interferometry methods, and provided a general overview of the theoretical basis of these techniques, addressing their limits and applications, with a particular eye to the use of these technologies for the processing of data coming from new-generation SAR satellites. Several SD techniques were presented for the co-registration of TOPS SAR Sentinel data, the estimation and the removal of the atmospheric phase screen in sequences of SAR data, as well as for the generation of 3-D ground deformation time-series. In the framework of Sentinel-1 SAR data processing, SD and ESD methods would require new adaptations for the analysis of non-stationary scenes where significant along-track displacement signals are present. Furthermore, a novel method combining MAI- and DInSAR-driven data of large deformations was presented, and the relevant results discussed. A quantitative assessment of the presented 3-D methods requires the processing of several datasets. This is a matter for future investigations.