On the Extension of the Minimum Cost Flow ... - Semantic Scholar
2374
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 9, SEPTEMBER 2006
On the Extension of the Minimum Cost Flow Algorithm for Phase Unwrapping of Multitemporal Differential SAR Interferograms Antonio Pepe and Riccardo Lanari, Senior Member, IEEE
Abstract—In this paper, an extension of the minimum cost flow (MCF) algorithm dealing with a sparse data grid, which allows the unwrapping of multitemporal differential synthetic aperture radar (SAR) interferograms for the generation of deformation time series, is presented. The proposed approach exploits both the spatial characteristics and the temporal relationships among multiple interferograms relevant to a properly chosen sequence. In particular, the presented solution involves two main steps: first of all, for each arc connecting neighboring pixels on the interferometric azimuth/range grid, the unwrapped phase gradients are estimated via the MCF technique applied in the temporal/ perpendicular baseline plane. Following this step, these estimates are used as a starting point for the spatial-unwrapping operation implemented again via the MCF approach but carried out in the azimuth/range plane. The presented results, achieved on simulated and real European Remote Sensing satellite SAR data, confirm the effectiveness of the extended MCF unwrapping algorithm. Index Terms—Differential synthetic aperture radar interferometry (DInSAR), minimum cost flow (MCF), phase unwrapping (PhU).
I. I NTRODUCTION
D
IFFERENTIAL synthetic aperture radar interferometry (DInSAR) is a remote-sensing technique with important applications in the investigation of several geophysical processes due to its capability to produce spatially dense deformation maps with centimeter to millimeter accuracy [1]. While the DInSAR approach has been applied first to the analysis of single deformation episodes [2]–[5], there is a growing interest on extending this technique to the study of the temporal evolution of the detected displacements via the generation of deformation time series. To achieve this task, the information available from each interferometric data pair must be properly related to those included in the other acquisitions via the generation, and a subsequent combination [6]–[12], of an appropriate sequence of DInSAR interferograms. In this context, a critical task is represented by the phase unwrapping
Manuscript received August 17, 2005; revised December 16, 2005. This work was supported in part by the Italian Space Agency and in part by the European Union on Provision 3.16 under the project of the Campania Regional Center of Competence “Analysis and Monitoring of the Environmental Risk.” A. Pepe is with the Dipartimento di Ingegneria Elettronica e delle Telecomunicazioni, Università di Napoli “Federico II,” Napoli 80125, Italy. R. Lanari is with the Istituto per il Rilevamento Elettromagnetico dell’Ambiente (IREA), National Research Council of Italy (CNR), Napoli 80124, Italy (e-mail: [email protected]). Digital Object Identifier 10.1109/TGRS.2006.873207
(PhU) operation, especially if we consider high temporal and/or spatial baseline DInSAR data pairs that are characterized by significant decorrelation phenomena [13]. In the SAR scenario, a rather popular PhU technique is based on the solution of a minimum cost flow (MCF) network [14], and this approach has been subsequently changed in order to deal with sparse data [15]. In this case, the grid of the investigated samples is relevant to the coherent pixels of the DInSAR interferograms, while the Delaunay triangulation can be used to define the neighboring points and the elementary cycles in the set of the identified coherent sparse pixels. The possibility to extend this approach to the three-dimensional case, in order to simultaneously unwrap an interferometric sequence, has been already investigated in [16]. However, in this case, the problem could not be formulated in terms of network MCF procedures. Accordingly, no computationally efficient codes are available, and therefore, the overall unwrapping process can be extremely time consuming, particularly for long interferogram sequences. In this paper, we propose an innovative unwrapping solution, shortly described first in [17], that, in addition to exploiting the spatial characteristics of each DInSAR interferogram, also exploits the temporal relationships among a properly selected interferogram sequence, thus allowing to improve the performance of the original MCF technique [15]. The starting point of our approach, mostly oriented to deformation time series generation, is the computation of two Delaunay triangulations. The former is relevant to the SAR data acquisitions distribution in the so-called “temporal/ perpendicular baseline” (T × B⊥ ) plane and allows us to identify the DInSAR interferogram sequence to be computed; the latter is carried out in the azimuth/range (AZ × RG ) plane and involves the coherent pixels common to the interferograms within the sequence. The unwrapping operation of the overall data set is then performed via a two-step processing procedure: first of all, we identify all the segments of the computed “spatial” triangles (i.e., the arcs connecting the coherent neighboring pixels in the AZ × RG plane), and for each of these, we carry out a “temporal” PhU step by applying the MCF technique to the grid relevant to the T × B⊥ plane. The second step uses, for each arc, the previously computed unwrapped phases as a starting point for the subsequent spatial-unwrapping operation performed on each single interferogram. Again, the standard MCF network-programming technique is applied, and in this case, we use the “costs” of the previous solutions (i.e., achieved in the T × B⊥ plane) to set the weights associated
0196-2892/$20.00 © 2006 IEEE
PEPE AND LANARI: PHASE UNWRAPPING OF MULTITEMPORAL DIFFERENTIAL SAR INTERFEROGRAMS
2375
Fig. 1. SAR data representation in the temporal/perpendicular baseline plane for the ERS SAR data analyzed in the following experiments (see Section IV). (a) SAR image distribution. (b) Delaunay triangulation. (c) Triangulation after removal of triangles with sides characterized by spatial and/or temporal baseline values exceeding the selected thresholds (corresponding in our experiments to 300 m and 1500 days, respectively).
with the single arcs involved in the spatial unwrapping. Overall, the basic rationale of the procedure is quite simple: to exploit the temporal relationships among the computed multitemporal interferograms to bootstrap the subsequent spatial-unwrapping operation. We remark that the presented approach is focused on multilook interferograms. Moreover, the possibility of applying the MCF network-programming algorithm to both the temporal and spatial-unwrapping steps leads to a computationally efficient procedure. In order to validate the presented approach, we have carried out a number of experiments on simulated and real SAR data, the latter involving 58 European Remote Sensing (ERS) satellite descending acquisitions relevant to the Abruzzi area (Italy), spanning the 1992–2002 time interval. The achieved results demonstrate the effectiveness of the proposed PhU approach. II. DI N SAR I NTERFEROGRAM G ENERATION Let us start our analysis by investigating the generation process of the interferograms needed for the algorithm implementation discussed in the following sections. Accordingly, we consider a set of N + 1 independent SAR images of the same area, acquired at the ordered (t0 , t1 , . . . , tN ) times. We also assume that each image is coregistered with respect to a reference one, with respect to which we may compute the temporal
and spatial (perpendicular) baseline components ti − tmaster , i = 0, . . . , N and b⊥i , i = 0, . . . , N , respectively, tmaster being the image reference time. Accordingly, each SAR image can be represented by a point in the T × B⊥ plane [see Fig. 1(a)], where we may also compute a Delaunay triangulation [see Fig. 1(b)]. Note that, in order to generate this triangulation, we need to define a ratio between the perpendicular and temporal baseline axis units. In our analysis, we assume that this ratio is equal to δT /δb⊥ , wherein δT and δb⊥ represent the maximum allowed temporal and spatial baselines, respectively, which have been introduced in order to avoid the generation of interferograms strongly corrupted by decorrelation phenomena [13]. Note that the selection of the δT /δb⊥ factor is not a very critical issue. Indeed, within a realistic range of variation of this ratio, say in the interval between one and five, the corresponding triangulation in the T × B⊥ plane typically does not significantly change, as shown in the analysis presented in the Appendix. Note also that each arc connecting two different points in the T × B⊥ plane, say Pi ≡ (ti , b⊥i ) and Pj ≡ (tj , b⊥j ), identifies a corresponding DInSAR data pair. Therefore, as a result of this triangulation, we finally identify a sequence of interferograms for which ∆t = (∆t0 , ∆t1 , . . . , ∆tM −1 ) and ∆b⊥ = (∆b⊥0 , ∆b⊥1 , . . . , ∆b⊥M −1 ) represent the associated temporal and perpendicular baseline vectors, respectively, while M is the overall number of interferograms. As a final issue, we have defined for each arc of our triangulation (i.e., for each selected
2376
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 9, SEPTEMBER 2006
interferometric data pair) a normalized length expressed as follows: 2 2 b⊥i − b⊥j ti − tj Li,j = + (1) δT δb⊥ wherein δT and δb⊥ are our normalization factors. Despite our constraints on the maximum allowed baseline extensions, we remark that the obtained triangulation may involve arcs relevant to data pairs whose baselines exceed the assumed maxima, thus potentially leading to generate interferograms with strong decorrelations. To avoid these effects without losing our triangulation representation, we may remove all the triangles involving arcs with too-large baselines, as shown in Fig. 1(c). Equivalently, we may also remove the triangles corresponding to interferograms including data pairs with large Doppler centroid [18] differences; this is often the case for interferograms involving ERS data acquired after 2000, i.e., following the gyroscope failure events [19]. Note that this triangle removal step may lead to discarding some acquisitions and/or to the generation of more than one independent subset of triangles, i.e., to a data representation consistent with the one described in the DInSAR small baseline subset (SBAS) procedure [7]. Therefore, the compatibility between this data organization and the one exploited in the SBAS technique is evident. Following the identification of the final triangulation in the T × B⊥ plane, the computation of the DInSAR interferogram sequence, hereafter referred to as ϕ = (ϕ0 , ϕ1 , . . . , ϕM −1 ), is performed. At this stage, we can generate the “mask” of the pixels in the AZ × RG planes that are considered coherent within the generated sequence; this mask can be obtained, for instance, by considering those pixels with an estimated coherence value greater than a selected threshold, which are common to a part or even to the entire interferogram sequence. Accordingly, as a final step, we compute a second Delaunay triangulation, which involves the arcs connecting the neighboring pixels of the computed mask in the AZ × RG plane, as shown in the example of Fig. 2. III. P H U A LGORITHM The presented PhU procedure is based on a two-step processing approach that benefits from the information available from both the T × B⊥ and AZ × RG grids. In particular, the key idea is to first carry out, for each arc connecting neighboring pixels of the grid (see Fig. 2), a “temporal” unwrapping operation, which implies the basic MCF approach. The second step relies on the use of these results as a starting point for the “spatial” unwrapping performed on each single interferogram. The key issues of these two processing steps are described in the following analysis, which is focused on the use of multilook interferograms. In particular, this section is organized as follows: first of all, we briefly describe the main aspects of the original MCF approach. Subsequently, we address the main characteristics of the temporal and spatial PhU steps, respectively. A. Basic MCF Approach Although an extensive analysis on the basic principles of the MCF PhU procedure can be found in [14], let us briefly
Fig. 2. Delaunay triangulation in the azimuth/range plane involving the set of spatially coherent pixels (in black) relevant to the SAR data analyzed in the following experiments presented in Section IV.
summarize the main characteristics of this algorithm, in order to better clarify the key concepts of the proposed PhU approach. For the sake of simplicity, we refer to a generic DInSAR pair and assume that ϕ˜ and ψ˜ represent the wrapped and unwrapped phase interferograms, respectively. The MCF PhU method benefits from the relationships existing between the phase differences of pixel pairs relevant to the wrapped and unwrapped signals. In particular, if we consider the spatial arc connecting the generic A and B pixel pair in the AZ × RG plane, the unknown unwrapped phase difference ∆ψ˜AB can be expressed as follows: ˜ ˜ ∆ψ˜AB = ψ(A) − ψ(B) = ϕ(A) ˜ − ϕ(B) ˜ −π,π + 2πHAB = ∆ϕ˜AB + 2πHAB
(2)
wherein the symbol ·−π,π represents the modulo-2π operation, ∆ϕ˜AB is the phase difference computed from the wrapped data, and HAB is the unknown integer number we want to estimate. Therefore, by referring to the MCF unwrapping approach relevant to a sparse data grid [15], we may compute a Delaunay triangulation in the AZ × RG plane (see, for instance, the example shown in Fig. 2) in order to define a set of elementary cycles relevant to the coherent pixels only, and we may impose the irrotational property for the phase gradient, in a discrete ˜ Referring to a generic triangle in the AZ × RG space, of ψ. plane (whose arcs are labeled as α, β, and γ, respectively), this is equivalent to imposing the following constraint: ∆ψ˜α + ∆ψ˜β + ∆ψ˜γ = 0
(3)
PEPE AND LANARI: PHASE UNWRAPPING OF MULTITEMPORAL DIFFERENTIAL SAR INTERFEROGRAMS
which can be also expressed, by taking into account of (2), as follows: ∆ϕ˜α + ∆ϕ˜β + ∆ϕ˜γ Hα + H β + H γ = − . (4) 2π At this stage, the PhU problem can be solved in a very efficient way by recognizing that a network structure underlines it and by searching for the network MCF solution [14]. This can be carried out by choosing the weighted L1 norm for the error criterion, via the following minimization problem: N −1 A wp · |Hp | (5) min Hp
p=0
subject to the constraints (4), wherein the min{·} symbol represents the minimization operation, NA represents the overall number of arcs relevant to the Delaunay triangulation, and w = (w0 , . . . , wp , . . . , wNA −1 ) sets the weights relevant to the weighted L1 norm. B. Temporal Unwrapping Let us start our discussion by reconsidering the spatial Delaunay triangulation in the AZ × RG plane (see, for example, the one in Fig. 2) and to assume that, for each given arc, δϕ = (δϕ0 , δϕ1 , . . . , δϕM −1 ) and δψ = (δψ0 , δψ1 , . . . , δψM −1 ) are the (measured) wrapped and (unknown) unwrapped DInSAR phase-gradient vectors of the DInSAR interferometric sequence, respectively. Moreover, in order to facilitate the procedure, we observe that a profitable model of the unknown unwrapped phase gradient vectors can be considered as follows [6], [7]: ∆b⊥ 4π · · ∆z + ∆t · ∆v (6) m(∆z, ∆v) ≈ λ r · sin(ϑ) wherein ∆z accounts for the error in the knowledge of the scene topography while ∆v accounts for the deformation velocity variations along the considered spatial arc. Moreover, r represents the sensor target distance, λ is the transmitted signal central wavelength, and ϑ is the incidence angle. Note that, because of the spatially correlated behavior of the atmospheric artifacts [6], we have assumed in (6) that they do not contribute to the model m(·). Since the wrapped phases differ from the unwrapped ones by 2π-integer multiples only, the DInSAR unwrapped vector can be expressed as follows: δψ = m + δϕ − m−π,π + 2πK = δχ + 2πK
(7)
wherein K is the 2π-integer multiple vector and δχ = m + δϕ − m−π,π
(8)
represents the phase component, for each considered arc, related to the computed model m(·) and to the measured phase contribution δϕ. Note also that, for the sake of simplicity, the explicit dependence of m(·) on the ∆z and ∆v factors has been neglected in (7) and (8); this notation will be maintained hereafter.
2377
At this stage, we can evaluate for each (∆z, ∆v) pair the unknown vector K in (7) by applying the basic MCF technique (summarized in the previous section) to the T × B⊥ grid. In particular, we search in this case for the solution of the following minimization problem: M −1 min C = |kj | (9) kj j=0
subject to the constraints
δχα + δχβ + δχγ kα + kβ + kγ = −round 2π
(10)
wherein kj is the jth element of the K vector in (7). Moreover, the round[·] symbol represents the operation of approximation to the nearest integer number, and, consistent with (3) and (4) in Section III-A, the constraints in (10) are expressed in terms of a generic triangle in the T × B⊥ plane whose arcs have been labeled as α, β, and γ, respectively. We remark that for each (∆z, ∆v) pair, we will have a different solution, which is characterized by its overall network cost C=
M −1
|kj |.
(11)
j=0
Accordingly, by exploring different values of the (∆z, ∆v) pair and the corresponding model m in (6), we may finally identify the unknown vector δψ opt = δχopt + 2πKopt
(12)
characterized by the “overall minimum” cost Cmin =
M −1
koptj
(13)
j=0
wherein δχopt and Kopt are the vectors (where koptj and δχoptj are the jth components, respectively) defined in (7) and (8) and evaluated in the minimum-cost condition. Accordingly, δψ opt in (12) will represent our estimate of the unwrapped vector for the considered arc. C. Spatial Unwrapping With regard to this operation, it exploits the unwrapped DInSAR phase differences computed via the previous temporal PhU step [see (12)] as a starting point for the spatialunwrapping operation. This second step is carried out on the single interferograms through the application of the MCF unwrapping technique [15] summarized in Section III-A. Moreover, in this case, we may use the minimum costs of the solutions achieved in the T × B⊥ plane, in order to set the weights wp [see (5)] associated to the single arcs. In particular, since the weights represent the confidence on the correctness of the achieved solution, we can reasonably assume that, for each spatial arc, the smaller is the temporal network cost, the larger
2378
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 9, SEPTEMBER 2006
TABLE I ERS 1/2 SAR DATA SET
Fig. 3. SAR products relevant to the investigated area and represented in the azimuth/range plane. (a) Multilook image. (b) Mask of the spatially coherent pixels (in black) defining the investigated spatial grid. (c) Corresponding spatial Delaunay triangulation. The latter is the same as that in Fig. 2 but it is repeated for the sake of completeness.
should be the weight applied within the minimization step [see (5)] involved in the spatial unwrapping. Accordingly, an inverse exponential relationship between the temporal network cost and the corresponding spatial weight wp has been chosen wp =
2S 2Cmin
1,
, Cmin < ρ Cmin ≥ ρ
(14)
wherein 2S represents the maximum allowed weight, Cmin is the minimum temporal cost defined in (13), with Cmin < S, and ρ is a threshold value that, based on experimental results, is typically set not greater than 5% of the total number of interferograms M . Following the implementation of this spatial-unwrapping step, which is carried out for each interferogram, the unwrapped phase can be easily retrieved via the integration of the computed spatial gradients. IV. R ESULTS In order to investigate the performance of the proposed unwrapping approach, a number of experiments, involving both real and simulated data, were carried out. The considered real data set is relevant to the Central Apennines (Abruzzi, Italy) area [see the multilook image shown in Fig. 3(a)], and is composed of 58 acquisitions collected by the ERS 1/2 sensors between August 1992 and October 2002 (see Table I), which are relevant to track 308 and frame 2755. The space/time data distribution is the one depicted in Fig. 1(a), while the final triangulation is shown in Fig. 1(c); the latter has been obtained by imposing a maximum spatial and temporal baseline value of 300 m and 1500 days, respectively, and a Doppler centroid separation, between the image pair, not larger than 1000 Hz. We further remark that, instead of imposing separate constraints on the maximum allowed spatial and temporal baselines, we could alternatively consider an upper limit for the normalized arc length defined in (1). Moreover, the value of the imposed constraints could be also changed depending on the characteristics of the investigated zone. For instance, larger values could be considered in the presence of highly coherent areas, as in the case of urbanized zones.
Following the triangulation step shown in Fig. 1(c), we have identified the sequence of DInSAR interferograms to be generated for implementing the proposed algorithm. In particular, we have computed 145 interferograms with four looks in the range direction and 20 looks in the azimuth one, having a pixel spacing of about 90 × 90 m. Note that the number of interferograms generated from each SAR image is dependent on the result of the Delaunay triangulation within the T × B⊥ plane (see Fig. 1). In this case, it is evident that acquisitions
PEPE AND LANARI: PHASE UNWRAPPING OF MULTITEMPORAL DIFFERENTIAL SAR INTERFEROGRAMS
2379
Fig. 4. Number of occurrences for each SAR image within the generated interferogram data set.
Fig. 6. Examples of the computed DInSAR interferograms and of the corresponding coherence maps. (a) and (d) have been computed from the ERS data pair acquired on July 15, 1997 and September 23, 1997, respectively (|∆b⊥ | = 39 m). (b) and (e) have been computed from the ERS data pair acquired on April 6, 1999 and September 28, 1999, respectively (|∆b⊥ | = 138 m). (c) and (f) have been computed from the ERS data pair acquired on October 28, 1997 and January 11, 2000, respectively (|∆b⊥ | = 207 m).
Fig. 5. Distribution of the temporal and spatial baselines relevant to the generated interferogram data set.
relevant to the triangulation boundaries are characterized by fewer arcs, i.e., are involved in a smaller number of interferograms; obviously, this effect becomes less significant when increasing the number of SAR acquisitions. With regard to our data set, we present in Fig. 4 the plot showing the number of appearances of each SAR image within the interferogram data set. Note that only two acquisitions contribute to two interferograms only. Moreover, we also present in Fig. 5 the distribution of the temporal and spatial baselines obtained by using our interferogram selection based on the Delaunay triangulation. We remark that, within the limits we defined as maximum baselines, the interferogram distribution is rather uniform, which is a relevant issue within the DInSAR scenarios. Based on the computed interferograms, we have generated the coherence mask shown in Fig. 3(b), i.e., the mask of the pixels that are considered coherent in our sequence of multilook interferograms. In our case, they have been selected by identifying the pixels that have a coherence, estimated in a box of four pixels in the range direction and 20 pixels in the azimuth one, greater than 0.35 in at least the 30% of the interferograms. Based on the mask shown in Fig. 3(b), the spatial triangulation has been generated [see Fig. 3(c)]. We remark that the selection of this Abruzzi test-site area is related to the difficulty in the unwrapping operation of the relevant interferograms caused
by the presence of strong decorrelation phenomena. To give an idea of these effects, we present in Fig. 6 a selection of computed interferograms. It is evident that the impact of the decorrelation effects can be very strong. Before discussing in detail the results achieved on the considered real SAR data, we prefer to start our analysis with some simulations that represent controlled experiments, allowing a direct assessment of the algorithm performance. To this end, we have retained both the AZ × RG and T × B⊥ triangulations shown in Figs. 1(c) and 3(c), respectively. Moreover, the same interferogram sequence generated for the real data set has been produced in the simulated experiments, but by exploiting the high-pass (HP) spatial components of a digital elevation model (DEM) of the area and by introducing a deformation parametric model. In particular, the jth interferogram of the simulated DInSAR sequence, say ϕsj , has the following expression: ϕsj
∆b⊥j 4π 1 2 · · h + ∆tj · vel + · ∆tj · acc = λ r · sin(ϑ) 2 j = 0, . . . , M − 1
(15)
wherein the first term represents the topographic phase contribution depending on the DEM height h, while the vel and acc factors are the velocity and acceleration model terms, respectively, the latter representing a significant deviation from the assumed model in (6). Moreover, in our experiments, we also assumed that both the simulated deformation components had a Gaussian spatial shape, as shown in Fig. 7, wherein both the masked and unmasked data are presented. In order to investigate the algorithm performance on the simulated data, we unwrapped the generated interferogram
2380
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 9, SEPTEMBER 2006
Fig. 7. Maps relevant to the three components of the model presented in (15) and used to generated the simulated DInSAR interferograms. (a) and (d) Unmasked and masked topography maps, respectively (hmin = −327 m, hmax = 428 m). (b) and (e) Unmasked and masked velocity patterns, respectively (velmin = −2.83 cm/year, velmax = −0.07 cm/year). (c) and (f) Unmasked and masked acceleration patterns, respectively (accmin = −14.15 cm/year2 , accmax = −0.06 cm/year2 ). Note that the position of the minima of the simulated velocity and acceleration patterns are relatively shifted by 100 pixels in both directions.
sequence both via the original and the extended MCF techniques; obviously, the former has been independently applied to each single interferogram. In the performed experiments, a comparison between the original and the unwrapped phases is directly possible. As a result of our analysis, we present in Fig. 8(a)–(c) the percentage, for each interferogram, of correctly unwrapped pixels of the mask shown in Fig. 3(b), obtained by applying the original MCF technique and considering on the horizontal axis the spatial and temporal baselines and the normalized arc length, respectively. Moreover, the results achieved via the extended MCF procedure are shown in Fig. 8(d), where the horizontal axis is relevant to the normalized arc length. The results presented in Fig. 8(a)–(c) clearly show that, for the original MCF approach, an increase of the spatial and/or of the temporal baseline (thus, of the normalized arc length) may lead to significant degradations of the achieved reconstruction. On the contrary, the extended MCF approach guarantees a very good retrieval, as shown in Fig. 8(d). With regard to the computational aspects, we remark that, on a personal computer equipped with an AMD Athlon MP2600+ processor, the basic MCF unwrapping procedure required about 75 min for unwrapping the overall interferogram sequence, while the extended MCF algorithm needed about 225 min. We remark that this computational increase is mostly related to the evaluation of the model function m within the temporalunwrapping step (see Section III-B). Indeed, this operation is carried out via an exhaustive search, with a limited spacing, of the ∆z and ∆v factors within the intervals (−100, 100) m and (−5, 5) cm/year, respectively. Obviously, the availability of an a priori information on the ∆z and ∆v terms would lead to
Fig. 8. Percentage of the coherent pixels of each interferogram that have been correctly unwrapped. (a) Results obtained by applying the original MCF approach with the horizontal axis relevant to the perpendicular baseline. (b) Same as (a) but with the horizontal axis relevant to the temporal baseline. (c) Same as (a) but with the horizontal axis relevant to the normalized arc length L defined in (1). (d) Results obtained by applying the extended MCF approach with the horizontal axis relevant to the normalized arc length L. (e) Results obtained by applying the extended MCF approach but without the model adjustment step; the horizontal axis is relevant to the normalized arc length L.
a significant computational improvement of the extended MCF PhU algorithm performance. As a final test on our simulated data, we investigated the quality of the retrieved solution if the model adjustment step, relevant to m in (6), is not implemented. In this case, the
PEPE AND LANARI: PHASE UNWRAPPING OF MULTITEMPORAL DIFFERENTIAL SAR INTERFEROGRAMS
2381
Fig. 9. Temporal coherence maps obtained from (a) the original and (b) the extended MCF unwrapping procedures, respectively.
required computing time reduced to about 80 min. However, the achieved results presented in Fig. 8(e) clearly show that the computational reduction is “paid” with a considerable degradation of the obtained unwrapped data with respect to those generated through the extended MCF approach, with the latter shown in Fig. 8(d). Let us now move to the real data analysis. In this case, the dimensions of the processed data set are the same as the simulated signals, thus the computing requirements remained practically unchanged; on the contrary, we now had to define a quality index for our results. To achieve this task, we decided to use the SBAS algorithm [7] that, as previously stated, is fully compatible with the presented data representation and allows us to generate DInSAR deformation time series from a sequence of unwrapped multitemporal interferograms. In particular, we considered the following strategy: first of all, we applied the SBAS procedure to both the unwrapped sequences computed via the original and the extended MCF algorithm; this allowed us to produce deformation time series for each pixel of the mask shown in Fig. 3(b). Subsequently, we used the computed time series (without any atmospheric filtering) and the estimated topographic errors to regenerate the original interferograms referred to as ϕ = (ϕ0 , ϕ1 , . . . , ϕM −1 ). At this stage, we use as quality index of the PhU retrieval the temporal coherence factor defined for each pixel as follows: M −1 exp j(ϕp − ϕp ) p=0 , 0 ≤ γ ≤ 1. (16) γ= M Note that for pixels where γ → 1, we expect that no unwrapping errors are present, since a nearly perfect retrieval of the original phase has been obtained. On the other hand, low values of γ will correspond to poorly unwrapped data. Following the application of the SBAS algorithm to the two unwrapped DInSAR sequences, obtained through the original and the extended MCF techniques, respectively, we have evaluated via (16) the temporal coherence map for both data sets. The achieved results are shown in Fig. 9(a) and (b), respectively. The improvement of the coherence obtained through the application
Fig. 10. Histograms of the temporal coherence maps. (a) The dashed line is related to the results shown in Fig. 9(a), which are obtained via the basic MCF approach. The continuous line is relevant to those achieved by using the extended MCF approach, which are presented in Fig. 9(b). (b) Same plots as in (a) with the superimposed diamond plot relevant to the results of the extended MCF approach implemented without 2π multiple corrections, i.e., by assuming Kopt = 0 in (12).
of the extended PhU algorithm is evident, particularly in the areas in the lower left and upper right corners of the map. To emphasize the achieved gain, the histograms relevant to these two results have been computed and superimposed in Fig. 10(a); again, the obtained improvement is clear. In particular, we remark that, in this case, we achieved an increase of about 60% with regard to the number of pixels with a temporal coherence value greater than the selected threshold γ = 0.7, which is a typical value in DInSAR applications [20], [21]. As a final point, we have investigated the impact, on the unwrapping algorithm performance, of the Kopt vector component in (12). To this end, we have repeated the overall processing operation based on the extended MCF approach but assuming no 2π integer correction, i.e., with Kopt = 0. This implies that, based on (12), we assumed δψ opt = δχopt . The histogram of the computed new temporal coherence has been superimposed [see Fig. 10(b)] to those shown in Fig. 10(a). The degradation of the achieved results is remarkable with respect to those of the overall extended MCF procedure. In particular, in this case, we had a reduction of the number of coherent pixels, with γ ≥ 0.7, of about 30% with respect to the extended MCF results. This clearly indicates the significant impact of the 2π integer correction component within the proposed unwrapping procedure.
2382
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 9, SEPTEMBER 2006
As a final remark, we underline that the overall analysis has been focused on multilook interferograms. However, the extension of this technique to the unwrapping of full-resolution DInSAR interferograms is under development. In this case, the noise phenomena affecting the single-look interferograms become extremely relevant; moreover, at the variance of the multilook case, we cannot benefit from the availability of an estimated coherence map generated via spatial averaging. Accordingly, an appropriate selection of temporally coherent pixels on the full-resolution interferograms is a crucial issue. We are presently investigating a solution based on the use of both multilook and single-look interferograms, which follows the lines of the full-resolution SBAS approach presented in [22]; this analysis can be a topic for future discussions. A PPENDIX
Fig. 11. Delaunay triangulations in the temporal/perpendicular baseline plane relevant to the ERS SAR data analyzed in the experiments of Section IV. The triangulations are obtained by assuming that the ratio between the perpendicular and temporal baseline axis units is equal to δT /δb⊥ , with (a) δT = 1500 days and δb⊥ = 300 m, (b) δT = 900 days and δb⊥ = 300 m, and (c) δT = 300 days and δb⊥ = 300 m.
V. C ONCLUSION We have proposed a solution for extending the MCF PhU algorithm dealing with a sparse data grid [15], to process multitemporal DInSAR interferograms for the generation of deformation time series. The approach involves the computation of a properly chosen interferogram sequence and is based on the cascade of two main steps, both involving the use of the basic MCF technique. In particular, we identify first the arcs connecting neighboring coherent pixels of the azimuth/range grid, and for each of these, we estimate the unwrapped phase gradients via the MCF technique applied in the temporal/perpendicular baseline plane. These results are used to bootstrap the subsequent spatial-unwrapping operation carried out on each single interferogram via the conventional MCF approach. The presented results obtained on simulated and real data confirm the effectiveness of the approach. Moreover, we underline that the proposed solution is naturally compatible with the SBAS algorithm for the generation of deformation time series. The implemented extension of the unwrapping procedure is based on the application of the network-programming algorithm for both the temporal- and spatial-unwrapping steps. Accordingly, the overall approach is computationally efficient.
Let us investigate, in the following, the impact of the different selection of the δT /δb⊥ factor, which represents the ratio between the temporal and perpendicular baseline axis units, which is needed in order to generate the triangulation in the T × B⊥ plane (see Section II). To achieve this task, we started from the triangulation relevant to our experiments presented in Section IV, where we have assumed δT = 1500 days and δb⊥ = 300 m. Subsequently, we maintained constant the value of δb⊥ and considered the two additional cases relevant to δT = 900 days and δT = 300 days, respectively. The results of the corresponding triangulations, presented in Fig. 11(a)–(c), respectively, show a strong similarity. More specifically, we found that more than 60% of the arcs (thus, of the corresponding interferograms) are common to the three triangulations. Accordingly, we may finally conclude that, within a realistic range of variation of the δT /δb⊥ ratio, the corresponding triangulations in the T × B⊥ plane do not significantly change. ACKNOWLEDGMENT Precise satellite orbits were provided by the University of Delft, The Netherlands. A Shuttle Radar Topography Mission DEM of the investigated zone has been used for the interferogram generation. The European Space Agency provided the ERS data relevant to the Abruzzi site through the CAT-1 proposal N. 313. The authors would like to thank the reviewers for their constructive remarks and suggestions. R EFERENCES [1] A. K. Gabriel, R. M. Goldstein, and H. A. Zebker, “Mapping small elevation changes over large areas: Differential interferometry,” J. Geophys. Res., vol. 94, no. B7, pp. 9183–9191, Jul. 1989. [2] D. Massonnet, M. Rossi, C. Carmona, F. Ardagna, G. Peltzer, K. Feigl, and T. Rabaute, “The displacement field of the Landers earthquake mapped by radar interferometry,” Nature, vol. 364, no. 6433, pp. 138– 142, Jul. 1993. [3] G. Peltzer and P. A. Rosen, “Surface displacement of the 17 Eureka valley, California, earthquake observed by SAR interferometry,” Science, vol. 268, no. 5215, pp. 1333–1336, Jun. 1995. [4] D. Massonnet, P. Briole, and A. Arnaud, “Deflation of Mount Etna monitored by spaceborne radar interferometry,” Nature, vol. 375, no. 6532, pp. 567–570, Jun. 1995. [5] E. Rignot, “Fast recession of a west Antarctic glacier,” Science, vol. 281, no. 5376, pp. 549–551, Jul. 1998.
PEPE AND LANARI: PHASE UNWRAPPING OF MULTITEMPORAL DIFFERENTIAL SAR INTERFEROGRAMS
[6] A. Ferretti, C. Prati, and F. Rocca, “Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry,” IEEE Trans. Geosci. Remote Sens., vol. 38, no. 5, pp. 2202–2212, Sep. 2000. [7] P. Berardino, G. Fornaro, R. Lanari, and E. Sansosti, “A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms,” IEEE Trans. Geosci. Remote Sens., vol. 40, no. 11, pp. 2375–2383, Nov. 2002. [8] O. Mora, J. J. Mallorquí, and A. Broquetas, “Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images,” IEEE Trans. Geosci. Remote Sens., vol. 41, no. 10, pp. 2243–2253, Oct. 2003. [9] S. Usai, “A least squares database approach for SAR interferometric data,” IEEE Trans. Geosci. Remote Sens., vol. 41, no. 4, pp. 753–760, Apr. 2003. [10] C. Werner, U. Wegmuller, T. Strozzi, and A. Wiesmann, “Interferometric point target analysis for deformation mapping,” in Proc. IGARSS, Toulouse, France, Jul. 2003, pp. 4362–4364. [11] A. Hooper, H. Zebker, P. Segall, and B. Kampes, “A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers,” Geophys. Res. Lett., vol. 31, no. 23, Dec. 2004, DOI: 10.1029/2004GL021737, CD-ROM. [12] M. Crosetto, B. Crippa, and E. Biescas, “Early detection and in-depth analysis of deformation phenomena by radar interferometry,” Eng. Geol., vol. 79, no. 1/2, pp. 81–91, Jun. 2005. [13] H. A. Zebker and J. Villasenor, “Decorrelation in interferometric radar echoes,” IEEE Trans. Geosci. Remote Sens., vol. 30, no. 5, pp. 950–959, Sep. 1992. [14] M. Costantini, “A novel phase unwrapping method based on network programming,” IEEE Trans. Geosci. Remote Sens., vol. 36, no. 3, pp. 813– 821, May 1998. [15] M. Costantini and P. A. Rosen, “A generalized phase unwrapping approach for sparse data,” in Proc. IGARSS, Hamburg, Germany, Jun. 1999, pp. 267–269. [16] M. Costantini, F. Malvarosa, F. Minati, L. Pietranera, and G. Milillo, “A three-dimensional phase unwrapping algorithm for processing of multitemporal SAR interferometric measurements,” in Proc. IGARSS, Toronto, ON, Canada, Jun. 2002, pp. 1741–1743. [17] A. Pepe and R. Lanari, “A space-time minimum cost flow phase unwrapping algorithm for the generation of DInSAR deformation time-series,” in Proc. IGARSS, Seoul, Korea, Jul. 2005, pp. 1979–1982. [18] G. Franceschetti and R. Lanari, Synthetic Aperture Radar Processing. Boca Raton, FL: CRC, Mar. 1999. [19] N. Miranda, B. Rosich, C. Santella, and M. Grion, “Review of the impact of ERS-2 piloting modes on the SAR Doppler stability,” in Proc. Fringe, Frascati, Italy, Dec. 2003, CD-ROM. [20] R. Lanari, P. Lundgren, M. Manzo, and F. Casu, “Satellite radar interferometry time series analysis of surface deformation for Los Angeles, California,” Geophys. Res. Lett., vol. 31, no. 23, Dec. 2004, DOI: 10.1029/2004GL021294, CD-ROM. [21] A. Borgia, P. Tizzani, G. Solaro, M. Manzo, F. Casu, G. Luongo, A. Pepe, P. Berardino, G. Fornaro, E. Sansosti, G. P. Ricciardi, N. Fusi, G. Di Donna, and R. Lanari, “Volcanic spreading of Vesuvius, a new paradigm for interpreting its volcanic activity,” Geophys. Res. Lett., vol. 32, no. 3, Feb. 2005, DOI: 10.1029/2004GL022155, CD-ROM. [22] R. Lanari, O. Mora, M. Manunta, J. J. Mallorquí, P. Berardino, and E. Sansosti, “A small baseline approach for investigating deformations on full resolution differential SAR interferograms,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 7, pp. 1377–1386, Jul. 2004.
2383
Antonio Pepe was born in Salerno, Italy, in 1974. He received the laurea degree in electronic engineering from the University of Napoli “Federico II,” Naples, Italy, in 2000, with a thesis on the synthetic aperture radar (SAR) image processing. He is currently working toward the Ph.D. degree at the same institution. He joined the Istituto per il Rilevamento Elettromagnetico dell’Ambiente (IREA), Institute of the Italian National Research Council (CNR), Naples, in 2002. Since 2003, he has been with the Department of Electronic and Telecommunication Engineering, University of Napoli. In 2005, he was a Visiting Research at the Aerospace Engineering and Engineering Mechanics, The University of Texas at Austin. His research interests concern the differential SAR interferometry applications for the monitoring of surface displacements, such as those produced by subsidence, volcano activity, and earthquakes.
Riccardo Lanari (M’91–SM’01) was born in 1964. He received the degree in electronic engineering (summa cum laude) from the University of Napoli, “Federico II,” Naples, Italy, in 1989. After graduation, he joined the Istituto di Ricerca per l’Elettromagnetismo e i Componenti Elletronici (IRECE), now the Istituto per il Rilevamento Elettromagnetico dell’Ambiente (IREA), a Research Institute of the Italian National Research Council (CNR), where he currently occupies the position of Senior Researcher. He has been a Visiting Scientist at different foreign research institutes. He was a Research Fellow at the Institute of Space and Astronautical Science, Japan in 1993; a Visiting Scientist at the German Aerospace Research Establishment (DLR), in 1991 and 1994, and at the Jet Propulsion Laboratory, in 1997. He has lectured in several national and foreign universities and research centers and is currently a Lecturer of the SAR module course of the International Master in Airborne Photogrammetry and Remote Sensing, offered by the Institute of Geomatics in Barcelona, Spain. He is the holder of two patents on SAR data-processing techniques. His main research interests are in the synthetic aperture radar data-processing field as well as in SAR interferometry techniques. On these fields, he has authored 40 international journal papers, and, in 1999, he authored a book entitled Synthetic Aperture Radar Processing (CRC Press, 1999). Mr. Lanari received a NASA recognition for the development of the ScanSAR processor used for the SRTM mission. He has been invited as Chairman and Cochairman at several international conferences and as a member of the Technical Program Committee for the IGARSS’2000 symposia.
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 9, SEPTEMBER 2006
On the Extension of the Minimum Cost Flow Algorithm for Phase Unwrapping of Multitemporal Differential SAR Interferograms Antonio Pepe and Riccardo Lanari, Senior Member, IEEE
Abstract—In this paper, an extension of the minimum cost flow (MCF) algorithm dealing with a sparse data grid, which allows the unwrapping of multitemporal differential synthetic aperture radar (SAR) interferograms for the generation of deformation time series, is presented. The proposed approach exploits both the spatial characteristics and the temporal relationships among multiple interferograms relevant to a properly chosen sequence. In particular, the presented solution involves two main steps: first of all, for each arc connecting neighboring pixels on the interferometric azimuth/range grid, the unwrapped phase gradients are estimated via the MCF technique applied in the temporal/ perpendicular baseline plane. Following this step, these estimates are used as a starting point for the spatial-unwrapping operation implemented again via the MCF approach but carried out in the azimuth/range plane. The presented results, achieved on simulated and real European Remote Sensing satellite SAR data, confirm the effectiveness of the extended MCF unwrapping algorithm. Index Terms—Differential synthetic aperture radar interferometry (DInSAR), minimum cost flow (MCF), phase unwrapping (PhU).
I. I NTRODUCTION
D
IFFERENTIAL synthetic aperture radar interferometry (DInSAR) is a remote-sensing technique with important applications in the investigation of several geophysical processes due to its capability to produce spatially dense deformation maps with centimeter to millimeter accuracy [1]. While the DInSAR approach has been applied first to the analysis of single deformation episodes [2]–[5], there is a growing interest on extending this technique to the study of the temporal evolution of the detected displacements via the generation of deformation time series. To achieve this task, the information available from each interferometric data pair must be properly related to those included in the other acquisitions via the generation, and a subsequent combination [6]–[12], of an appropriate sequence of DInSAR interferograms. In this context, a critical task is represented by the phase unwrapping
Manuscript received August 17, 2005; revised December 16, 2005. This work was supported in part by the Italian Space Agency and in part by the European Union on Provision 3.16 under the project of the Campania Regional Center of Competence “Analysis and Monitoring of the Environmental Risk.” A. Pepe is with the Dipartimento di Ingegneria Elettronica e delle Telecomunicazioni, Università di Napoli “Federico II,” Napoli 80125, Italy. R. Lanari is with the Istituto per il Rilevamento Elettromagnetico dell’Ambiente (IREA), National Research Council of Italy (CNR), Napoli 80124, Italy (e-mail: [email protected]). Digital Object Identifier 10.1109/TGRS.2006.873207
(PhU) operation, especially if we consider high temporal and/or spatial baseline DInSAR data pairs that are characterized by significant decorrelation phenomena [13]. In the SAR scenario, a rather popular PhU technique is based on the solution of a minimum cost flow (MCF) network [14], and this approach has been subsequently changed in order to deal with sparse data [15]. In this case, the grid of the investigated samples is relevant to the coherent pixels of the DInSAR interferograms, while the Delaunay triangulation can be used to define the neighboring points and the elementary cycles in the set of the identified coherent sparse pixels. The possibility to extend this approach to the three-dimensional case, in order to simultaneously unwrap an interferometric sequence, has been already investigated in [16]. However, in this case, the problem could not be formulated in terms of network MCF procedures. Accordingly, no computationally efficient codes are available, and therefore, the overall unwrapping process can be extremely time consuming, particularly for long interferogram sequences. In this paper, we propose an innovative unwrapping solution, shortly described first in [17], that, in addition to exploiting the spatial characteristics of each DInSAR interferogram, also exploits the temporal relationships among a properly selected interferogram sequence, thus allowing to improve the performance of the original MCF technique [15]. The starting point of our approach, mostly oriented to deformation time series generation, is the computation of two Delaunay triangulations. The former is relevant to the SAR data acquisitions distribution in the so-called “temporal/ perpendicular baseline” (T × B⊥ ) plane and allows us to identify the DInSAR interferogram sequence to be computed; the latter is carried out in the azimuth/range (AZ × RG ) plane and involves the coherent pixels common to the interferograms within the sequence. The unwrapping operation of the overall data set is then performed via a two-step processing procedure: first of all, we identify all the segments of the computed “spatial” triangles (i.e., the arcs connecting the coherent neighboring pixels in the AZ × RG plane), and for each of these, we carry out a “temporal” PhU step by applying the MCF technique to the grid relevant to the T × B⊥ plane. The second step uses, for each arc, the previously computed unwrapped phases as a starting point for the subsequent spatial-unwrapping operation performed on each single interferogram. Again, the standard MCF network-programming technique is applied, and in this case, we use the “costs” of the previous solutions (i.e., achieved in the T × B⊥ plane) to set the weights associated
0196-2892/$20.00 © 2006 IEEE
PEPE AND LANARI: PHASE UNWRAPPING OF MULTITEMPORAL DIFFERENTIAL SAR INTERFEROGRAMS
2375
Fig. 1. SAR data representation in the temporal/perpendicular baseline plane for the ERS SAR data analyzed in the following experiments (see Section IV). (a) SAR image distribution. (b) Delaunay triangulation. (c) Triangulation after removal of triangles with sides characterized by spatial and/or temporal baseline values exceeding the selected thresholds (corresponding in our experiments to 300 m and 1500 days, respectively).
with the single arcs involved in the spatial unwrapping. Overall, the basic rationale of the procedure is quite simple: to exploit the temporal relationships among the computed multitemporal interferograms to bootstrap the subsequent spatial-unwrapping operation. We remark that the presented approach is focused on multilook interferograms. Moreover, the possibility of applying the MCF network-programming algorithm to both the temporal and spatial-unwrapping steps leads to a computationally efficient procedure. In order to validate the presented approach, we have carried out a number of experiments on simulated and real SAR data, the latter involving 58 European Remote Sensing (ERS) satellite descending acquisitions relevant to the Abruzzi area (Italy), spanning the 1992–2002 time interval. The achieved results demonstrate the effectiveness of the proposed PhU approach. II. DI N SAR I NTERFEROGRAM G ENERATION Let us start our analysis by investigating the generation process of the interferograms needed for the algorithm implementation discussed in the following sections. Accordingly, we consider a set of N + 1 independent SAR images of the same area, acquired at the ordered (t0 , t1 , . . . , tN ) times. We also assume that each image is coregistered with respect to a reference one, with respect to which we may compute the temporal
and spatial (perpendicular) baseline components ti − tmaster , i = 0, . . . , N and b⊥i , i = 0, . . . , N , respectively, tmaster being the image reference time. Accordingly, each SAR image can be represented by a point in the T × B⊥ plane [see Fig. 1(a)], where we may also compute a Delaunay triangulation [see Fig. 1(b)]. Note that, in order to generate this triangulation, we need to define a ratio between the perpendicular and temporal baseline axis units. In our analysis, we assume that this ratio is equal to δT /δb⊥ , wherein δT and δb⊥ represent the maximum allowed temporal and spatial baselines, respectively, which have been introduced in order to avoid the generation of interferograms strongly corrupted by decorrelation phenomena [13]. Note that the selection of the δT /δb⊥ factor is not a very critical issue. Indeed, within a realistic range of variation of this ratio, say in the interval between one and five, the corresponding triangulation in the T × B⊥ plane typically does not significantly change, as shown in the analysis presented in the Appendix. Note also that each arc connecting two different points in the T × B⊥ plane, say Pi ≡ (ti , b⊥i ) and Pj ≡ (tj , b⊥j ), identifies a corresponding DInSAR data pair. Therefore, as a result of this triangulation, we finally identify a sequence of interferograms for which ∆t = (∆t0 , ∆t1 , . . . , ∆tM −1 ) and ∆b⊥ = (∆b⊥0 , ∆b⊥1 , . . . , ∆b⊥M −1 ) represent the associated temporal and perpendicular baseline vectors, respectively, while M is the overall number of interferograms. As a final issue, we have defined for each arc of our triangulation (i.e., for each selected
2376
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 9, SEPTEMBER 2006
interferometric data pair) a normalized length expressed as follows: 2 2 b⊥i − b⊥j ti − tj Li,j = + (1) δT δb⊥ wherein δT and δb⊥ are our normalization factors. Despite our constraints on the maximum allowed baseline extensions, we remark that the obtained triangulation may involve arcs relevant to data pairs whose baselines exceed the assumed maxima, thus potentially leading to generate interferograms with strong decorrelations. To avoid these effects without losing our triangulation representation, we may remove all the triangles involving arcs with too-large baselines, as shown in Fig. 1(c). Equivalently, we may also remove the triangles corresponding to interferograms including data pairs with large Doppler centroid [18] differences; this is often the case for interferograms involving ERS data acquired after 2000, i.e., following the gyroscope failure events [19]. Note that this triangle removal step may lead to discarding some acquisitions and/or to the generation of more than one independent subset of triangles, i.e., to a data representation consistent with the one described in the DInSAR small baseline subset (SBAS) procedure [7]. Therefore, the compatibility between this data organization and the one exploited in the SBAS technique is evident. Following the identification of the final triangulation in the T × B⊥ plane, the computation of the DInSAR interferogram sequence, hereafter referred to as ϕ = (ϕ0 , ϕ1 , . . . , ϕM −1 ), is performed. At this stage, we can generate the “mask” of the pixels in the AZ × RG planes that are considered coherent within the generated sequence; this mask can be obtained, for instance, by considering those pixels with an estimated coherence value greater than a selected threshold, which are common to a part or even to the entire interferogram sequence. Accordingly, as a final step, we compute a second Delaunay triangulation, which involves the arcs connecting the neighboring pixels of the computed mask in the AZ × RG plane, as shown in the example of Fig. 2. III. P H U A LGORITHM The presented PhU procedure is based on a two-step processing approach that benefits from the information available from both the T × B⊥ and AZ × RG grids. In particular, the key idea is to first carry out, for each arc connecting neighboring pixels of the grid (see Fig. 2), a “temporal” unwrapping operation, which implies the basic MCF approach. The second step relies on the use of these results as a starting point for the “spatial” unwrapping performed on each single interferogram. The key issues of these two processing steps are described in the following analysis, which is focused on the use of multilook interferograms. In particular, this section is organized as follows: first of all, we briefly describe the main aspects of the original MCF approach. Subsequently, we address the main characteristics of the temporal and spatial PhU steps, respectively. A. Basic MCF Approach Although an extensive analysis on the basic principles of the MCF PhU procedure can be found in [14], let us briefly
Fig. 2. Delaunay triangulation in the azimuth/range plane involving the set of spatially coherent pixels (in black) relevant to the SAR data analyzed in the following experiments presented in Section IV.
summarize the main characteristics of this algorithm, in order to better clarify the key concepts of the proposed PhU approach. For the sake of simplicity, we refer to a generic DInSAR pair and assume that ϕ˜ and ψ˜ represent the wrapped and unwrapped phase interferograms, respectively. The MCF PhU method benefits from the relationships existing between the phase differences of pixel pairs relevant to the wrapped and unwrapped signals. In particular, if we consider the spatial arc connecting the generic A and B pixel pair in the AZ × RG plane, the unknown unwrapped phase difference ∆ψ˜AB can be expressed as follows: ˜ ˜ ∆ψ˜AB = ψ(A) − ψ(B) = ϕ(A) ˜ − ϕ(B) ˜ −π,π + 2πHAB = ∆ϕ˜AB + 2πHAB
(2)
wherein the symbol ·−π,π represents the modulo-2π operation, ∆ϕ˜AB is the phase difference computed from the wrapped data, and HAB is the unknown integer number we want to estimate. Therefore, by referring to the MCF unwrapping approach relevant to a sparse data grid [15], we may compute a Delaunay triangulation in the AZ × RG plane (see, for instance, the example shown in Fig. 2) in order to define a set of elementary cycles relevant to the coherent pixels only, and we may impose the irrotational property for the phase gradient, in a discrete ˜ Referring to a generic triangle in the AZ × RG space, of ψ. plane (whose arcs are labeled as α, β, and γ, respectively), this is equivalent to imposing the following constraint: ∆ψ˜α + ∆ψ˜β + ∆ψ˜γ = 0
(3)
PEPE AND LANARI: PHASE UNWRAPPING OF MULTITEMPORAL DIFFERENTIAL SAR INTERFEROGRAMS
which can be also expressed, by taking into account of (2), as follows: ∆ϕ˜α + ∆ϕ˜β + ∆ϕ˜γ Hα + H β + H γ = − . (4) 2π At this stage, the PhU problem can be solved in a very efficient way by recognizing that a network structure underlines it and by searching for the network MCF solution [14]. This can be carried out by choosing the weighted L1 norm for the error criterion, via the following minimization problem: N −1 A wp · |Hp | (5) min Hp
p=0
subject to the constraints (4), wherein the min{·} symbol represents the minimization operation, NA represents the overall number of arcs relevant to the Delaunay triangulation, and w = (w0 , . . . , wp , . . . , wNA −1 ) sets the weights relevant to the weighted L1 norm. B. Temporal Unwrapping Let us start our discussion by reconsidering the spatial Delaunay triangulation in the AZ × RG plane (see, for example, the one in Fig. 2) and to assume that, for each given arc, δϕ = (δϕ0 , δϕ1 , . . . , δϕM −1 ) and δψ = (δψ0 , δψ1 , . . . , δψM −1 ) are the (measured) wrapped and (unknown) unwrapped DInSAR phase-gradient vectors of the DInSAR interferometric sequence, respectively. Moreover, in order to facilitate the procedure, we observe that a profitable model of the unknown unwrapped phase gradient vectors can be considered as follows [6], [7]: ∆b⊥ 4π · · ∆z + ∆t · ∆v (6) m(∆z, ∆v) ≈ λ r · sin(ϑ) wherein ∆z accounts for the error in the knowledge of the scene topography while ∆v accounts for the deformation velocity variations along the considered spatial arc. Moreover, r represents the sensor target distance, λ is the transmitted signal central wavelength, and ϑ is the incidence angle. Note that, because of the spatially correlated behavior of the atmospheric artifacts [6], we have assumed in (6) that they do not contribute to the model m(·). Since the wrapped phases differ from the unwrapped ones by 2π-integer multiples only, the DInSAR unwrapped vector can be expressed as follows: δψ = m + δϕ − m−π,π + 2πK = δχ + 2πK
(7)
wherein K is the 2π-integer multiple vector and δχ = m + δϕ − m−π,π
(8)
represents the phase component, for each considered arc, related to the computed model m(·) and to the measured phase contribution δϕ. Note also that, for the sake of simplicity, the explicit dependence of m(·) on the ∆z and ∆v factors has been neglected in (7) and (8); this notation will be maintained hereafter.
2377
At this stage, we can evaluate for each (∆z, ∆v) pair the unknown vector K in (7) by applying the basic MCF technique (summarized in the previous section) to the T × B⊥ grid. In particular, we search in this case for the solution of the following minimization problem: M −1 min C = |kj | (9) kj j=0
subject to the constraints
δχα + δχβ + δχγ kα + kβ + kγ = −round 2π
(10)
wherein kj is the jth element of the K vector in (7). Moreover, the round[·] symbol represents the operation of approximation to the nearest integer number, and, consistent with (3) and (4) in Section III-A, the constraints in (10) are expressed in terms of a generic triangle in the T × B⊥ plane whose arcs have been labeled as α, β, and γ, respectively. We remark that for each (∆z, ∆v) pair, we will have a different solution, which is characterized by its overall network cost C=
M −1
|kj |.
(11)
j=0
Accordingly, by exploring different values of the (∆z, ∆v) pair and the corresponding model m in (6), we may finally identify the unknown vector δψ opt = δχopt + 2πKopt
(12)
characterized by the “overall minimum” cost Cmin =
M −1
koptj
(13)
j=0
wherein δχopt and Kopt are the vectors (where koptj and δχoptj are the jth components, respectively) defined in (7) and (8) and evaluated in the minimum-cost condition. Accordingly, δψ opt in (12) will represent our estimate of the unwrapped vector for the considered arc. C. Spatial Unwrapping With regard to this operation, it exploits the unwrapped DInSAR phase differences computed via the previous temporal PhU step [see (12)] as a starting point for the spatialunwrapping operation. This second step is carried out on the single interferograms through the application of the MCF unwrapping technique [15] summarized in Section III-A. Moreover, in this case, we may use the minimum costs of the solutions achieved in the T × B⊥ plane, in order to set the weights wp [see (5)] associated to the single arcs. In particular, since the weights represent the confidence on the correctness of the achieved solution, we can reasonably assume that, for each spatial arc, the smaller is the temporal network cost, the larger
2378
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 9, SEPTEMBER 2006
TABLE I ERS 1/2 SAR DATA SET
Fig. 3. SAR products relevant to the investigated area and represented in the azimuth/range plane. (a) Multilook image. (b) Mask of the spatially coherent pixels (in black) defining the investigated spatial grid. (c) Corresponding spatial Delaunay triangulation. The latter is the same as that in Fig. 2 but it is repeated for the sake of completeness.
should be the weight applied within the minimization step [see (5)] involved in the spatial unwrapping. Accordingly, an inverse exponential relationship between the temporal network cost and the corresponding spatial weight wp has been chosen wp =
2S 2Cmin
1,
, Cmin < ρ Cmin ≥ ρ
(14)
wherein 2S represents the maximum allowed weight, Cmin is the minimum temporal cost defined in (13), with Cmin < S, and ρ is a threshold value that, based on experimental results, is typically set not greater than 5% of the total number of interferograms M . Following the implementation of this spatial-unwrapping step, which is carried out for each interferogram, the unwrapped phase can be easily retrieved via the integration of the computed spatial gradients. IV. R ESULTS In order to investigate the performance of the proposed unwrapping approach, a number of experiments, involving both real and simulated data, were carried out. The considered real data set is relevant to the Central Apennines (Abruzzi, Italy) area [see the multilook image shown in Fig. 3(a)], and is composed of 58 acquisitions collected by the ERS 1/2 sensors between August 1992 and October 2002 (see Table I), which are relevant to track 308 and frame 2755. The space/time data distribution is the one depicted in Fig. 1(a), while the final triangulation is shown in Fig. 1(c); the latter has been obtained by imposing a maximum spatial and temporal baseline value of 300 m and 1500 days, respectively, and a Doppler centroid separation, between the image pair, not larger than 1000 Hz. We further remark that, instead of imposing separate constraints on the maximum allowed spatial and temporal baselines, we could alternatively consider an upper limit for the normalized arc length defined in (1). Moreover, the value of the imposed constraints could be also changed depending on the characteristics of the investigated zone. For instance, larger values could be considered in the presence of highly coherent areas, as in the case of urbanized zones.
Following the triangulation step shown in Fig. 1(c), we have identified the sequence of DInSAR interferograms to be generated for implementing the proposed algorithm. In particular, we have computed 145 interferograms with four looks in the range direction and 20 looks in the azimuth one, having a pixel spacing of about 90 × 90 m. Note that the number of interferograms generated from each SAR image is dependent on the result of the Delaunay triangulation within the T × B⊥ plane (see Fig. 1). In this case, it is evident that acquisitions
PEPE AND LANARI: PHASE UNWRAPPING OF MULTITEMPORAL DIFFERENTIAL SAR INTERFEROGRAMS
2379
Fig. 4. Number of occurrences for each SAR image within the generated interferogram data set.
Fig. 6. Examples of the computed DInSAR interferograms and of the corresponding coherence maps. (a) and (d) have been computed from the ERS data pair acquired on July 15, 1997 and September 23, 1997, respectively (|∆b⊥ | = 39 m). (b) and (e) have been computed from the ERS data pair acquired on April 6, 1999 and September 28, 1999, respectively (|∆b⊥ | = 138 m). (c) and (f) have been computed from the ERS data pair acquired on October 28, 1997 and January 11, 2000, respectively (|∆b⊥ | = 207 m).
Fig. 5. Distribution of the temporal and spatial baselines relevant to the generated interferogram data set.
relevant to the triangulation boundaries are characterized by fewer arcs, i.e., are involved in a smaller number of interferograms; obviously, this effect becomes less significant when increasing the number of SAR acquisitions. With regard to our data set, we present in Fig. 4 the plot showing the number of appearances of each SAR image within the interferogram data set. Note that only two acquisitions contribute to two interferograms only. Moreover, we also present in Fig. 5 the distribution of the temporal and spatial baselines obtained by using our interferogram selection based on the Delaunay triangulation. We remark that, within the limits we defined as maximum baselines, the interferogram distribution is rather uniform, which is a relevant issue within the DInSAR scenarios. Based on the computed interferograms, we have generated the coherence mask shown in Fig. 3(b), i.e., the mask of the pixels that are considered coherent in our sequence of multilook interferograms. In our case, they have been selected by identifying the pixels that have a coherence, estimated in a box of four pixels in the range direction and 20 pixels in the azimuth one, greater than 0.35 in at least the 30% of the interferograms. Based on the mask shown in Fig. 3(b), the spatial triangulation has been generated [see Fig. 3(c)]. We remark that the selection of this Abruzzi test-site area is related to the difficulty in the unwrapping operation of the relevant interferograms caused
by the presence of strong decorrelation phenomena. To give an idea of these effects, we present in Fig. 6 a selection of computed interferograms. It is evident that the impact of the decorrelation effects can be very strong. Before discussing in detail the results achieved on the considered real SAR data, we prefer to start our analysis with some simulations that represent controlled experiments, allowing a direct assessment of the algorithm performance. To this end, we have retained both the AZ × RG and T × B⊥ triangulations shown in Figs. 1(c) and 3(c), respectively. Moreover, the same interferogram sequence generated for the real data set has been produced in the simulated experiments, but by exploiting the high-pass (HP) spatial components of a digital elevation model (DEM) of the area and by introducing a deformation parametric model. In particular, the jth interferogram of the simulated DInSAR sequence, say ϕsj , has the following expression: ϕsj
∆b⊥j 4π 1 2 · · h + ∆tj · vel + · ∆tj · acc = λ r · sin(ϑ) 2 j = 0, . . . , M − 1
(15)
wherein the first term represents the topographic phase contribution depending on the DEM height h, while the vel and acc factors are the velocity and acceleration model terms, respectively, the latter representing a significant deviation from the assumed model in (6). Moreover, in our experiments, we also assumed that both the simulated deformation components had a Gaussian spatial shape, as shown in Fig. 7, wherein both the masked and unmasked data are presented. In order to investigate the algorithm performance on the simulated data, we unwrapped the generated interferogram
2380
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 9, SEPTEMBER 2006
Fig. 7. Maps relevant to the three components of the model presented in (15) and used to generated the simulated DInSAR interferograms. (a) and (d) Unmasked and masked topography maps, respectively (hmin = −327 m, hmax = 428 m). (b) and (e) Unmasked and masked velocity patterns, respectively (velmin = −2.83 cm/year, velmax = −0.07 cm/year). (c) and (f) Unmasked and masked acceleration patterns, respectively (accmin = −14.15 cm/year2 , accmax = −0.06 cm/year2 ). Note that the position of the minima of the simulated velocity and acceleration patterns are relatively shifted by 100 pixels in both directions.
sequence both via the original and the extended MCF techniques; obviously, the former has been independently applied to each single interferogram. In the performed experiments, a comparison between the original and the unwrapped phases is directly possible. As a result of our analysis, we present in Fig. 8(a)–(c) the percentage, for each interferogram, of correctly unwrapped pixels of the mask shown in Fig. 3(b), obtained by applying the original MCF technique and considering on the horizontal axis the spatial and temporal baselines and the normalized arc length, respectively. Moreover, the results achieved via the extended MCF procedure are shown in Fig. 8(d), where the horizontal axis is relevant to the normalized arc length. The results presented in Fig. 8(a)–(c) clearly show that, for the original MCF approach, an increase of the spatial and/or of the temporal baseline (thus, of the normalized arc length) may lead to significant degradations of the achieved reconstruction. On the contrary, the extended MCF approach guarantees a very good retrieval, as shown in Fig. 8(d). With regard to the computational aspects, we remark that, on a personal computer equipped with an AMD Athlon MP2600+ processor, the basic MCF unwrapping procedure required about 75 min for unwrapping the overall interferogram sequence, while the extended MCF algorithm needed about 225 min. We remark that this computational increase is mostly related to the evaluation of the model function m within the temporalunwrapping step (see Section III-B). Indeed, this operation is carried out via an exhaustive search, with a limited spacing, of the ∆z and ∆v factors within the intervals (−100, 100) m and (−5, 5) cm/year, respectively. Obviously, the availability of an a priori information on the ∆z and ∆v terms would lead to
Fig. 8. Percentage of the coherent pixels of each interferogram that have been correctly unwrapped. (a) Results obtained by applying the original MCF approach with the horizontal axis relevant to the perpendicular baseline. (b) Same as (a) but with the horizontal axis relevant to the temporal baseline. (c) Same as (a) but with the horizontal axis relevant to the normalized arc length L defined in (1). (d) Results obtained by applying the extended MCF approach with the horizontal axis relevant to the normalized arc length L. (e) Results obtained by applying the extended MCF approach but without the model adjustment step; the horizontal axis is relevant to the normalized arc length L.
a significant computational improvement of the extended MCF PhU algorithm performance. As a final test on our simulated data, we investigated the quality of the retrieved solution if the model adjustment step, relevant to m in (6), is not implemented. In this case, the
PEPE AND LANARI: PHASE UNWRAPPING OF MULTITEMPORAL DIFFERENTIAL SAR INTERFEROGRAMS
2381
Fig. 9. Temporal coherence maps obtained from (a) the original and (b) the extended MCF unwrapping procedures, respectively.
required computing time reduced to about 80 min. However, the achieved results presented in Fig. 8(e) clearly show that the computational reduction is “paid” with a considerable degradation of the obtained unwrapped data with respect to those generated through the extended MCF approach, with the latter shown in Fig. 8(d). Let us now move to the real data analysis. In this case, the dimensions of the processed data set are the same as the simulated signals, thus the computing requirements remained practically unchanged; on the contrary, we now had to define a quality index for our results. To achieve this task, we decided to use the SBAS algorithm [7] that, as previously stated, is fully compatible with the presented data representation and allows us to generate DInSAR deformation time series from a sequence of unwrapped multitemporal interferograms. In particular, we considered the following strategy: first of all, we applied the SBAS procedure to both the unwrapped sequences computed via the original and the extended MCF algorithm; this allowed us to produce deformation time series for each pixel of the mask shown in Fig. 3(b). Subsequently, we used the computed time series (without any atmospheric filtering) and the estimated topographic errors to regenerate the original interferograms referred to as ϕ = (ϕ0 , ϕ1 , . . . , ϕM −1 ). At this stage, we use as quality index of the PhU retrieval the temporal coherence factor defined for each pixel as follows: M −1 exp j(ϕp − ϕp ) p=0 , 0 ≤ γ ≤ 1. (16) γ= M Note that for pixels where γ → 1, we expect that no unwrapping errors are present, since a nearly perfect retrieval of the original phase has been obtained. On the other hand, low values of γ will correspond to poorly unwrapped data. Following the application of the SBAS algorithm to the two unwrapped DInSAR sequences, obtained through the original and the extended MCF techniques, respectively, we have evaluated via (16) the temporal coherence map for both data sets. The achieved results are shown in Fig. 9(a) and (b), respectively. The improvement of the coherence obtained through the application
Fig. 10. Histograms of the temporal coherence maps. (a) The dashed line is related to the results shown in Fig. 9(a), which are obtained via the basic MCF approach. The continuous line is relevant to those achieved by using the extended MCF approach, which are presented in Fig. 9(b). (b) Same plots as in (a) with the superimposed diamond plot relevant to the results of the extended MCF approach implemented without 2π multiple corrections, i.e., by assuming Kopt = 0 in (12).
of the extended PhU algorithm is evident, particularly in the areas in the lower left and upper right corners of the map. To emphasize the achieved gain, the histograms relevant to these two results have been computed and superimposed in Fig. 10(a); again, the obtained improvement is clear. In particular, we remark that, in this case, we achieved an increase of about 60% with regard to the number of pixels with a temporal coherence value greater than the selected threshold γ = 0.7, which is a typical value in DInSAR applications [20], [21]. As a final point, we have investigated the impact, on the unwrapping algorithm performance, of the Kopt vector component in (12). To this end, we have repeated the overall processing operation based on the extended MCF approach but assuming no 2π integer correction, i.e., with Kopt = 0. This implies that, based on (12), we assumed δψ opt = δχopt . The histogram of the computed new temporal coherence has been superimposed [see Fig. 10(b)] to those shown in Fig. 10(a). The degradation of the achieved results is remarkable with respect to those of the overall extended MCF procedure. In particular, in this case, we had a reduction of the number of coherent pixels, with γ ≥ 0.7, of about 30% with respect to the extended MCF results. This clearly indicates the significant impact of the 2π integer correction component within the proposed unwrapping procedure.
2382
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 44, NO. 9, SEPTEMBER 2006
As a final remark, we underline that the overall analysis has been focused on multilook interferograms. However, the extension of this technique to the unwrapping of full-resolution DInSAR interferograms is under development. In this case, the noise phenomena affecting the single-look interferograms become extremely relevant; moreover, at the variance of the multilook case, we cannot benefit from the availability of an estimated coherence map generated via spatial averaging. Accordingly, an appropriate selection of temporally coherent pixels on the full-resolution interferograms is a crucial issue. We are presently investigating a solution based on the use of both multilook and single-look interferograms, which follows the lines of the full-resolution SBAS approach presented in [22]; this analysis can be a topic for future discussions. A PPENDIX
Fig. 11. Delaunay triangulations in the temporal/perpendicular baseline plane relevant to the ERS SAR data analyzed in the experiments of Section IV. The triangulations are obtained by assuming that the ratio between the perpendicular and temporal baseline axis units is equal to δT /δb⊥ , with (a) δT = 1500 days and δb⊥ = 300 m, (b) δT = 900 days and δb⊥ = 300 m, and (c) δT = 300 days and δb⊥ = 300 m.
V. C ONCLUSION We have proposed a solution for extending the MCF PhU algorithm dealing with a sparse data grid [15], to process multitemporal DInSAR interferograms for the generation of deformation time series. The approach involves the computation of a properly chosen interferogram sequence and is based on the cascade of two main steps, both involving the use of the basic MCF technique. In particular, we identify first the arcs connecting neighboring coherent pixels of the azimuth/range grid, and for each of these, we estimate the unwrapped phase gradients via the MCF technique applied in the temporal/perpendicular baseline plane. These results are used to bootstrap the subsequent spatial-unwrapping operation carried out on each single interferogram via the conventional MCF approach. The presented results obtained on simulated and real data confirm the effectiveness of the approach. Moreover, we underline that the proposed solution is naturally compatible with the SBAS algorithm for the generation of deformation time series. The implemented extension of the unwrapping procedure is based on the application of the network-programming algorithm for both the temporal- and spatial-unwrapping steps. Accordingly, the overall approach is computationally efficient.
Let us investigate, in the following, the impact of the different selection of the δT /δb⊥ factor, which represents the ratio between the temporal and perpendicular baseline axis units, which is needed in order to generate the triangulation in the T × B⊥ plane (see Section II). To achieve this task, we started from the triangulation relevant to our experiments presented in Section IV, where we have assumed δT = 1500 days and δb⊥ = 300 m. Subsequently, we maintained constant the value of δb⊥ and considered the two additional cases relevant to δT = 900 days and δT = 300 days, respectively. The results of the corresponding triangulations, presented in Fig. 11(a)–(c), respectively, show a strong similarity. More specifically, we found that more than 60% of the arcs (thus, of the corresponding interferograms) are common to the three triangulations. Accordingly, we may finally conclude that, within a realistic range of variation of the δT /δb⊥ ratio, the corresponding triangulations in the T × B⊥ plane do not significantly change. ACKNOWLEDGMENT Precise satellite orbits were provided by the University of Delft, The Netherlands. A Shuttle Radar Topography Mission DEM of the investigated zone has been used for the interferogram generation. The European Space Agency provided the ERS data relevant to the Abruzzi site through the CAT-1 proposal N. 313. The authors would like to thank the reviewers for their constructive remarks and suggestions. R EFERENCES [1] A. K. Gabriel, R. M. Goldstein, and H. A. Zebker, “Mapping small elevation changes over large areas: Differential interferometry,” J. Geophys. Res., vol. 94, no. B7, pp. 9183–9191, Jul. 1989. [2] D. Massonnet, M. Rossi, C. Carmona, F. Ardagna, G. Peltzer, K. Feigl, and T. Rabaute, “The displacement field of the Landers earthquake mapped by radar interferometry,” Nature, vol. 364, no. 6433, pp. 138– 142, Jul. 1993. [3] G. Peltzer and P. A. Rosen, “Surface displacement of the 17 Eureka valley, California, earthquake observed by SAR interferometry,” Science, vol. 268, no. 5215, pp. 1333–1336, Jun. 1995. [4] D. Massonnet, P. Briole, and A. Arnaud, “Deflation of Mount Etna monitored by spaceborne radar interferometry,” Nature, vol. 375, no. 6532, pp. 567–570, Jun. 1995. [5] E. Rignot, “Fast recession of a west Antarctic glacier,” Science, vol. 281, no. 5376, pp. 549–551, Jul. 1998.
PEPE AND LANARI: PHASE UNWRAPPING OF MULTITEMPORAL DIFFERENTIAL SAR INTERFEROGRAMS
[6] A. Ferretti, C. Prati, and F. Rocca, “Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry,” IEEE Trans. Geosci. Remote Sens., vol. 38, no. 5, pp. 2202–2212, Sep. 2000. [7] P. Berardino, G. Fornaro, R. Lanari, and E. Sansosti, “A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms,” IEEE Trans. Geosci. Remote Sens., vol. 40, no. 11, pp. 2375–2383, Nov. 2002. [8] O. Mora, J. J. Mallorquí, and A. Broquetas, “Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images,” IEEE Trans. Geosci. Remote Sens., vol. 41, no. 10, pp. 2243–2253, Oct. 2003. [9] S. Usai, “A least squares database approach for SAR interferometric data,” IEEE Trans. Geosci. Remote Sens., vol. 41, no. 4, pp. 753–760, Apr. 2003. [10] C. Werner, U. Wegmuller, T. Strozzi, and A. Wiesmann, “Interferometric point target analysis for deformation mapping,” in Proc. IGARSS, Toulouse, France, Jul. 2003, pp. 4362–4364. [11] A. Hooper, H. Zebker, P. Segall, and B. Kampes, “A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers,” Geophys. Res. Lett., vol. 31, no. 23, Dec. 2004, DOI: 10.1029/2004GL021737, CD-ROM. [12] M. Crosetto, B. Crippa, and E. Biescas, “Early detection and in-depth analysis of deformation phenomena by radar interferometry,” Eng. Geol., vol. 79, no. 1/2, pp. 81–91, Jun. 2005. [13] H. A. Zebker and J. Villasenor, “Decorrelation in interferometric radar echoes,” IEEE Trans. Geosci. Remote Sens., vol. 30, no. 5, pp. 950–959, Sep. 1992. [14] M. Costantini, “A novel phase unwrapping method based on network programming,” IEEE Trans. Geosci. Remote Sens., vol. 36, no. 3, pp. 813– 821, May 1998. [15] M. Costantini and P. A. Rosen, “A generalized phase unwrapping approach for sparse data,” in Proc. IGARSS, Hamburg, Germany, Jun. 1999, pp. 267–269. [16] M. Costantini, F. Malvarosa, F. Minati, L. Pietranera, and G. Milillo, “A three-dimensional phase unwrapping algorithm for processing of multitemporal SAR interferometric measurements,” in Proc. IGARSS, Toronto, ON, Canada, Jun. 2002, pp. 1741–1743. [17] A. Pepe and R. Lanari, “A space-time minimum cost flow phase unwrapping algorithm for the generation of DInSAR deformation time-series,” in Proc. IGARSS, Seoul, Korea, Jul. 2005, pp. 1979–1982. [18] G. Franceschetti and R. Lanari, Synthetic Aperture Radar Processing. Boca Raton, FL: CRC, Mar. 1999. [19] N. Miranda, B. Rosich, C. Santella, and M. Grion, “Review of the impact of ERS-2 piloting modes on the SAR Doppler stability,” in Proc. Fringe, Frascati, Italy, Dec. 2003, CD-ROM. [20] R. Lanari, P. Lundgren, M. Manzo, and F. Casu, “Satellite radar interferometry time series analysis of surface deformation for Los Angeles, California,” Geophys. Res. Lett., vol. 31, no. 23, Dec. 2004, DOI: 10.1029/2004GL021294, CD-ROM. [21] A. Borgia, P. Tizzani, G. Solaro, M. Manzo, F. Casu, G. Luongo, A. Pepe, P. Berardino, G. Fornaro, E. Sansosti, G. P. Ricciardi, N. Fusi, G. Di Donna, and R. Lanari, “Volcanic spreading of Vesuvius, a new paradigm for interpreting its volcanic activity,” Geophys. Res. Lett., vol. 32, no. 3, Feb. 2005, DOI: 10.1029/2004GL022155, CD-ROM. [22] R. Lanari, O. Mora, M. Manunta, J. J. Mallorquí, P. Berardino, and E. Sansosti, “A small baseline approach for investigating deformations on full resolution differential SAR interferograms,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 7, pp. 1377–1386, Jul. 2004.
2383
Antonio Pepe was born in Salerno, Italy, in 1974. He received the laurea degree in electronic engineering from the University of Napoli “Federico II,” Naples, Italy, in 2000, with a thesis on the synthetic aperture radar (SAR) image processing. He is currently working toward the Ph.D. degree at the same institution. He joined the Istituto per il Rilevamento Elettromagnetico dell’Ambiente (IREA), Institute of the Italian National Research Council (CNR), Naples, in 2002. Since 2003, he has been with the Department of Electronic and Telecommunication Engineering, University of Napoli. In 2005, he was a Visiting Research at the Aerospace Engineering and Engineering Mechanics, The University of Texas at Austin. His research interests concern the differential SAR interferometry applications for the monitoring of surface displacements, such as those produced by subsidence, volcano activity, and earthquakes.
Riccardo Lanari (M’91–SM’01) was born in 1964. He received the degree in electronic engineering (summa cum laude) from the University of Napoli, “Federico II,” Naples, Italy, in 1989. After graduation, he joined the Istituto di Ricerca per l’Elettromagnetismo e i Componenti Elletronici (IRECE), now the Istituto per il Rilevamento Elettromagnetico dell’Ambiente (IREA), a Research Institute of the Italian National Research Council (CNR), where he currently occupies the position of Senior Researcher. He has been a Visiting Scientist at different foreign research institutes. He was a Research Fellow at the Institute of Space and Astronautical Science, Japan in 1993; a Visiting Scientist at the German Aerospace Research Establishment (DLR), in 1991 and 1994, and at the Jet Propulsion Laboratory, in 1997. He has lectured in several national and foreign universities and research centers and is currently a Lecturer of the SAR module course of the International Master in Airborne Photogrammetry and Remote Sensing, offered by the Institute of Geomatics in Barcelona, Spain. He is the holder of two patents on SAR data-processing techniques. His main research interests are in the synthetic aperture radar data-processing field as well as in SAR interferometry techniques. On these fields, he has authored 40 international journal papers, and, in 1999, he authored a book entitled Synthetic Aperture Radar Processing (CRC Press, 1999). Mr. Lanari received a NASA recognition for the development of the ScanSAR processor used for the SRTM mission. He has been invited as Chairman and Cochairman at several international conferences and as a member of the Technical Program Committee for the IGARSS’2000 symposia.
