- Email: [email protected]
Generation of three-dimensional deformation maps from InSAR data using spectral diversity techniques
ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
Contents lists available at ScienceDirect
ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs
Generation of three-dimensional deformation maps from InSAR data using spectral diversity techniques E. Erten a,b,∗ , A. Reigber b , O. Hellwich a a
Computer Vision and Remote Sensing, Technical University of Berlin, 10587 Berlin, Germany
b
Microwaves and Radar Institute, German Aerospace Center, 82234 Oberpfaffenhofen, Germany
article
info
Article history: Received 5 June 2007 Received in revised form 21 April 2010 Accepted 22 April 2010 Available online 1 June 2010 Keywords: Spectral diversity DInSAR 3D deformation analysis Weighted least squares
abstract The capability of DInSAR (Differential Interferometric SAR) for precise large-scale deformation analysis has been shown in various case studies. Generally, DInSAR possesses a high potential for monitoring deformation, but only the velocity component parallel to the line-of-sight direction can be measured. An alternative approach, capable to retrieve the deformation velocity in both range and azimuth direction, is the so-called spectral diversity technique. Spectral diversity is based on a phase comparison between different sub-aperture interferograms of the scene and can generally be regarded as a high-performance technique for estimating the mis-registration between complex SAR images. In this paper, the following questions will be discussed: how to implement the spectral diversity technique for achieving the most accurate results; how to extract the full 3D deformation vector from a combination of ascending/descending passes and how to extract a surface deformation map if the data sets are not perfectly coherent. Finally, a statistical analysis of every individual processing step and an error propagation analysis is undertaken. In order to make a quantitative analysis of the technique, ENVISAT data sets of the Bam earthquake in Iran are used. © 2010 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.
1. Introduction In recent years, a special interest in operational deformation monitoring has emerged. Up until now, two main approaches have been used for the estimation of ground deformation with remotely sensed data: image matching for all types of remotely sensed data (Erten et al., 2009; Evans, 2000; Leberl and Maurice, 1994; Sarti et al., 2006; Strozzi et al., 2002), and interferometric techniques using SAR data (Hanssen, 2001; Bamler and Hartl, 1998). Intensity-based matching approaches define the displacement as the shift that yields the best fit between different images in time and they are not limited by phase stability problems and can be reliably acquired on a regular basis but the object-dependent phase information is not taken into account. In contrast, DInSAR can detect ground deformation with up to millimeter accuracy by considering the phase information contained in complex SAR images. However, one main limitation of deformation measurements made with interferometry is that an interferogram delivers only
∗ Corresponding author at: Computer Vision and Remote Sensing, Technical University of Berlin, 10587 Berlin, Germany. Tel.: +49 8153 28 3350; fax: +49 8153 28 1449. E-mail address: [email protected] (E. Erten).
the component of the deformation in line-of-sight (L.O.S.) direction of the sensor. In other words, DInSAR is not sensitive to detect the components of the deformation vector which are perpendicular to the L.O.S. direction, i.e. occurring in along-track (azimuth) direction. Nowadays, there is a growing interest in mapping surface deformation phenomena in three dimensions (3D). Regarding this demand, DInSAR and the matching algorithms have been combined to retrieve deformation information in the L.O.S. and in the azimuth direction. For example, normalised incoherent crosscorrelation (NICC) has been used to find horizontal displacement in the Hector Mine earthquake by Fialko and Simons (2001) and to detect glacier movement in Alps with ERS data by Strozzi et al. (2002). Sarti et al. (2006) performed the co-seismic fault rupture detection using coherence maximisation in the Bam earthquake area. Although ascending and descending pairs were used to obtain 2D and 3D deformation vectors in all these papers, the direct (without appropriate weighting) combination of different techniques, i.e., DInSAR and incoherent cross-correlation, considering their different accuracies was not performed. In this paper another interferometric coregistration technique called Spectral Diversity (SD), or split-band interferometry, is used to extract a 3D deformation map. Originally, this technique was developed for the coregistration of SAR interferograms and for the detection of residual motion errors in airborne SAR data (Scheiber and Moreira, 2000; Prats et al., 2004). Since spectral diversity allows to precisely coregister SAR image pairs, this technique
0924-2716/$ – see front matter © 2010 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved. doi:10.1016/j.isprsjprs.2010.04.005
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
may also be applied to precisely estimate the local image offsets caused by surface movement in range and azimuth directions. Although the spectral diversity technique is well suited for deformation monitoring, this technique also requires a coherent SAR data set as in the case of DInSAR applications. It is less accurate than DInSAR, since the sensibility of deformation detection is linked to the resolution of the sensor and not to its wavelength as in DInSAR. However, it is known that spectral diversity can acquire a significantly higher accuracy than cross-correlation approaches (Scheiber and Moreira, 2000) and is computationally very efficient. In this paper, two-dimensional shift measurements based on spectral diversity technique of ascending and descending satellite passes were combined with the aim of extracting of a 3D deformation map. A weighted least squares approach was used for homogenising the measurements. In order to obtain a general view of surface displacement, Delaunay triangulation and smooth interpolation were performed. Functional and stochastic relations allow the statistical analysis of every individual processing step and an error propagation analysis regarding the final results. Finally, it is shown that two interferometric pairs are sufficient to infer a complete 3D surface displacement by the use of amplitude and phase information of interferometric SAR data. 2. Methodology to retrieve 3D deformation vector with spectral diversity This section presents the methodology to retrieve the 3D deformation extraction using spectral diversity technique. In particular, firstly the spectral diversity technique and the image acquisition geometry will be briefly introduced. Then, a model-based inversion, i.e., weighted least squares approach, is carried out to retrieve the deformation vectors using Delaunay triangulation. 2.1. Spectral diversity The spectral diversity technique is applied to an interferometric pair by calculating the phase differences between the interferograms generated from two bandpass-filtered images (Bamler and Eineder, 2005; Madsen and Zebker, 1992; Massonnet et al., 1996; Scheiber and Moreira, 2000). This corresponds to precisely measuring the mis-registration between image amplitudes. Shortly, SD includes the generation of two low resolution interferogram Ď Ď I1 = M1 S1 and I2 = M2 S2 having the interferometric phase obtained from the first sub-look, i.e., M1 and S1 , and the second sublook, i.e., M2 and S2 , of the master and the slave images indicated by M and S, respectively. By combining the first (I1 ) and the second look (I2 ) interferograms in the following manner
φsd = M1 S1Ď (M2 S2Ď )Ď
(1)
a differential interferogram is retrieved. Here, the phase information is related to the coregistration error as a fractional number of resolution as
∆t =
arg{φsd } 2π (fc1 − fc2 )
(2)
where fc1 and fc2 are the central frequencies of the two sub-looks. One also sees that in contrast to interferometric approaches, this technique can be performed both in range and azimuth, delivering mis-registration maps in both directions. For this purpose, sub-aperture images and the corresponding interferograms are computed both in range and azimuth direction. It has to be mentioned that the coregistration accuracy of SD is significantly higher than the one of other techniques (Bamler and Eineder, 2005) and can reach precisions of up to chem1/100 pixels (Scheiber and Moreira, 2000). Further details about the spectral diversity approach can be found in Scheiber and Moreira (2000). For the extraction of deformation vectors, spectral diversity has to be applied after compensation of mis-registration caused by
389
Fig. 1. The relationship between interferometric phase observations in L.O.S. dr and in azimuth da direction and ground deformation considering the right looking ascending satellite geometry.
the scene’s topography. This can be accomplished by geometrical transforms using precise orbit information and a Digital Elevation Model (DEM) as explained by Fornaro et al. (2005) and Reigber et al. (2006). 2.2. Retrieval of the 3D deformation vector Up to now, it has been explained how two-dimensional (2D) deformation vectors can be extracted by means of spectral diversity. It remains the question how to construct the complete 3D deformation vectors. For this aim, the relationship between offset maps produced in range and azimuth direction and the ground deformation has to be defined. The projection of a 3D displacement
− →
vector d = [dv dn de ]T onto line-of-sight direction and satellite heading direction can be expressed as (compare Fig. 1),
− →A d
r
− →A d
a
= δr + dv cos(θinc ) − sin(θinc )[dn sin(H ) + de cos(H )] = δa + dn cos(H ) − de sin(H )
(3)
where θinc and H denote the incidence angle and the heading angle
− →
(positive clockwise from the North), respectively, d = [dv dn de ]T indicates the deformation vector in vertical, north and east direction, and the superscripts of A shows ascending geometry. Since there are three unknowns for the displacement, another interferometric observation from a descending orbit is required: data from ascending and descending orbits have different L.O.S. and azimuth directions, providing two linearly independent measurements. With D indicating descending geometry, one obtains
− →D d
r
− →D d
a
= δr + dv cos(θinc ) − sin(θinc )[dn sin(H ) − de cos(H )] = δa + dn cos(H ) − de sin(H ).
(4)
Although these equations are directly solvable with simple linear equation techniques, this is not necessarily optimal when the observations are obtained with different geometry and have not the same phase accuracies. Instead, a Weighted Least Squares (WLS) approach can be used. Weighting is performed by the phase standard deviation, which is highly related to the interferometric coherence. In this way, observations with high accuracy (high coherence value = low standard deviation) get higher weights and will, therefore, have a stronger influence on the estimated
390
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
Table 1 Co-seismic interferometric pairs. Co-seismic Pairs
Baseline (m)
Descending
Ascending
03.12.2003–11.02.2004 03.12.2003–07.01.2004 03.12.2003–07.01.2004 03.12.2003–11.02.2004
16.11.2003–29.02.2004 16.11.2003–29.02.2004 16.11.2003–29.02.2004 16.11.2003–29.02.2004
Temporal resolution
0.6–2 570–2 570–30 0.6–30
70–105 35–105 35–70 70–70
Fig. 2. Block diagram of the implemented processing chain to retrieve the 3D deformation map. 58°E
unknowns and vice versa. The weighting matrix P, which is inversely proportional to the covariance matrix of observations, is defined as in the following: 29°30'N
(5) 29°15'N
1/σ
2 Da
.
σ = √
B
2 N B−b
1 − γ2
1 28°45'N
1
29°N
Here the subscripts of D and A refer to the descending and ascending geometries, respectively, while r and a refer to the observation at the range and the azimuth directions, respectively. I.e., σAr indicates the standard deviation of the observation at the ascending geometry in range direction. The phase standard deviations defined by Bamler and Eineder (2005) can be retrieved using
29°N
1/σ 0
2 Dr
0 0 0
59°E
03.12.2003 - 11.02.2004 0 30 km
29°15'N
2 1/σAa 0 0
0 0
58°45'E
28°45'N
0
58°30'E
29°30'N
2 1/σAr 0 P = 0 0
58°15'E
r p B b
πγ
0
(6)
where N is the number of independent samples used in estimation of the SD phase, γ is the interferometric coherence, B is the bandwidth of the full-resolution signal and b is the new bandwidth after bandpass filtering. It may noted that b = 3B has been reported as a optimum bandwidth ratio by Bamler and Eineder (2005) for the estimation of differential shifts. − → An expression for d , which minimises the estimation error, can be derived as in the following
58°E
58°15'E
58°30'E
58°45'E
59°E
Fig. 3. The coherence image retrieved from the dataset having a 2 m spatial baseline and 70 days temporal resolution.
where Q is the matrix including the coefficients in Eqs. (3) and (4) regarding the direction of the observations, ∆ is the vector including observations from Eq. (2) and P is the weighting matrix formed as in Eq. (5).
was employed using SRTM-3 DEM. In the backward geocoding approach, for each DEM element, the image pixel with the nearest range-Doppler coordinate is calculated (Reigber et al., 2006). As reported by Sansosti et al. (2006), the accuracy of 1/8th pixel can be achieved with backward geocoding using SRTM-3 in the case of ENVISAT data. Fig. 2 shows the processing chain to retrieve the deformation map from the SD technique. As it can be observed in Fig. 3, the effect of the earthquake is already visible in the coherence map, where the city center (approx. in the middle) appears decorrelated due to the severe building damages caused by the earthquake; while normally urban areas are known to show a very high correlation. Moreover, the earthquake’s fault line can also be detected in the coherence map.
3. Experimental results
3.1. Deformation map over the city of Bam
The proposed SD technique with the aim of detection of 3D deformation was tested on co-seismic interferometric ENVISAT ASAR images, which are summarised in Table 1, over the city of Bam in the south-eastern part of Iran, which suffered from a strong earthquake (M v = 6.5) in winter 2003/2004. It is worth again mentioning that the precise geometric coregistration of SAR data is a strict requirement for deformation monitoring. In order to perform a precise coregistration, backward geocoding
This section presents and discuss the detection of ground deformation based on DInSAR and spectral diversity methods from co-seismic interferometric data sets. For this purpose, 3 descending orbits and 3 ascending orbits have been processed. In Fig. 4(a), a descending co-seismic interferogram, depicting the slant-range displacement due to the earthquake is shown, generated from the 03 December 2003 (before earthquake) and 07 January 2004 (after earthquake) ASAR products. Because of
− →
d =
Q T P∆ Q T PQ
(7)
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
a
b
-180°
(degree)
180°
391
c
1.20
(m)
-1.30
-2.50
(m)
-2.90
Fig. 4. (a) Line-of-sight displacement map with differential interferogram. (b) Along-track displacement (azimuth direction) map with spectral diversity. (c) Line-of-sight displacement map with spectral diversity. Sub-look images have been obtained by selecting 1/3 of the full bandwidth.
Fig. 5. Profiles from azimuth offset map in the left–right direction in meters.
the large displacement, the deformation pattern appears wrapped and its gradient could not be measured directly by DInSAR due to phase unwrapping and noise problems in the area close to the fault. In this area, the extrapolation of fringes or similar techniques could be considered. From the same data pair, spectral diversity deformation maps in azimuth and range direction, respectively, were derived, shown in Fig. 4(b–c). In contrast to Fig. 4(a), the fault line, which is located in the north–south direction, can now easily be detected without any a priori knowledge. In order to derive a quantitative analysis, the deformations along six profiles, marked as black lines in Fig. 4(b–c), were analysed. As a result of the earthquake, surface movements of various amounts could be found (Figs. 5–7). A discontinuity and a sign change at the fault crossing can be found in azimuth profiles, see Fig. 5. An analysis of the azimuth displacements also shows that the amount of deformation decreases from north to south
direction. In the first horizontal profile, which is close to the city center, the estimated amount of azimuth deformation changes by about 70 cm from the beginning of the profile to the end of the profile. In the last azimuth profile, this change is only between −0.01 and −0.22 cm (21 cm along the profile line). However, the SNR (signal to noise ratio) in range is much lower than the one in azimuth. The discontinuity and the increase of deformation along the range profiles can been observed in Fig. 6. In addition to these ‘‘left-to-right’’ profile analysis, a ‘‘south–north’’ profile analysis was performed on the azimuth deformation map, shown in Fig. 7. This figure supports the assumption that both sides of the fault move away from each other: one of them tends more towards the north direction, the other more towards the south direction. This fact is for example not visible in a descending differential interferogram, because the surface motion happens almost perfectly perpendicular to the radar line-of-sight, where the radar has no sensitivity.
392
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
Fig. 6. Profiles from range offset map in the left–right direction in meters. Note that the displacement map in range direction is close to noise level, its weight is lower in 3D vector extraction than that in azimuth direction.
Fig. 7. Profiles from azimuth offset map in the south–north direction in meters.
After applying the processing chain expounded in Fig. 2 to each group of data sets (see Table 1), several SD deformation maps (after spatial averaging) were combined using weighted least squares approach to deduce the full 3D deformation field
− →
of the Bam earthquake. The unknown parameters, i.e., d = [dv dn de ]T , were not solved for each pixel but rather for pixels with sufficient interferometric coherence, i.e., γ = 0.7. With this criterion, only the pixels containing reliable phase information were used in the calculations. Delaunay triangulation was used for interpolation and visual display of irregularly gridded data considering the selected pixels. After triangulation, a (smoothed) view of 3D deformation and surface motion in all directions can been obtained for each pixel. Fig. 8(a) shows the vertical deformation, derived using the proposed approach, which fits very well with DInSAR. In this figure it can be seen how the vertical deformation changes its sign through the fault line. This is consistent with the DInSAR since, from west to east, the colour cycle has a negative phase value in the southern part of the earthquake, which shows that the surface
motion is towards the sensor. The case of northern part, colour cycle shows that the surface deformation is away from the sensor (see Fig. 4(a)). Fig. q 8(b) shows the magnitude of the horizontal deformation, i.e.,
− → 2
− →
dn + d2e .
3.2. Validation The accuracy of the estimated deformation is directly related to the estimated coherence of the interferometric pair, which indicate the phase quality. According to this, the accuracy assessment of the proposed technique was performed based on the coherence level of data set. Fig. 9 shows the standard deviation of the splitband interferometric phase in azimuth and in range direction and the interferometric phase as a function of coherence, for three multilook levels regarding Eq. (6) and Hanssen (2001, eq. 4.2.27) respectively. As expected, the accuracy of the interferometric measurement (without consideration of phase unwrapping) is more accurate than SD ones for an equal number of samples with
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
a
393
b 10
50
(cm)
(cm)
-10
0
Fig. 8. (a) Vertical and (b) horizontal components of surface displacement obtained from the data set in Table 1 having 0.6 m baseline with descending acquisition and 30 m baseline with ascending acquisition.
a
b
c
Along-Track offsets (meters)
a
0.6
b
0.6
0.4
Slant-Range offsets (meters)
Fig. 9. Standard deviations of the split-band interferometric phase in azimuth (a) and in range (b) direction and the interferometric phase (c) as a function of coherence, for three multilook levels.
0.4
0.2
0.0
data set with 0.6 m baseline data set with 570 m naseline
-0.2
0.0
data set with 570 m naseline
-0.2
data set with 0.6 m baseline
-0.4
-0.6
0.2
-0.4
0
50
100
150
200
250
samples
-0.6 0
50
100
150
200
250
samples
Fig. 10. The averaged profiles of displacement and their error bars in azimuth (a) and range (b).
fixed correlation. Moreover, it can be seen that the accuracy in range direction by SD is much lower than in azimuth direction by SD. With the threshold of the coherence of γ = 0.7 using 49 × 49 samples, the worst accuracy of 9 cm in azimuth and 16 cm in range direction have been predicted. In order to validate our results regarding the different combination of data sets, we plotted azimuth and range offsets and their error bars along the profiles in the area of maximum deformation. As expected and can be clearly seen from Fig. 10 that offsets and their variances obtained from data sets having various spatial and temporal baselines are different. Due to this the final
product, i.e., Fig. 8, has been retrieved from the ascending and the descending data set that have a maximum coherence. 4. Conclusions and discussions A comparison of spectral diversity and interferometry techniques for detecting of ground deformation has been carried out. Two-dimensional offset measurements based on spectral diversity technique from two complementary ascending and descending passes were combined using weighted least squares approach with the aim of extracting a fully 3D deformation map. In order to obtain
394
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
a general view of surface deformation map, Delaunay triangulation and smooth interpolation were performed. It is demonstrated that the spectral diversity technique can not only be used for precise coregistration, but also for resolving 3D movement vectors in vertical, north–south and east–west directions. It is also demonstrated that 3D deformation map can be extracted from only two interferometric pairs, which is not possible by means of DInSAR. However, with spectral diversity technique, only decimeter accuracy is possible in case of ENVISAT, in contrast to millimeter accuracy in case of DInSAR. However, at least for analysing large co-seismic motion, this seems to be sufficient. Additionally, new sensors, like for example TerraSAR-X and RADARSAT-2, which provide a much higher spatial resolution, will significantly enhance the precision of the proposed technique. Acknowledgements This work has been supported by the German Academic Exchange Service (DAAD). The authors would also like to thank the European Space Agency (ESA) for providing the ENVISAT data sets of the Bam region. References Bamler, R., Hartl, P., 1998. Synthetic aperture radar interferometry. Inverse Problems 14 (4), 1–54. Bamler, R., Eineder, M., 2005. Accuracy of differential shift estimation by correlation and split-bandwidth interferometry for wideband and delta-k SAR systems. IEEE Geoscience and Remote Sensing Letters 2 (2), 151–155. Erten, E., Reigber, A., Hellwich, O., Prats, P., 2009. Glacier velocity monitoring by maximum likelihood texture tracking. IEEE Transactions on Geoscience and Remote Sensing 47 (2), 394–405.
Evans, A.N., 2000. Glacier surface motion computation from digital image sequences. IEEE Transactions on Geoscience and Remote Sensing 38 (2), 1064–1072. Fialko, Y., Simons, M., 2001. The complete (3-D) surface displacement field in the epicentre area of the 1999 Mw 7.1 Hector Mine earthquake, California, from space geodetic observations. Geophysical Research Letters 28 (16), 3063–3066. Fornaro, G., Manunta, M., Serafino, F., Berardino, P., Sansosti, E., 2005. Advances in multipass SAR image registration. In: Proceedings of IGARSS’05, Seoul, SouthKorea, 25–29 July, pp. 4832–4835. Hanssen, F.R., 2001. Radar Interferometry: Data Interpretation and Error Analysis. Kluwer Academic Publishers. Leberl, F.W., Maurice, K., 1994. Automated radar image matching experiment. ISPRS Journal of Photogrammetry and Remote Sensing 49 (3), 19–33. Madsen, S.N., Zebker, H.A., 1992. Automated absolute phase retrieval in acrosstrack interferometry. In: Proceedings of IGARSS’92, Houston, TX, 26–29 May, 2, pp. 1582–1584. Massonnet, D., Vadon, H., Rossi, M., 1996. Reduction of the need for phase unwrapping in radar interferometry. IEEE Transactions on Geoscience and Remote Sensing 34 (2), 489–497. Prats, P., Reigber, A., Mallorqui, J.J., Broquetas, A., 2004. Interpolation-free coregistration and phase-correction of airborne SAR interferograms. IEEE Geoscience and Remote Sensing Letters 2 (2), 206–210. Reigber, A., Erten, E., Guillaso, S., Hellwich, O., 2006. I.D.I.O.T.: a free and easyto-use software tool for DInSAR analysis. 26. Wissenschaftlich-Technische Jahrestagung der DGPF. (on CDROM), pp. 287–293. Sansosti, E., Berardino, P., Manunta, M., Serafino, F., Fornaro, G., Irea, N., 2006. Geometrical SAR image registration. IEEE Transactions on Geoscience and Remote Sensing 44 (10-2), 2861–2870. Sarti, F., Briole, P., Pirri, M., 2006. Co-seismic fault rupture detection and slip measurement by ASAR precise correlation using coherence maximization: application to a north-south blind fault in the vicinity of Bam (Iran). IEEE Geoscience and Remote Sensing Letters 3 (2), 1–5. Scheiber, R., Moreira, A., 2000. Coregistration of interferometric SAR images using spectral diversity. IEEE Transactions on Geoscience and Remote Sensing 38 (5), 2179–2191. Strozzi, T., Kuckman, A., Nurray, T., Wegmüller, U., Werner, C.L., 2002. Glacier motion estimation using SAR offset-tracking procedures. IEEE Transactions on Geoscience and Remote Sensing 40 (11), 2384–2390.
Contents lists available at ScienceDirect
ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs
Generation of three-dimensional deformation maps from InSAR data using spectral diversity techniques E. Erten a,b,∗ , A. Reigber b , O. Hellwich a a
Computer Vision and Remote Sensing, Technical University of Berlin, 10587 Berlin, Germany
b
Microwaves and Radar Institute, German Aerospace Center, 82234 Oberpfaffenhofen, Germany
article
info
Article history: Received 5 June 2007 Received in revised form 21 April 2010 Accepted 22 April 2010 Available online 1 June 2010 Keywords: Spectral diversity DInSAR 3D deformation analysis Weighted least squares
abstract The capability of DInSAR (Differential Interferometric SAR) for precise large-scale deformation analysis has been shown in various case studies. Generally, DInSAR possesses a high potential for monitoring deformation, but only the velocity component parallel to the line-of-sight direction can be measured. An alternative approach, capable to retrieve the deformation velocity in both range and azimuth direction, is the so-called spectral diversity technique. Spectral diversity is based on a phase comparison between different sub-aperture interferograms of the scene and can generally be regarded as a high-performance technique for estimating the mis-registration between complex SAR images. In this paper, the following questions will be discussed: how to implement the spectral diversity technique for achieving the most accurate results; how to extract the full 3D deformation vector from a combination of ascending/descending passes and how to extract a surface deformation map if the data sets are not perfectly coherent. Finally, a statistical analysis of every individual processing step and an error propagation analysis is undertaken. In order to make a quantitative analysis of the technique, ENVISAT data sets of the Bam earthquake in Iran are used. © 2010 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.
1. Introduction In recent years, a special interest in operational deformation monitoring has emerged. Up until now, two main approaches have been used for the estimation of ground deformation with remotely sensed data: image matching for all types of remotely sensed data (Erten et al., 2009; Evans, 2000; Leberl and Maurice, 1994; Sarti et al., 2006; Strozzi et al., 2002), and interferometric techniques using SAR data (Hanssen, 2001; Bamler and Hartl, 1998). Intensity-based matching approaches define the displacement as the shift that yields the best fit between different images in time and they are not limited by phase stability problems and can be reliably acquired on a regular basis but the object-dependent phase information is not taken into account. In contrast, DInSAR can detect ground deformation with up to millimeter accuracy by considering the phase information contained in complex SAR images. However, one main limitation of deformation measurements made with interferometry is that an interferogram delivers only
∗ Corresponding author at: Computer Vision and Remote Sensing, Technical University of Berlin, 10587 Berlin, Germany. Tel.: +49 8153 28 3350; fax: +49 8153 28 1449. E-mail address: [email protected] (E. Erten).
the component of the deformation in line-of-sight (L.O.S.) direction of the sensor. In other words, DInSAR is not sensitive to detect the components of the deformation vector which are perpendicular to the L.O.S. direction, i.e. occurring in along-track (azimuth) direction. Nowadays, there is a growing interest in mapping surface deformation phenomena in three dimensions (3D). Regarding this demand, DInSAR and the matching algorithms have been combined to retrieve deformation information in the L.O.S. and in the azimuth direction. For example, normalised incoherent crosscorrelation (NICC) has been used to find horizontal displacement in the Hector Mine earthquake by Fialko and Simons (2001) and to detect glacier movement in Alps with ERS data by Strozzi et al. (2002). Sarti et al. (2006) performed the co-seismic fault rupture detection using coherence maximisation in the Bam earthquake area. Although ascending and descending pairs were used to obtain 2D and 3D deformation vectors in all these papers, the direct (without appropriate weighting) combination of different techniques, i.e., DInSAR and incoherent cross-correlation, considering their different accuracies was not performed. In this paper another interferometric coregistration technique called Spectral Diversity (SD), or split-band interferometry, is used to extract a 3D deformation map. Originally, this technique was developed for the coregistration of SAR interferograms and for the detection of residual motion errors in airborne SAR data (Scheiber and Moreira, 2000; Prats et al., 2004). Since spectral diversity allows to precisely coregister SAR image pairs, this technique
0924-2716/$ – see front matter © 2010 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved. doi:10.1016/j.isprsjprs.2010.04.005
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
may also be applied to precisely estimate the local image offsets caused by surface movement in range and azimuth directions. Although the spectral diversity technique is well suited for deformation monitoring, this technique also requires a coherent SAR data set as in the case of DInSAR applications. It is less accurate than DInSAR, since the sensibility of deformation detection is linked to the resolution of the sensor and not to its wavelength as in DInSAR. However, it is known that spectral diversity can acquire a significantly higher accuracy than cross-correlation approaches (Scheiber and Moreira, 2000) and is computationally very efficient. In this paper, two-dimensional shift measurements based on spectral diversity technique of ascending and descending satellite passes were combined with the aim of extracting of a 3D deformation map. A weighted least squares approach was used for homogenising the measurements. In order to obtain a general view of surface displacement, Delaunay triangulation and smooth interpolation were performed. Functional and stochastic relations allow the statistical analysis of every individual processing step and an error propagation analysis regarding the final results. Finally, it is shown that two interferometric pairs are sufficient to infer a complete 3D surface displacement by the use of amplitude and phase information of interferometric SAR data. 2. Methodology to retrieve 3D deformation vector with spectral diversity This section presents the methodology to retrieve the 3D deformation extraction using spectral diversity technique. In particular, firstly the spectral diversity technique and the image acquisition geometry will be briefly introduced. Then, a model-based inversion, i.e., weighted least squares approach, is carried out to retrieve the deformation vectors using Delaunay triangulation. 2.1. Spectral diversity The spectral diversity technique is applied to an interferometric pair by calculating the phase differences between the interferograms generated from two bandpass-filtered images (Bamler and Eineder, 2005; Madsen and Zebker, 1992; Massonnet et al., 1996; Scheiber and Moreira, 2000). This corresponds to precisely measuring the mis-registration between image amplitudes. Shortly, SD includes the generation of two low resolution interferogram Ď Ď I1 = M1 S1 and I2 = M2 S2 having the interferometric phase obtained from the first sub-look, i.e., M1 and S1 , and the second sublook, i.e., M2 and S2 , of the master and the slave images indicated by M and S, respectively. By combining the first (I1 ) and the second look (I2 ) interferograms in the following manner
φsd = M1 S1Ď (M2 S2Ď )Ď
(1)
a differential interferogram is retrieved. Here, the phase information is related to the coregistration error as a fractional number of resolution as
∆t =
arg{φsd } 2π (fc1 − fc2 )
(2)
where fc1 and fc2 are the central frequencies of the two sub-looks. One also sees that in contrast to interferometric approaches, this technique can be performed both in range and azimuth, delivering mis-registration maps in both directions. For this purpose, sub-aperture images and the corresponding interferograms are computed both in range and azimuth direction. It has to be mentioned that the coregistration accuracy of SD is significantly higher than the one of other techniques (Bamler and Eineder, 2005) and can reach precisions of up to chem1/100 pixels (Scheiber and Moreira, 2000). Further details about the spectral diversity approach can be found in Scheiber and Moreira (2000). For the extraction of deformation vectors, spectral diversity has to be applied after compensation of mis-registration caused by
389
Fig. 1. The relationship between interferometric phase observations in L.O.S. dr and in azimuth da direction and ground deformation considering the right looking ascending satellite geometry.
the scene’s topography. This can be accomplished by geometrical transforms using precise orbit information and a Digital Elevation Model (DEM) as explained by Fornaro et al. (2005) and Reigber et al. (2006). 2.2. Retrieval of the 3D deformation vector Up to now, it has been explained how two-dimensional (2D) deformation vectors can be extracted by means of spectral diversity. It remains the question how to construct the complete 3D deformation vectors. For this aim, the relationship between offset maps produced in range and azimuth direction and the ground deformation has to be defined. The projection of a 3D displacement
− →
vector d = [dv dn de ]T onto line-of-sight direction and satellite heading direction can be expressed as (compare Fig. 1),
− →A d
r
− →A d
a
= δr + dv cos(θinc ) − sin(θinc )[dn sin(H ) + de cos(H )] = δa + dn cos(H ) − de sin(H )
(3)
where θinc and H denote the incidence angle and the heading angle
− →
(positive clockwise from the North), respectively, d = [dv dn de ]T indicates the deformation vector in vertical, north and east direction, and the superscripts of A shows ascending geometry. Since there are three unknowns for the displacement, another interferometric observation from a descending orbit is required: data from ascending and descending orbits have different L.O.S. and azimuth directions, providing two linearly independent measurements. With D indicating descending geometry, one obtains
− →D d
r
− →D d
a
= δr + dv cos(θinc ) − sin(θinc )[dn sin(H ) − de cos(H )] = δa + dn cos(H ) − de sin(H ).
(4)
Although these equations are directly solvable with simple linear equation techniques, this is not necessarily optimal when the observations are obtained with different geometry and have not the same phase accuracies. Instead, a Weighted Least Squares (WLS) approach can be used. Weighting is performed by the phase standard deviation, which is highly related to the interferometric coherence. In this way, observations with high accuracy (high coherence value = low standard deviation) get higher weights and will, therefore, have a stronger influence on the estimated
390
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
Table 1 Co-seismic interferometric pairs. Co-seismic Pairs
Baseline (m)
Descending
Ascending
03.12.2003–11.02.2004 03.12.2003–07.01.2004 03.12.2003–07.01.2004 03.12.2003–11.02.2004
16.11.2003–29.02.2004 16.11.2003–29.02.2004 16.11.2003–29.02.2004 16.11.2003–29.02.2004
Temporal resolution
0.6–2 570–2 570–30 0.6–30
70–105 35–105 35–70 70–70
Fig. 2. Block diagram of the implemented processing chain to retrieve the 3D deformation map. 58°E
unknowns and vice versa. The weighting matrix P, which is inversely proportional to the covariance matrix of observations, is defined as in the following: 29°30'N
(5) 29°15'N
1/σ
2 Da
.
σ = √
B
2 N B−b
1 − γ2
1 28°45'N
1
29°N
Here the subscripts of D and A refer to the descending and ascending geometries, respectively, while r and a refer to the observation at the range and the azimuth directions, respectively. I.e., σAr indicates the standard deviation of the observation at the ascending geometry in range direction. The phase standard deviations defined by Bamler and Eineder (2005) can be retrieved using
29°N
1/σ 0
2 Dr
0 0 0
59°E
03.12.2003 - 11.02.2004 0 30 km
29°15'N
2 1/σAa 0 0
0 0
58°45'E
28°45'N
0
58°30'E
29°30'N
2 1/σAr 0 P = 0 0
58°15'E
r p B b
πγ
0
(6)
where N is the number of independent samples used in estimation of the SD phase, γ is the interferometric coherence, B is the bandwidth of the full-resolution signal and b is the new bandwidth after bandpass filtering. It may noted that b = 3B has been reported as a optimum bandwidth ratio by Bamler and Eineder (2005) for the estimation of differential shifts. − → An expression for d , which minimises the estimation error, can be derived as in the following
58°E
58°15'E
58°30'E
58°45'E
59°E
Fig. 3. The coherence image retrieved from the dataset having a 2 m spatial baseline and 70 days temporal resolution.
where Q is the matrix including the coefficients in Eqs. (3) and (4) regarding the direction of the observations, ∆ is the vector including observations from Eq. (2) and P is the weighting matrix formed as in Eq. (5).
was employed using SRTM-3 DEM. In the backward geocoding approach, for each DEM element, the image pixel with the nearest range-Doppler coordinate is calculated (Reigber et al., 2006). As reported by Sansosti et al. (2006), the accuracy of 1/8th pixel can be achieved with backward geocoding using SRTM-3 in the case of ENVISAT data. Fig. 2 shows the processing chain to retrieve the deformation map from the SD technique. As it can be observed in Fig. 3, the effect of the earthquake is already visible in the coherence map, where the city center (approx. in the middle) appears decorrelated due to the severe building damages caused by the earthquake; while normally urban areas are known to show a very high correlation. Moreover, the earthquake’s fault line can also be detected in the coherence map.
3. Experimental results
3.1. Deformation map over the city of Bam
The proposed SD technique with the aim of detection of 3D deformation was tested on co-seismic interferometric ENVISAT ASAR images, which are summarised in Table 1, over the city of Bam in the south-eastern part of Iran, which suffered from a strong earthquake (M v = 6.5) in winter 2003/2004. It is worth again mentioning that the precise geometric coregistration of SAR data is a strict requirement for deformation monitoring. In order to perform a precise coregistration, backward geocoding
This section presents and discuss the detection of ground deformation based on DInSAR and spectral diversity methods from co-seismic interferometric data sets. For this purpose, 3 descending orbits and 3 ascending orbits have been processed. In Fig. 4(a), a descending co-seismic interferogram, depicting the slant-range displacement due to the earthquake is shown, generated from the 03 December 2003 (before earthquake) and 07 January 2004 (after earthquake) ASAR products. Because of
− →
d =
Q T P∆ Q T PQ
(7)
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
a
b
-180°
(degree)
180°
391
c
1.20
(m)
-1.30
-2.50
(m)
-2.90
Fig. 4. (a) Line-of-sight displacement map with differential interferogram. (b) Along-track displacement (azimuth direction) map with spectral diversity. (c) Line-of-sight displacement map with spectral diversity. Sub-look images have been obtained by selecting 1/3 of the full bandwidth.
Fig. 5. Profiles from azimuth offset map in the left–right direction in meters.
the large displacement, the deformation pattern appears wrapped and its gradient could not be measured directly by DInSAR due to phase unwrapping and noise problems in the area close to the fault. In this area, the extrapolation of fringes or similar techniques could be considered. From the same data pair, spectral diversity deformation maps in azimuth and range direction, respectively, were derived, shown in Fig. 4(b–c). In contrast to Fig. 4(a), the fault line, which is located in the north–south direction, can now easily be detected without any a priori knowledge. In order to derive a quantitative analysis, the deformations along six profiles, marked as black lines in Fig. 4(b–c), were analysed. As a result of the earthquake, surface movements of various amounts could be found (Figs. 5–7). A discontinuity and a sign change at the fault crossing can be found in azimuth profiles, see Fig. 5. An analysis of the azimuth displacements also shows that the amount of deformation decreases from north to south
direction. In the first horizontal profile, which is close to the city center, the estimated amount of azimuth deformation changes by about 70 cm from the beginning of the profile to the end of the profile. In the last azimuth profile, this change is only between −0.01 and −0.22 cm (21 cm along the profile line). However, the SNR (signal to noise ratio) in range is much lower than the one in azimuth. The discontinuity and the increase of deformation along the range profiles can been observed in Fig. 6. In addition to these ‘‘left-to-right’’ profile analysis, a ‘‘south–north’’ profile analysis was performed on the azimuth deformation map, shown in Fig. 7. This figure supports the assumption that both sides of the fault move away from each other: one of them tends more towards the north direction, the other more towards the south direction. This fact is for example not visible in a descending differential interferogram, because the surface motion happens almost perfectly perpendicular to the radar line-of-sight, where the radar has no sensitivity.
392
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
Fig. 6. Profiles from range offset map in the left–right direction in meters. Note that the displacement map in range direction is close to noise level, its weight is lower in 3D vector extraction than that in azimuth direction.
Fig. 7. Profiles from azimuth offset map in the south–north direction in meters.
After applying the processing chain expounded in Fig. 2 to each group of data sets (see Table 1), several SD deformation maps (after spatial averaging) were combined using weighted least squares approach to deduce the full 3D deformation field
− →
of the Bam earthquake. The unknown parameters, i.e., d = [dv dn de ]T , were not solved for each pixel but rather for pixels with sufficient interferometric coherence, i.e., γ = 0.7. With this criterion, only the pixels containing reliable phase information were used in the calculations. Delaunay triangulation was used for interpolation and visual display of irregularly gridded data considering the selected pixels. After triangulation, a (smoothed) view of 3D deformation and surface motion in all directions can been obtained for each pixel. Fig. 8(a) shows the vertical deformation, derived using the proposed approach, which fits very well with DInSAR. In this figure it can be seen how the vertical deformation changes its sign through the fault line. This is consistent with the DInSAR since, from west to east, the colour cycle has a negative phase value in the southern part of the earthquake, which shows that the surface
motion is towards the sensor. The case of northern part, colour cycle shows that the surface deformation is away from the sensor (see Fig. 4(a)). Fig. q 8(b) shows the magnitude of the horizontal deformation, i.e.,
− → 2
− →
dn + d2e .
3.2. Validation The accuracy of the estimated deformation is directly related to the estimated coherence of the interferometric pair, which indicate the phase quality. According to this, the accuracy assessment of the proposed technique was performed based on the coherence level of data set. Fig. 9 shows the standard deviation of the splitband interferometric phase in azimuth and in range direction and the interferometric phase as a function of coherence, for three multilook levels regarding Eq. (6) and Hanssen (2001, eq. 4.2.27) respectively. As expected, the accuracy of the interferometric measurement (without consideration of phase unwrapping) is more accurate than SD ones for an equal number of samples with
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
a
393
b 10
50
(cm)
(cm)
-10
0
Fig. 8. (a) Vertical and (b) horizontal components of surface displacement obtained from the data set in Table 1 having 0.6 m baseline with descending acquisition and 30 m baseline with ascending acquisition.
a
b
c
Along-Track offsets (meters)
a
0.6
b
0.6
0.4
Slant-Range offsets (meters)
Fig. 9. Standard deviations of the split-band interferometric phase in azimuth (a) and in range (b) direction and the interferometric phase (c) as a function of coherence, for three multilook levels.
0.4
0.2
0.0
data set with 0.6 m baseline data set with 570 m naseline
-0.2
0.0
data set with 570 m naseline
-0.2
data set with 0.6 m baseline
-0.4
-0.6
0.2
-0.4
0
50
100
150
200
250
samples
-0.6 0
50
100
150
200
250
samples
Fig. 10. The averaged profiles of displacement and their error bars in azimuth (a) and range (b).
fixed correlation. Moreover, it can be seen that the accuracy in range direction by SD is much lower than in azimuth direction by SD. With the threshold of the coherence of γ = 0.7 using 49 × 49 samples, the worst accuracy of 9 cm in azimuth and 16 cm in range direction have been predicted. In order to validate our results regarding the different combination of data sets, we plotted azimuth and range offsets and their error bars along the profiles in the area of maximum deformation. As expected and can be clearly seen from Fig. 10 that offsets and their variances obtained from data sets having various spatial and temporal baselines are different. Due to this the final
product, i.e., Fig. 8, has been retrieved from the ascending and the descending data set that have a maximum coherence. 4. Conclusions and discussions A comparison of spectral diversity and interferometry techniques for detecting of ground deformation has been carried out. Two-dimensional offset measurements based on spectral diversity technique from two complementary ascending and descending passes were combined using weighted least squares approach with the aim of extracting a fully 3D deformation map. In order to obtain
394
E. Erten et al. / ISPRS Journal of Photogrammetry and Remote Sensing 65 (2010) 388–394
a general view of surface deformation map, Delaunay triangulation and smooth interpolation were performed. It is demonstrated that the spectral diversity technique can not only be used for precise coregistration, but also for resolving 3D movement vectors in vertical, north–south and east–west directions. It is also demonstrated that 3D deformation map can be extracted from only two interferometric pairs, which is not possible by means of DInSAR. However, with spectral diversity technique, only decimeter accuracy is possible in case of ENVISAT, in contrast to millimeter accuracy in case of DInSAR. However, at least for analysing large co-seismic motion, this seems to be sufficient. Additionally, new sensors, like for example TerraSAR-X and RADARSAT-2, which provide a much higher spatial resolution, will significantly enhance the precision of the proposed technique. Acknowledgements This work has been supported by the German Academic Exchange Service (DAAD). The authors would also like to thank the European Space Agency (ESA) for providing the ENVISAT data sets of the Bam region. References Bamler, R., Hartl, P., 1998. Synthetic aperture radar interferometry. Inverse Problems 14 (4), 1–54. Bamler, R., Eineder, M., 2005. Accuracy of differential shift estimation by correlation and split-bandwidth interferometry for wideband and delta-k SAR systems. IEEE Geoscience and Remote Sensing Letters 2 (2), 151–155. Erten, E., Reigber, A., Hellwich, O., Prats, P., 2009. Glacier velocity monitoring by maximum likelihood texture tracking. IEEE Transactions on Geoscience and Remote Sensing 47 (2), 394–405.
Evans, A.N., 2000. Glacier surface motion computation from digital image sequences. IEEE Transactions on Geoscience and Remote Sensing 38 (2), 1064–1072. Fialko, Y., Simons, M., 2001. The complete (3-D) surface displacement field in the epicentre area of the 1999 Mw 7.1 Hector Mine earthquake, California, from space geodetic observations. Geophysical Research Letters 28 (16), 3063–3066. Fornaro, G., Manunta, M., Serafino, F., Berardino, P., Sansosti, E., 2005. Advances in multipass SAR image registration. In: Proceedings of IGARSS’05, Seoul, SouthKorea, 25–29 July, pp. 4832–4835. Hanssen, F.R., 2001. Radar Interferometry: Data Interpretation and Error Analysis. Kluwer Academic Publishers. Leberl, F.W., Maurice, K., 1994. Automated radar image matching experiment. ISPRS Journal of Photogrammetry and Remote Sensing 49 (3), 19–33. Madsen, S.N., Zebker, H.A., 1992. Automated absolute phase retrieval in acrosstrack interferometry. In: Proceedings of IGARSS’92, Houston, TX, 26–29 May, 2, pp. 1582–1584. Massonnet, D., Vadon, H., Rossi, M., 1996. Reduction of the need for phase unwrapping in radar interferometry. IEEE Transactions on Geoscience and Remote Sensing 34 (2), 489–497. Prats, P., Reigber, A., Mallorqui, J.J., Broquetas, A., 2004. Interpolation-free coregistration and phase-correction of airborne SAR interferograms. IEEE Geoscience and Remote Sensing Letters 2 (2), 206–210. Reigber, A., Erten, E., Guillaso, S., Hellwich, O., 2006. I.D.I.O.T.: a free and easyto-use software tool for DInSAR analysis. 26. Wissenschaftlich-Technische Jahrestagung der DGPF. (on CDROM), pp. 287–293. Sansosti, E., Berardino, P., Manunta, M., Serafino, F., Fornaro, G., Irea, N., 2006. Geometrical SAR image registration. IEEE Transactions on Geoscience and Remote Sensing 44 (10-2), 2861–2870. Sarti, F., Briole, P., Pirri, M., 2006. Co-seismic fault rupture detection and slip measurement by ASAR precise correlation using coherence maximization: application to a north-south blind fault in the vicinity of Bam (Iran). IEEE Geoscience and Remote Sensing Letters 3 (2), 1–5. Scheiber, R., Moreira, A., 2000. Coregistration of interferometric SAR images using spectral diversity. IEEE Transactions on Geoscience and Remote Sensing 38 (5), 2179–2191. Strozzi, T., Kuckman, A., Nurray, T., Wegmüller, U., Werner, C.L., 2002. Glacier motion estimation using SAR offset-tracking procedures. IEEE Transactions on Geoscience and Remote Sensing 40 (11), 2384–2390.