Refining Low-Quality Digital Elevation Models Using Synthetic Aperture Radar Interferometry by Michael S. Seymour B.Sc. Hons., Queen's University, 1988 M.Eng., McMaster University, 1991 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF Doctor of Philosophy in T H E FACULTY OF GRADUATE STUDIES (Department of Electrical and Computer Engineering) We accept this thesis as conforming to the required standard The University of British Columbia July 1999 © Michael S. Seymour, 1999 In presenting this degree at the thesis in partial fulfilment of University of British Columbia, I agree freely available for reference copying of this department publication or and study. of this his or her requirements that the I further agree thesis for scholarly purposes by the representatives. may be It thesis for financial gain shall not is permission. Department of glgOjegfl. AM» 6>*tvrf£ The University of British Columbia Vancouver, Canada Date DE-6 (2/88) AUG 10/0/0* faw-jfirttQirr an advanced Library shall make it that permission for extensive granted by the understood be for allowed that without head of my copying or my written Abstract Two-pass synthetic aperture radar (SAR) interferometry (InSAR) is a technique for processing the phase difference between coincident SAR images to obtain the range difference from the two radars to a common point on the earth's surface. The accuracy of the range difference measurement is in the order of one millimeter, and this range information can be processed to obtain digital elevation models (DEMs) of the surface topography. The digital processing required to make the D E M is quite complicated, mainly due to two factors. Firstly, the phase information is obtained from complexvalued data and therefore lies between — n and + 7 T whereas the complete phase information is needed. To obtain this, the phase must be "unwrapped" where the missing integer number of 2n are estimated for each data sample. Secondly, the geometry of the satellite passes relative to each other must be known to an accuracy of a few millimeters in order to obtain the surface height values to the required accuracy (about 10 m). Both of these steps require supplemental information and manual guidance to be performed correctly. Phase unwrapping is difficult because of noise and undersampling inherent in the measurements. The geometry estimates are difficult to make because the orbit is only known to an accuracy of a few meters, and the received phase data is a non-linear function of the satellite geometry. In the past, ii the geometry estimates have been made using known ground control points (GCPs), which requires a considerable manual effort, and has its own set of errors. The objective of this thesis is to use supplemental information in the form of a coarse D E M to make the InSAR processing more accurate and more automatic. We achieve this objective by developing a new algorithm which incorporates the coarse D E M directly into the processing stream, with the result that phase unwrapping and geometry estimation are performed accurately and reliably. In effect, the input D E M points serve as a large, dense set of GCPs. While the accuracy of each input D E M point is not very high, the large number of them provide adequate geometric accuracy, particularly as an automatic algorithm can register them directly to the radar data. There are two key steps in the new algorithm, which are interwoven in an iterative framework. First of all, the satellite geometry is estimated from the D E M and interferometric phase. This is done with a non-linear, iterative optimization algorithm without having to unwrap the phase. Avoiding phase unwrapping is important, as phase unwrapping errors can significantly bias the geometry estimates. Second, the input D E M along with the refined satellite geometry are used to create a model of the unwrapped interferogram phase that should be received from the two satellite passes. When this phase is wrapped, and compared with the measured phase, a differential interferogram is obtained which represents the difference between the coarse input D E M and the topography as measured by the satellite. This differential interferogram has a relatively low bandwidth, which means that it can be filtered and unwrapped reliably and accurately. Finally, the information in the unwrapped interferogram is used to refine the grid spacing and vertical accuracy of the coarse D E M . We have used mathematical analysis and simulation to develop the algorithm, to obtain statistical quality measures and to understand what system parameters iii affect the accuracy of the D E M results. We find that the main factors affecting accuracy are the interferometer's sensitivity of phase to height and the number of available D E M points, including the size and variability of the input DEMs' errors. We have successfully applied the D E M refinement algorithm to ERS Tandem Mission and RADARS AT-1 data. The generated InSAR DEMs had standard deviations of 12 to 20 meters compared to a control D E M with approximately 3 meters standard deviation. The output InSAR-enhanced DEMs had two to four times improvement in height accuracy compared with the input DEMs. In this way we have demonstrated that one can generate reliable estimates of topography for standard SAR scenes without having access to precision orbit data. Thus we have shown that the processing bottlenecks in dealing with repeatpass satellite InSAR data can be overcome, and useful topographic information can be obtained from the vast supply of existing InSAR data sets. iv Contents Abstract ii Contents v List of Tables x List of Figures xii Acronyms xviii List of Symbols xx Acknowledgements xxii Dedication 1 xxiii Introduction 1 1.1 Motivation for Research 1 1.2 Problem Definition 3 v 1.3 Thesis Scope 4 1.4 Methodology 5 1.5 Thesis Contributions 5 1.6 Organization of the Thesis 6 2 Interferometric SAR Background 8 2.1 Overview of SAR Interferometry 9 2.1.1 InSAR Phase Data 9 2.1.2 Interferometric SAR Coherence Magnitude 13 2.1.3 Summary 15 2.2 Interferometric SAR Data Model 15 2.3 Interferometric SAR Geometric Phase Model 17 2.4 Baseline Accuracy Requirements 24 2.5 Satellite InSAR Algorithms 32 2.6 Baseline Estimation 34 2.6.1 35 2.6.2 2.7 Baseline Estimation By Identification of Flat Areas Calibration Using Registration Parameters [1] 35 2.6.3 DIAPASON [2, 3] 36 2.6.4 JPL's Baseline Estimation Method [4] 37 2.6.5 Image Modeling Approach [5] 38 2.6.6 Review of Baseline Estimation 38 Summary 39 vi 3 D E M Refinement Algorithm 40 3.1 Overview of the D E M Refinement algorithm 40 3.2 Development of the D E M Flattening Algorithm 45 3.3 Development of the D E M Updating Algorithm 50 3.4 D E M Refinement Algorithm Implementation 51 3.4.1 D E M Flattening Algorithm Implementation Issues 51 3.4.2 Phase Unwrapping Implementation Issues 3.4.3 D E M Updating Implementation Issues 52 3.4.4 Computational Resource Requirements 52 3.5 3.6 Algorithm Performance Issues 55 3.5.1 D E M Flattening Algorithm Performance Issues 55 3.5.2 D E M Updating Performance Issues 66 3.5.3 Effects of D E M Trend Errors 67 3.5.4 Expected Baseline and Height Accuracy 68 Summary 72 4 Simulations of D E M Refinement Algorithm 4.1 51 74 Random D E M Error Simulations 74 4.1.1 D E M Flattening Simulations 78 4.1.2 D E M Updating Simulations 83 4.2 Effect of Trend Errors 89 4.3 Summary 92 vii 5 D E M Refinement Processing Experiments 5.1 Test Site Overview 95 5.2 D E M Data Overview 97 5.3 Experimental Procedure 103 5.4 ERS D E M Refinement Experiments 106 5.4.1 ERS D E M Flattening Experiments 110 5.4.2 ERS D E M Updating Results 122 RADARSAT D E M Refinement Experiments 130 5.5.1 RADARSAT D E M Flattening Experiment 134 5.5.2 RADARSAT D E M Updating Results 141 5.5 5.6 6 94 Summary of Experimental Results Conclusions 149 151 6.1 Thesis Contributions 153 6.2 Future Work 153 Bibliography 156 Appendix A Spectral Shift Theorem 164 Appendix B Outline of Satellite InSAR Algorithm to Obtain Topography 168 Appendix C Phase Unwrapping 175 viii Appendix D D E M Refinement Algorithm Details D.l Minimization Algorithm 182 182 D.2 Baseline Accuracy Requirements 185 D.3 Cramer-Rao Bounds for 2-D Frequency Estimates 187 D.4 Conservation of Mean Values in LS Problems 189 D.5 Plane Fits to Random Noise 190 Appendix E D E M Refinement Algorithm Definition 195 Appendix F ERS D E M Histograms 200 Appendix G R A D A R S A T D E M Histograms 210 ix L i s t o f Tables 3.1 Data stores required for the D E M refinement algorithm 3.2 Estimated mean and across swath height errors in planar fit to D E M errors for N = N = 100 r 4.1 53 70 a Statistics and lower bound on output D E M accuracy for degraded D E M height errors in meters 76 4.2 D E Mflatteningsimulation results 79 4.3 Input D E M error statistics for D E M updating experiments 86 4.4 Output D E M error statistics for D E M updating experiments 87 4.5 Output D E M error statistics using unwrapped phase with D E M flattening baseline parameters 88 4.6 Baseline estimates and height errors in the trend error experiments. 91 4.7 Baseline estimates and height errors in the trend error experiments. 91 5.1 Summary of D E M characteristics for the Chilcotin test site 99 5.2 Nominal ERS SAR parameters [6] 106 x 5.3 Summary of estimated baseline parameter errors contributions to output InSAR D E M accuracy for ERS data 107 5.4 Summary of results for ERS topography estimation 124 5.5 Nominal RADARSAT Fine Beam 3 SAR parameters 130 5.6 Summary of estimated baseline parameter errors contributions to out- 5J put InSAR D E M accuracy for RADARSAT data 131 Summary of RADARSAT-1 D E M error statistics 143 xi List of Figures 2.1 Phase data for the Death Valley RADARSAT-1 SAR image data. . . 11 2.2 RADARSAT-1 Death Valley interferogram 12 2.3 Coherence magnitude calculation 14 2.4 Two dimensional projection of interferometer geometry 18 2.5 Half-ambiguity height in meters as a function of normal baseline for ERS and RADARSAT 23 2.6 Effects of parameter errors on height estimates 27 2.7 Baseline length estimation requirements in meters for ERS data. . . 29 2.8 Baseline length estimation requirements in meters for RADARSAT-1 Fine 3 beam data 29 Sample ERS baseline parameter errors 31 2.10 Sample ERS baseline parameter errors 31 3.1 Algorithm for updating DEMs 42 3.2 Overview of the D E M "flattening" algorithm 43 2.9 xii 3.3 Overview of the D E M "updating" algorithm, following phase unwrapping 3.4 44 Example of residual phase term after initial estimate of geometry for 25m data set 3.5 49 Relationship between chip size, number of chips and accuracy of estimate of mean of 8(r) 3.6 58 Cramer-Rao bounds calculated for Nr — Na for the model parameters of the residual interferogram phase 3.7 60 Threshold coherence magnitude (//) and equivalent SNR for ./V = N =N 61 3.8 Sample spectra for D E M errors 65 3.9 Range of baseline parameter errors for GTOPO30 input D E M for a x y standard ERS scene 71 4.1 Pruned D E M examples 77 4.2 Comparison of residual interferogram spectra plotted in dB for azimuth,range, and azimuth range frequency components 82 5.1 Aerial photos of the Chilcotin test site 96 5.2 TRIM D E M of the Chilcotin test site. 100 5.3 DTED-1 D E M of the Chilcotin test site 100 5.4 GTOPO30 D E M of the Chilcotin test site 101 5.5 Zoomed plot of TRIM data 101 5.6 Zoomed plot of D T E D data 102 xiii 5.7 Zoomed plot of GTOPO30 data 102 5.8 Summary of ERS Tandem mission processed data set 108 5.9 Expected parameter errors for input DEMs for ERS data 109 5.10 Input DEMs for ERS interferograms 114 5.11 Results of the baseline estimation algorithm applied to ERS Tandem mission data 115 5.12 Model interferograms using estimated baseline parameters from ERS processing 116 5.13 Raw residual interferograms from ERS processing 117 5.14 Filtered residual interferograms from ERS processing 118 5.15 Comparison between DTED, trend corrected D T E D , and TRIM baseline estimates for the ERS dataset 119 5.16 Full resolution comparison of input TRIM D E M and filtered residual interferogram phase 120 5.17 Full resolution comparison of input DTED-1 D E M and filtered residual interferogram phase 121 5.18 Summary of baseline parameter estimates 125 5.19 Results of InSAR-TRIM D E M Experiment with ERS Tandem data. 126 5.20 Results of InSAR-DTED-1 D E M experiment with ERS Tandem data. 127 5.21 Results of InSAR-Modified DTED-1 D E M experiment with ERS Tandem data 127 5.22 Results of InSAR-GTOPO30 D E M experiment with ERS Tandem data 128 xiv 5.23 Results of InSAR-Modified GTOPO30 D E M experiment with ERS Tandem data 128 5.24 Errors of input and output GTOPO30 DEMs 129 5.25 Comparison of output modified InSAR DEMs with output InSARTRIM DEMs. Grey scale goes from -50 m to 50 m 5.26 Summary of RADARSAT processed data sets 129 132 5.27 Estimated parameter errors for input DEMs for RADARSAT data. . 133 5.28 Input DEMs for RADARSAT data 135 5.29 Results of the flattening algorithm applied to RADARSAT data. . . 136 5.30 Model interferograms using estimated baseline parameters from RADARSAT1 processing • 5.31 Residual interferograms from processed RADARSAT data 137 138 5.32 Filtered residual interferograms from processed RADARSAT data. . 139 5.33 Azimuth spectra of residual RADARSAT interferograms from each of the processed input D E M datasets 140 5.34 Summary of baseline parameter estimates for RADARSAT-1 data. . 144 5.35 Results of InSAR-TRIM D E M Experiment with RADARSAT-1 data. 145 5.36 Results of InSAR-DTED-1 D E M experiment with RADARSAT-1 data. 146 5.37 Results of InSAR-Modified DTED-1 D E M experiment with RADARSAT1 data 146 5.38 Results of InSAR-GTOPO30 D E M experiment with RADARSAT-1 data 147 xv 5.39 Results of InSAR-Modified GTOPO30 D E M experiment with RADARSAT1 data 147 5.40 Errors of input and output D T E D DEMs 148 5.41 Difference of output InSAR DEMs for RADARSAT-1 data 148 A.l C.l Diagram showing geometrical relation between ground wavenumbers and slant range wavenumbers 167 Example of residue calculation 176 C.2 Absolute phase surface with phase inconsistency 177 C.3 Wrapped phase and residues 178 E. l Top level block diagram for the D E M refinement algorithm 196 F. l InSAR-TRIM D E M error analysis for ERS Tandem data 201 F.2 Input DTED-1 D E M error analysis F.3 202 Output InSAR-DTED-1 D E M error analysis for ERS Tandem Mission data 203 F.4 Modified DTED-1 D E M error analysis 204 F.5 InSAR-Modified DTED-1 D E M error analysis 205 F.6 Input GTOPO30 D E M error analysis 206 F.7 Output InSAR-GTOPO30 D E M error analysis 207 F.8 Modified GTOPO30 D E M error analysis 208 F.9 209 InSAR-Modified GTOPO30 D E M error analysis xvi G.l Output InSAR-TRIM D E M error analysis for RADARSAT-1 data. . 211 G.2 Input DTED-1 D E M error analysis 212 G.3 Output InSAR-DTED-1 D E M error analysis for RADARSAT-1 data. 213 G.4 Modified DTED-1 input D E M error analysis 214 G.5 Output InSAR-Modified DTED-1 D E M error analysis for RADARS AT1 data 215 G.6 Input GTOPO30 D E M error analysis 216 G.7 Output InSAR-GTOPO30 D E M error analysis for RADARSAT-1 data.217 G.8 Modified input GTOPO30 D E M error analysis 218 G.9 Output InSAR-Modified GTOPO30 D E M error analysis for RADARSAT1 data 219 xvii Acronyms ATI along-track interferometry CCRS Canada Center for Remote Sensing CE90 circular error at 90% confidence CNES Centre National d'Etudes Spatial, French space research center DCW Digital Chart of the World DEM digital elevation model DFT discrete Fourier transform DLR Deutsches Zentrum fur Luft-und Raumfahrt, German aerospace research DTED digital terrain elevation data DOSAR Dornier airborne SAR DTM digital terrain model EMISAR Danish airborne remote sensing platform ERS Earth Remote Sensing satellite FFT fast Fourier transform GHz gigaherz, 10 Hz InSAR Interferometric Synthetic Aperture Radar, also known as IfSAR ISAR Inverse Synthetic Aperture Radar JPL Jet Propulsion Laboratory JERS Japanese Earth Remote Sensing satellite LE90 linear error at 90% confidence 9 xviii MHz megahertz, 10 Hz POLIMI Politecnico di Milano radar RAdio Detection And Ranging RADARSAT Canadian radar remote sensing satellite RMSE root mean square error SAR Synthetic Aperture Radar SLC single look complex-valued (image) SRTM Shuttle Radar Topographic Mission TOPSAR JPL airborne interferometry system and processor TRIM Terrain Resource Information Map TUD Technical University of Denmark TUDelft Technical University of Delft XTI cross-track interferometry 6 xix List of Symbols a satellite radial distance a ground slope HP) ground backscatter coefficient B baseline magnitude AB error in baseline magnitude, B BW parallel baseline, Bcoscf) B- perpendicular or normal baseline, B sin <j> c speed of light in a vacuum Ar the slant range sampling interval V along track time fo SAR center frequency fr the range fringe frequency term 1 Ah mean height error Ah ifi height error across the swath h(r) the ground radial distance H the satellite altitude above a reference plane at the Earth's surface I SAR image magnitude k 2TT/X, A SAR center frequency wavelength mean swa m v wavenumber orbit convergence rate XX m(r) an integer number accounting for the number of 2TT wraps P coherence magnitude w a azimuth radian frequency w range radian frequency W, azimuth-range cross fringe radian frequency r R A InSAR interior angle $ wrapped interferogram phase, n > ip > —x phase offset unwrapped interferogram phase, VP = ip + 2 m(r) n + ip 0 p scatterer coordinates SAR impulse response in range and azimuth Pg ground range sampling interval Pr slant range sampling interval r slant range from reference satellite antenna's phase center r2 slant range from candidate satellite antenna's phase center "(P) E [6(p)6(p)*], the effective volume scattering coefficient 0 off-nadir angle 0 baseline orientation with respect to nadir line AQ error in baseline orientation, 0 T range time X ground range c SAR image phase xxi Acknowledgements I would like to thank my thesis supervisor Dr. Ian Cumming for his SAR data processing knowledge and for providing me the possibility of exploring SAR interferometry in more depth. Drs. Pasquale, Prati, and Rocca of Politecnico di Milano provided me with a stimulating research visit that shaped my research and gave me a good introduction to the Italian language. Lab life was stimulating and enjoyable because of Gordon Davidson, Dave Michelson, Juan-Luis Valero and Sandor Albrecht and other lab members of the Radar Remote Sensing Group. Last but not least, Dr. Frank Wong for sacrificing the last days of summer to give me feedback on my thesis. My parents Diana and Vince have always supported my academic endeavors and have always had an open door when I needed to get away from it all. Enfin, il faut remercier ma femme Karine qui m'a attendu si long. Sans doute, sans elle je serais toujours en train de finir ma these. I would also like to thank MacDonald Dettwiler and the BC Science Council for financial support and the Canadian Space Agency and European Space Agency for supplying data used during the course of the thesis research. xxii For my father and friend Vince. xxiii Chapter 1 Introduction 1.1 M o t i v a t i o n for R e s e a r c h Synthetic aperture radar interferometry (InSAR) is a technique for deriving digital elevation models (DEMs) or surface motion maps (differential InSAR). The technique is based on processing the phase difference between coincident synthetic aperture radar (SAR) images and can yield accurate digital elevation models (for stable ground targets) [4, 7]. Digital elevation models are useful products for those dealing with topography such as urban planners, foresters, and remote sensing researchers. In SAR remote sensing, DEMs are used to resample the SAR images to well known coordinate systems and to correct for the distortions of the SAR imaging sensor [8]. There are many sources of InSAR data coming from a variety of sensor configurations (see [9, 10, 11, 12] for example). However, we focus on two-pass satellite interferometers as they are sensors with the most ground coverage and number of scenes currently available. In this form of interferometry, the scene of interest is imaged at least twice by the interferometer system at different times. 1 Despite the large amount of data, much of the two-pass satellite SAR data remains unexploited due to two major problems: 1. low SNR, 2. lack of precise interferometer geometry estimates. The low SNR is due mostly to temporal change between the two separate imaging passes [13] which adds noise to the interferometric phase. System limitations due to the need for operation in space (i.e. relatively low bandwidth compared with airborne InSAR sensors) also limit the utility of the data. In addition, two-pass satellite InSAR data can suffer from atmospheric artifacts [14, 15] that can be filtered if enough images are available [16]. The generally poor signal to noise ratio of two-pass satellite InSAR data can complicate the data processing significantly. One technique used to pre-condition the data for filtering is to "flatten" the data by removing a phase trend which models the response of the flat Earth in the interferogram phase. The resulting phase difference image looks much like a height contour map. Flattening the data facilitates low-pass filtering and subsampling of the interferogram to increase the SNR of the interferogram phase. In addition, the relative position of the two sensors during the imaging passes needs to be known precisely to produce accurate topographic data [4]. Currently, no satellite SAR system is equipped with sufficiently precise orbit data to produce accurate topography estimates directly [17]. In particular, Canada's RADARSAT1 SAR has.poor orbit data [18] which precludes accurate topographic estimation without further processing. Most production systems for topographic mapping rely on manually "tieing" selected identifiable points in the SAR image to points of accurately known position and height to estimate the interferometer geometry [4, 19, 20]. The requirement for manual matching of SAR image features and accurate 2 height control points and collection of elevation height control points themselves limits the feasibility of producing topographic maps in a production setting. Recently, a coarse, low-quality D E M has become available with coverage over most of the Earth's land surface. The GTOPOP30 [21] is currently available pub1 licly at minimal cost from the United States Geological Survey (USGS). GTOPO30 has low resolution with nominal data posting of 1 km and nominal 90% significance level of 160 meters. Further testing the accuracy of the GTOPO30 data set [21] suggests that the global RMSE accuracy of the dataset is actually « 70 meters; better than the nominal standard deviation of w 100 meters. However, the local quality of the GTOPO30 dataset varies as a function of the underlying data source and may be subject to outliers and local error trends [21]. The main theme of the thesis research is investigating the use of coarse, lowquality DEMs such as GTOPO30 in two-pass satellite InSAR data processing when there is no accurate satellite position data available. We want to take advantage of the large amount of low-quality D E M data available to ease the difficulties in two-pass satellite InSAR data processing as suggested in [22] and to minimize the need for manual intervention in the processing stream. 1.2 Problem Definition There are two obvious places where the input coarse D E M can be included in the processing: • "flattening", where a model of the interferogram phase derived from the coarse D E M is removed from the interferogram phase [23, 3]; • "updating", where the final estimate of InSAR D E M is made. 'See http://edcwww.cr.usgs.gov/landdaac/gtopo30/gtopo30.html. 3 The key to inclusion of D E M data in the processing stream is estimation of the geometry of the interferometric system (i.e. the baseline parameters) as the interferometer geometry defines the transformation between topography and measured interferogram phase. This problem is often complicated because of the lack of good orbit data and the low-quality of the input DEMs. It is important to note that in the flattening stage, we are dealing with the raw interferogram phase and do not have access to the unwrapped phase. In the updating step we must use the unwrapped interferogram phase in concert with the whole D E M to produce the output InSAR D E M . 1.3 Thesis Scope The purpose of this thesis is to develop, analyze, and test "flattening" and "updating" algorithms for inclusion of coarse low-quality D E M in the InSAR processing stream. The specific objectives are: • Develop a model of measured interferometric phase. • Review current methods of two-pass satellite interferometric SAR processing with emphasis on methods of baseline estimation. Develop the requirements for the baseline estimate accuracies. • Develop algorithms for "flattening" interferograms and "updating" DEMs. • Analyze the proposed algorithms' performances. • Implement the algorithms and apply them to simulated data, ERS Tandem mission data, and to RADARSAT-1 data. 4 • Analyze the results of the tests with SAR data and draw conclusions about the performance of the proposed algorithms with real data. 1.4 Methodology The "flattening" and "updating" algorithms were developed by expressing the geometry of the interferometer in terms of measurable parameters available at the time of processing. The baseline parameters are estimated in both the flattening and updating algorithm using a non-linear optimization scheme that fits measured InSAR data to geometrical models that are a function of baseline parameters. After establishing the algorithms, we analyzed the algorithms' performances in terms of baseline accuracy and height accuracy. The analysis was performed by examining the effects of parameter errors on the estimated baseline and topographic height. Special attention was paid to the effect that mean offsets and bi-linear trend errors in the input DEMs have on the output InSAR DEMs. Simulations were performed to verify the algorithms' analyses under controlled conditions. Finally, the whole D E M refinement algorithm was tested on satellite InSAR datasets from RADARSAT-1 and the ERS Tandem mission using three different input DEMs of varying qualities. The RADARSAT-1 dataset is an example of a poor satellite InSAR data set with relatively low SNR on average and some areas where there is no interferometric phase data at all. The RADARSAT-1 data also presents the additional difficulty of lack of accurate orbit data. In contrast, the ERS dataset has relatively good SNR and relatively precise orbit data. 1.5 Thesis Contributions The contributions of this thesis are: 5 • The development of a novel algorithm for two-pass satellite interferometric SAR processing using coarse, low-quality DEMs that consists of two subalgorithms: - A new algorithm for "flattening" the phase of an interferogram using the input low-quality D E M in the absence of accurate orbit data. - An algorithm for accurately "updating" the input low-quality DEM. • An analysis of the algorithms' performances which yielded - A method for determining when the baseline parameters used in the flattening algorithm are accurate. - An analysis of the effect of D E M errors on the determination of the interferometer geometry. - A better understanding of the aspects of parameter errors (including baseline parameter errors) with respect to the accuracy of output InSAR DEMs. - A method for calculating phase unwrapping group offsets to equalize independently unwrapped groups of interferogram phase. - A method for predicting the lower bound on topographic error due to baseline error when using coarse low-quality DEMs in the processing. • An analysis of the algorithm results for RADARSAT-1 and ERS Tandem interferometric SAR data and possible areas for improvement. 1.6 Organization of the Thesis The thesis is organized as follows: Chapter 2 gives detailed background information for the developed algorithm including 6 • mathematical relations between InSAR geometry and interferogram phase, • current difficulties with standard InSAR processing that would be alleviated by processing with a coarse low-quality D E M , • a review of baseline estimation techniques. Chapter 3 defines a method for incorporating coarse DEMs in the InSAR algorithm. A detailed analysis of the algorithm performance is presented. Chapter 4 reports the results of simulations that are used to verify the expected performance of the algorithms. Chapter 5 provides examples of the D E M refinement algorithm applied to ERS Tandem and RADARSAT-1 interferometric SAR data. The results are analyzed and discussed with respect to the performance predictions generated in Chapter 3. Chapter 6 concludes the thesis and summarizes the results of the work. Some possible directions for future work are given. Background material is given in the appendices as follows: • Appendix A describes the spectral shift phenomenon of InSAR data. • Appendix B gives an overview of conventional satellite interferometric SAR systems. • Appendix C describes the need for phase unwrapping and some of the common algorithms for performing it. • Appendix D gives detailed mathematical developments for the D E M refinement algorithm. • Appendix E gives a concise description of the D E M refinement algorithm. • Appendix F and Appendix G provide histograms of the errors of the DEMs generated during the experiments. 7 Chapter 2 Interferometric S A R Background This chapter gives a detailed review of the InSAR data model and current methods of InSAR processing. After giving an overview of SAR interferometry in Section 2.1, a mathematical model of InSAR data is reviewed in Section 2.2. The geometrical relationship between interferogram phase and topography is then derived in Section 2.3. Since the unwrapped phase to height conversion depends on the baseline parameters, we discuss the required accuracy for the parameters in Section 2.4. InSAR processing has been the focus of very active research in recent years and there is a number of different processing strategies used. These strategies are reviewed in Section 2.5. Special emphasis is placed on baseline estimation in this section. Finally, we summarize the background information from the view of using coarse DEMs in the interferometric SAR processing chain. 8 2.1 Overview of S A R Interferometry SAR interferometry (InSAR) is a technique for exploiting the phase difference between two registered SAR images generated from data collected from nearly the same or the same antenna phase center positions. Two types of information can be extracted from the images directly: • phase difference or interferogram phase, • magnitude of complex-valued correlation or coherence magnitude. The interferogram phase (assuming the images are registered and well correlated) is directly related to the difference in slant range distance from the two radar systems to the common ground patch imaged by the radar. The coherence magnitude is a sensitive indicator of the "similarity" of the two SAR images which takes into account both magnitude and phase information. In the following subsections, a brief overview of interferogram phase and coherence magnitude and the information they provide is given. 2.1.1 InSAR Phase Data The phase of a single SAR images also contains useful information which is not directly evident unless another correlated SAR image is available. Well-focused complex-valued SAR images contain phase and magnitude components due to the ground reflectivity and also a travel time phase component related to the distance of closest approach of the SAR platform and the ground patch. If a second correlated image which is coincident with the initial image is available, the travel time phase difference between the two images can be extracted by differencing their phases. An example of the SAR image phases is shown in Figure (2.1). Without more analysis, the phase of each image appears to be random noise. However, the 9 smoothed phase difference between these two SAR images shown in Figure (2.2) a) has a striking pattern which is directly related to topographic height (see Section 2.3). Because SAR amplitude images are sometimes sensitive to topography, one can also observe some similarity between interferogram phase and amplitude as shown in Figure (2.2). Graham [24] published the first well known work defining InSAR as a method for estimating topography from the phase difference between two SAR images. SAR interferometry can be used to create a digital elevation model (DEM) which enhances SAR image exploitation by facilitating: • correction of SAR imaging geometry by projecting the data to the ground range plane, • correction for radiometric variation in SAR images due to terrain slope, • geocoding of data sets for same and multi-sensor processing. Differential SAR interferometry [11] is an extension to SAR interferometry. In differential interferometry, the phase difference between two coincident interferograms is compared. If the target has changed positions between the two SAR imaging passes without significantly reducing the correlation of the two SAR images, the interferogram phase contains a term related to the amount of target motion. If the topographic interferogram phase can be subtracted from the interferogram, then the residual phase is related to the amount of target motion. The second interferogram may either be generated from an external D E M or another interferometric pair. Differential interferometry is extremely sensitive to terrain movement because the measuring stick is a fraction of the wavelength of the center frequency of the SAR. For current satellite SAR systems, this length is on the order of a centimeter. The sensitivity of differential interferometry coupled with the dense coverage of SAR data is providing new insights into geophysical processes [23]. This technique has 10 been successfully applied to observing earthquake displacement [25] and ice flow in ice streams and glaciers [26, 27]. (a) Phase of reference SAR image in rads (Sep. 24/96) (b) Phase of candidate SAR image in rads (Oct. 18/96) Figure 2.1: Phase data for the Death Valley RADARSAT-1 SAR image data. 11 (a) "Flattened", smoothed phase difference or interferogram between the phase images of Figure (2.1) in rads (b) Interferogram magnitude. Figure 2.2: RADARSAT-1 Death Valley interferogram and SAR image magnitude. 12 2.1.2 Interferometric SAR Coherence Magnitude Coherence is essentially a measure of the complex valued correlation coefficient between two datasets [28]. The maximum likelihood estimate of coherence magnitude assuming constant interferogram phase is [29]: -j* fie = !<=i I (2.1) N N where: p = the coherence magnitude, ip = the maximum likelihood estimate of interferogram phase, h,k e^''* = complex amplitude of the k th pixel of the i th image. The numerator of the coherence expression is a coherent sum of the interferogram data while the denominator is a measure of standard deviation of the two images when the images are assumed to have different standard deviations. The numerator may be represented as a vector sum as shown in Figure (2.3). Low coherence data tends to have vectors oriented in widely varying direction leading to a smaller magnitude of the complete sum and smaller numerator in the coherence magnitude calculation. For high coherence magnitude data the vectors tend to be oriented in the same direction leading to a larger net sum in the numerator of the coherence magnitude calculation. The coherence magnitude is a useful measure of sources of decorrelation such as: system noise [30, 31] Additive system noise is a cause of phase noise in the estimates. The majority of additive noise is usually receiver noise. 13 Im Im / Re Re High coherence Low coherence Figure 2.3: Schematic diagram of the numerator sum of coherence magnitude calculation. temporal decorrelation [13] To date in two-pass satellite SAR interferometry, the minimum time lapse between successive coincident data collections is one day. For changing landscapes such as forests, high correlation between the returns in the focused data is not likely. Some degradation of the interferometric phase estimates will occur over changing targets. spectral decorrelation [32, 33, 34] Spectral decorrelation is a function of the two-dimensional spectral overlap of the two coincident SAR images. In azimuth, the degree of overlap is controlled by the beam pointing directions of the radar systems. In range, the different cross-track positions of the two platforms also causes spectral mis-match between the SAR images. See Appendix A for more information. quality of focusing and registration [31, 35] The precision of registration and the quality of the SAR processing used to produce the SAR images for the interferogram affects the coherence magnitude. The variance of the interferogram phase estimate is a function of the coherence 14 magnitude of the data [32, 29] Thus, an estimate of coherence magnitude yields an estimate of interferogram phase variance. Since interferogram phase is related directly to InSAR height, the estimated coherence magnitude can serve as a quality control check on the terrain height estimates. The connection between coherence magnitude and the system parameters above can also be used to optimize interferometric SAR processing [35]. In addition, coherence magnitude has been shown to be useful for terrain classification [36]. 2.1.3 Summary SAR interferometry produces two useful quantities: interferogram phase and coherence magnitude. The unwrapped interferogram phase can be related through trigonometry and some external parameters to terrain height assuming reasonable data quality and stability of the scene. Such elevation models are a useful product for rectifying SAR images as well as for processing other remote sensing data. If the scene is slowly moving in a cohesive manner as in a glacier or an area where tectonic movement has taken place, the interferometric phase contains a component due to the surface motion which may be extracted to measure the surface velocity. Coherence magnitude is a measure of the correlation of the images used to estimate the interferogram phase. Coherence magnitude can be used as an input to a classifier as well as a data quality measure. The data quality can be used to estimate the accuracy of an interferometric SAR DEM. 2.2 Interferometric S A RD a t a Model A complex-valued SAR image (commonly called a single look complex or SLC image) can be modeled as the convolution between the system impulse response function 15 and the ground backscatter function: I( ,r)e = f q( -r ',T-T{p))b{p)e-^ ^dV' Jv K (2.2) foT v V ) where: I = the complex-valued SAR image, / ... dV' Jv = volume integral over the scattering elements, q(r), T) = SAR impulse response in range r and azimuth 77, f = SAR center frequency, b(p) = ground backscatter coefficient, and p = scatterer coordinates. 0 An interferogram is formed by a complex conjugate multiplication of two coincident SAR images. The expected value (E [•]) of the interferogram assuming statistically independent and identically distributed targets in the target volume is [37] E [7(77, r)e < /2(r;2, T2) e~ j = K2 f q(n - 77', T - Hp)) ?*(T?2 - T?', r2 Jv r2(p)) a(p) e ^/(r2(p)-T( )) P d T / / ^ where: .72(772, r2) e~^ = 2 the second or candidate SAR image, a(p) = E[6(p)6(p)T, the effective volume scattering coefficient. Assuming a constant propagation velocity c, one can convert the time delays (r, r2) into slant range (r, r2) from the antenna phase center to the common imaged ground patch: E [/(T?, r)e~ t /2(772, r2)e * ] = j / q{n Jv _ - r/, r - 2 Hp)) q*{v2 - 77', r2 16 r2(p)) a(p) e ]2k ( <P)- (p)W, r2 r (2.4) where: T 2r/c, 2r2/c, c speed of light, k 2TT/X, A radar wavelength. and The complex-valued interferogram is therefore an average of all scatterer responses weighted by the scatterer strength and the SAR impulse response function. The phase of the interferogram is the weighted difference of the travel time phase components in the two images, weighted by the impulse response. 2.3 Interferometric S A R Geometric Phase Model The interferometric SAR data model has two different interpretations which depend on whether or not a finite bandwidth (Appendix A) or single frequency approach is taken to characterize the InSAR data. The single frequency approach that directly connects the geometry of the satellite SAR systems used to collect the data with the interferogram phase is discussed below. The finite bandwidth approach reveals some of the inherent properties of InSAR data which lead to algorithmic steps in the InSAR processing algorithm and is discussed in detail in Appendix A. In this section, we develop a two-dimensional model of the connection between InSAR geometry and interferogram phase and then show how one can go from unwrapped phase to topography. One can model the interferogram phase as a function of the slant range difference between the two images. If we consider a single effective scatterer at some sampled slant range from the reference SAR antenna phase center (r) with a (most 17 scatterer B baseline magnitude 0 baseline orientation e off-nadir angle <> i interior angle r reference slant range r2 second slant range h ground radial distance a satellite radial distance S reference satellite S2 second satellite Origin (center of earth) Figure 2.4: Two dimensional projection of interferometer geometry. likely) different slant range (r2) from the scatterer to the second SAR antenna phase center, the phase difference between these two images is: (2.5) where = the absolute value of unwrapped phase. For a simple two-dimensional satellite geometry (see Figure (2.4)) 5(r) is 18 defined through the Law of Cosines as: 6(r) = r2 — r, = yV + B -2rB = ^/r + B - 2 r B cos(0 - 0(r)) - r, 2 cos(#(r)) - r, 2 2 (2.6) 2 where, for spherical geometry: r — the slant range of the reference image, = the slant range of the second image, = the baseline length, 0 = the baseline orientation with respect to nadir 9(r) = arccos r2 fa + r - h \ 2 a = h(r) = \ 2 2 , )' 2ar the satellite radius, the ground patch radius. The measured interferogram phase (ift) is a nonlinear function of the slant range difference S(r) from both sensors to the common imaged ground patch: -0 = arg(exp(jy5(r))). (2.7) where: ip = the wrapped interferogram phase, TT > tp > — w, A = the SAR wavelength. The arg(-) function maps the value of ^ S(r) to the principal values of (—7r,7r] and therefore defines the requirement for estimating the absolute number of 27T "wraps" or phase unwrapping. Under the assumption of no clock errors or atmospheric perturbations, 8(r) can be recovered from the phase difference between the two images after phase 19 unwrapping, scaling, and addition of a constant: S(r) = A[^( ) r + = m ( )27r+V ], r (2-8) 0 477 -^(r); (2.9) where: m(r) ip 0 = an integer number accounting for the number of 27T wraps, = phase offset, ip(r) — wrapped phase (TT, TT], <J/(r) = unwrapped phase. Inverting the above relationships, 5(r) can be used in turn to calculate the topographic radial distance h(r). Firstly, the interior angle <f> is calculated using the Law of Cosines (see Figure (2.4)): r COS<^> 2 2 B + _ r 2 2 2rB r + B - (r + 5) 2rB 2 substituting for r2 yields , 2 2 ,, expanding and simplifying ( 10 cos = Since h = (r + a -2racos6y , fB 2~rB~ -2rS-S— \ , , 6 = 2 Q-4, 2 2 /2 (2.12) 2 = arccos I the topographic height can be calculated as: (f> h = \r +a — 2 2 010 j ' \]Y 0 — arccos 'B „ \. -2rS-S 2rB 2 2ra cos 20 (2.1 2 /2 (2.14) ^ ' This final equation relates the unwrapped interferogram phase (\J/) through S(r) to the terrain height (h) of the target area and is the basis for cross-track interferometry. Note that the ground range position (x) of the target patch can also be calculated as x = r cos# = r cos(@ — S). (2.15) The along-track position of the satellite defines the other coordinate of the geometry. The derivative of the S(r) in range is interesting as the local interferogram frequency is a scalar multiple of it. This derivative allows one to relate height change to interferogram phase change. Evaluating the derivative with respect to the reference slant range r yields: o_m or = _ or ( 2 . 1 6 ) Expanding r2(r) and evaluating the derivative yields 8S(r) r — r2 — B cos S(r) dr r2 r B sin S(r) dO(r) r2 dr (2.17) The value of B cos <j>(r) is the projection of the baseline to the reference slant range direction and is very nearly equal to the difference r - r2. Cancelling the first term in the equation yields r B sin S(r) dO(r) r2 dr dS(r) dr (2.18) Substituting for the derivative of 6(r) and simplifying yields B sin S(r) dS(r) dr h(r) dh{r) + cos d(r) dr w r2 (1 - cos 0{r)) 2 1/2 w (2.19) which may be re-written as h(r) dh{r) J- ' « B sin Sir) dr Lar2sin^(r) dr 1_ r2tan#(r). (2.20) The unwrapped phase derivative is a scalar multiple of this equation. An important parameter from the geometry is the baseline magnitude projected to the direction 21 normal to the radar line of sight: B 1 = B sin (j>(r); called the normal or perpendic- ular baseline. This value multiplies two other terms: the scaled height derivative in the range direction and an approximately constant term. The approximately constant term comes from the increase in differential slant range as one moves along the Earth's surface. For a topographically flat scene, the raw interferogram phase is approximately a phase ramp with the alternating phase function of a single sinusoid. The nearly constant term is like the center frequency component of the interferogram signal. It is common practice for this term to be modeled and removed as part of the InSAR data processing algorithm; a process called flattening. In its simplest form, flattening uses an estimate of the predominant range frequency component in the interferogram to shift the center frequency of the interferogram to zero frequency. The scaled height derivative term contains information about topographic variation with respect to the "flat-earth" term. This term provides the "bandwidth" of the interferometric signal and acts to modulate the approximately constant frequency of the flat earth term. The spectral width of this term is a function of the variability of the topography and the size of the normal baseline. Generally speaking, larger perpendicular baselines mean more sensitivity to inter-pixel topographic height change and less sensitivity of estimated height to phase noise. For estimation of topography, it is generally desirable to have larger normal baselines. However, processing interferograms with larger normal baselines is generally more difficult than the smaller normal baseline case because the data is in general noisier (see Appendix A) and because the interferogram phase changes more rapidly requiring more work in the phase unwrapping processing (see Appendix C). For differential interferometry, small normal baselines are the norm. The sensitivity of a particular perpendicular baseline is usually described in terms of the ambiguity or half-ambiguity height. The height change causing ir 22 rads change in the measured phase is the half-ambiguity height. See Figure (2.5) for some sample half-ambiguity heights for different normal baselines for ERS and RADARSAT data. :::!:::::::::!:::::::::i::::::::i:::: ! i j J near range" ; ....i rrr-^r. m i d - r a n g e . \.....[ • : - far range ; : • \ v - - J l 0 i 50 i i 100 150 i I — 1 . 200 250 300 350 Normal Baseline in meters i 1 1 400 450 500 800 900 1000 (a) ERS Tandem Mission 0 100 200 300 400 500 600 Norma! Baseline in meters 700 (b) R A D A R S A T Fine 3 Beam Figure 2.5: Half-ambiguity height in meters as a function of normal baseline for ERS and RADARSAT. 23 2.4 Baseline Accuracy Requirements The desired performance of the baseline estimation algorithm depends on the application of the data. For differential InSAR using ERS data for earthquake and glacier measurement/monitoring, the accuracy of the ERS precision orbit product for estimation of the baseline is sufficient [17]. For higher accuracy differential applications such as tidal loading or monitoring of land subsidence where sub-centimeter accuracy is required, accurate external estimates of the baseline must be provided. Probably the most demanding application is estimation of topography because small errors in the baseline magnitude or orientation can result in substantial terrain height errors which vary with ground range [4]. In the following we calculate the sensitivity of the height estimates to baseline parameter errors to define the required accuracy of the baseline parameters for topographic mapping. The height error at any point may be calculated by taking a dot product between the gradient of the height function (see eqn. (2.14)) and the perturbations of the baseline magnitude (B), baseline orientation (0), altitude (a) and unwrapped phase (5) variables: AB dh dh dh dh Ah A0 (2.21) d~B <90 ~da, ~dSAa AS (2.22) Neglecting the cases where Sill d) rZi 0 as these cases have small perpendicular baselines which are more suited to differential interferometry as opposed to topographic mapping, the expression for the height error due to an error in the parameters is Ah = asm(d)(-B h y/1 + cos(<f>)r) ^ -cos(^) J3 a 24 rasin(fl) h A Q + ^ s w A „ ; "W; s i + h h - 2 AS. (2.23) cos(6) B 2 The expression for each error source can be simplified somewhat using the following approximations: • a/h ?s 1: the fractional difference between the satellite altitude a and the topographic height is on the order of 800 kms for satellite sensors while the radius of the Earth is on the order of 6600 kms. • y/1 — cos( 6) B = | B 2 x |: some of the error terms vary inversely with the modulus of the normal baseline. Using these approximations, one can get a better appreciation for the height error sensitivities of the Q, a, S variables as • ©: The error in the baseline orientation is approximately scaled by the ground . range variable [4] H«*(r). (- ) 2 24 Since the ground range distance is on the order of several hundreds of kilometers for satellite InSAR systems, the height error sensitivity to errors in the baseline orientation is large and requires precise calibration. • a: The error in altitude translates almost directly into errors in topography (2.25) as the numerator of expression is equal to the topography projected to the satellite radial direction. • 6: The differential for 5 is The sensitivity to errors in S varies across the swath as a function of the ground range variable scaled by the inverse of the normal baseline. 25 It is informative to consider the height errors generated for different combinations of AB errors for different combinations of 0 and B parameters. Graphs of these height errors are shown in Figure (2.6) for a baseline magnitude (B) of 1000 meters and baseline orientations (©) which produces normal baseline (fi ) values of x approximately 0, 707, and 1000 meters respectively. The interferometer parameters were modeled after the parameters of RADARSAT-l's Fine 3 beam. For B - fa 0, 1 the height errors are very large, ranging between -40 km and -90 km from near to far range. For B - fa 707 meters, the height errors are more reasonable, ranging 1 from approximately -307 meters to -309 meters. For B - fa 1000 meters, the errors 1 vary from approximately 1 to -1 meters from near to far range. The plotted errors suggest that height errors due to errors in B produce approximately linear trends and offsets in the output height. Given that the height errors due to baseline parameter errors are linear trends and biases across the swath, one must analyze the limits on baseline parameter errors by estimating the error in terrain at near and far ranges. One can start estimating the required accuracy of the baseline parameters under the assumption that the output mean error in height across the range swath is 0. 26 0 10 20 30 40 50 60 70 80 90 100 Range Index Ah for (f> fa 0 1 1 1 1 1 1 1 1 1 5 cm, B - 707 Ah lor AB = -3071 -307.2 \ ; 1 ; i > - \ - -307.4 -307.6 e ; -307.8 ^ : ; -308 ; '• '• '• | \ > >i -308.6 -308.8 -309' 0 \ ; iTT^^^^. g -308.2 -308.4 \ .>S4^.: ; ~ ~ - \ . 10 1 ; ; ; ; ; : i ; - 90 100 - rj>^..i : >....T^^^j .^...^^^.j • : \ ' 1 30 - ; ; ': 1 20 ; > : 1 ; 40 1 SO 1 80 1 70 1 80 Range Index Ah for AhtorAB <j> « 7r/4 = 5 cm, B 1 - 1000 Range Index Ah for <p « 7r/2 Figure 2.6: Effects of parameter errors on height estimates for different interferometer geometries (normal baselines). The interferometry model was generated for a flat scene using RADARSAT-1 Fine 3 beam parameters for A S = 5 centimeters. 27 . Assuming that the a and 6 variables are specified accurately and that the output height error is zero mean, the residual height error across the range swath of the data is approximately (see Appendix D.2): sin 6 (r - B ) a sin 6 „„ , » IX £ — — AB Aj> , 2 (sin^ (j)) ' (r - B cos <t>) B h 2 Ah swath 2 A swath , 2.27 n 3 2 where A<f> h is the elevation angle subtended by the range swath and the range swat varying parameters (8, <f>, r) are evaluated at the far range of the range swath. Given a desired error in height across the swath, the limits on baseline error can be calculated from this relation. The baseline magnitude error limits are shown in Figure (2.7) for ERS data and Figure (2.8) for RADARSAT-1 data using eqn. (2.27). The baseline error is calculated to limit the error across the swath to ± 5 meters. The results plotted are restricted to normal baselines between 50 meters and 500 meters for the ERS case and 50 and 2500 meters for the RADARSAT-1 case. These maximum perpendicular baselines are calculated as approximately half the maximum perpendicular baseline [10] for flat terrain. Generally the requirements on the baseline magnitude accuracy are more stringent for the ERS case than for the RADARSAT-1 case. The previous analysis of height error due to baseline parameter error allowed us to eliminate one error variable because of the assumed condition of zero mean height error. However, a more complete characterization of the limits on baseline parameter accuracy can be constructed by considering the range of baseline parameters which give a D E M which falls within a specified range of both mean error and error across the swath. One can analyze the general case using the height error at near and far ranges respectively: Ahf , Ah . ar near 28 Limits o n Baseline Error Li - i 3 - 2 i - i 1 6 B i 0 in r a d s i 1 L 2 3 Figure 2.7: Baseline length estimation requirements in meters as a function of baseline length and baseline orientation for nominal ERS parameters. Output D E M has zero mean height error and a maximum height error of ± 5 meters across 100 km ground range swath. - 3 - 2 - 1 0_ i n r a d s 0 1 2 3 Figure 2.8: Baseline length estimation requirements in meters as a function of baseline length and baseline orientation for nominal RADARSAT-1 Fine 3 beam parameters. Output D E M has zero mean height error and a maximum height error of ± 5 meters across 55 km ground range swath. 29 We can define two performance criterion: • error across the swath, Ah th swa — Ahj Ah ar near which is a measure of the slope error of the output InSAR D E M , and • mean D E M error », _ Ahf ar ^"-mean — + 2 Ah near ' which approximates the mean D E M error with average of the two swath errors. Given a specified range of mean height errors and slope errors across the swath, a region of the baseline parameter space where the output InSAR D E M has the desired quality can be defined and then plotted. Figures 2.7 and 2.8 are examples of this type of plot for ERS standard images and RADARSAT-1 Fine 3 beam images. Both simulations assumed a normal baseline of approximately 200 meters with nominal parameters for each type of images (see Table 5.2 and Table 5.5). The dark region of these plots denotes the acceptable baseline parameter errors to limit the across swath error (Ah th) swa and the mean error (Ah ) mean to ± 5 meters to less than 1 meter. Roughly speaking, the baseline magnitude errors for the ERS data set are limited to ± 0.05 meters while the orientation errors are limited to ± 10 /^rads for this set of D E M error constraints. For the RADARSAT-1 data, the range of valid parameters is at least twice as large but is much "thinner". 30 -0.1 -0.08 -0.06 -0.04 -0.02 §5 E 0 AO in radians Figure 2.9: Range of baseline parameter errors for maximum 5 meter error between far and near range and maximum mean error of 1 meter for ERS parameters. Figure 2.10: Range of baseline parameter errors for maximum 5 meter error between far and near range and maximum mean error of 1 meter for RADARSAT-1 Fine 3 beam parameters. 31 The differences in the size and shape of the areas of valid baseline parameters come from the different range of off-nadir angles and hence ground ranges that parameterize the two instruments RADARSAT-1 Fine 3 beam interferometry is about four times as sensitive to baseline orientation errors as the ERS data because the ground range distances are much larger. For a given baseline magnitude error in the RADARSAT-1 dataset, the feasible range of baseline orientation errors are smaller than that of the ERS data which leads to "thinner" area of valid baseline parameters. However, the RADARSAT-1 swath width is half that of the ERS data so the the swath error criterion is less restrictive for the RADASAT-1 parameters. Consequently the range of valid parameters is larger though more tightly constrained in the RADARSAT-1 case than for the ERS case. 2.5 Satellite I n S A R Algorithms The basic processing for SAR interferometry can be described as follows (a concise description for topographic estimation is given in Appendix B): 1. Identify suitable correlated SAR signal data sets. 2. Process the SAR data to yield complex images with good phase fidelity. 3. Register the overlapping imagery to an accuracy of about 1/8 or 1/10 pixel. 4. Form the interferometric image by multiplying one image pixel for pixel with the complex conjugate of the other registered image. 5. Remove a model of the flat Earth interferogram phase. 6. Generate coherence magnitude estimates. 7. Process the interferometric phase to yield the desired data product: topographic height or surface motion estimates. 32 The removal of the flat Earth model of interferometric phase can be performed in a number of different ways (e.g. orbit modeling [23, 3] or demodulation by largest frequency component). In general, it is an approximate correction which reduces the local phase variability in the interferogram to facilitate processing such as coherence magnitude estimation, filtering (and hence phase unwrapping) and visualization of the topographic phase. These first stages of processing are standard and have been implemented in Matlab for research purposes. The last stage of information extraction is where the bulk of the recent research in SAR interferometry is on-going. For differential SAR interferometry, processing the interferogram phase tends to be simpler compared to the topographic height estimation problem. Generation of the surface velocity map involves: 1. Removal of topographic phase. 2. Filtering the residual interferogram phase. 3. Projection of unwrapped interferogram phase to surface motion vectors. Usually, though not always, differential InSAR datasets are chosen with a small normal baseline to minimize the interferometric data's sensitivity to topographic height variation. In this case, the flattening procedure above in the basic processing model yields phase fringes due to motion. If the test site happens to be flat, removal of the topographic fringes again can also performed by the flat Earth correction. However if an existing D E M is available and the topography is variable, the D E M can be used to model the geometric phase based on precision orbit data [17] or using other post-processing techniques [3]. Finally, since the radar interferometer only observes data in the range dimension, the data must be transformed to put the data in the ground coordinate system [27]. To generate InSAR topographic maps, the following steps must be implemented: 33 1. Unwrapping the interferogram phase. 2. Estimating the baseline parameters or geometry of the interferometer. 3. Resampling the computed heights to a standard cartographic representation. The effects of local interferometric fringe frequency, system and decorrelation noise and terrain discontinuities in the interferometric observations require fairly sophisticated processing to maintain the phase fidelity of the data for successful phase unwrapping. Phase unwrapping is a field of active research and there are currently many possible phase unwrapping methods (see Appendix C for a short review). Some methods combine the filtering and phase unwrapping procedure [38] while others use multiple interferograms [39, 40]. It has been noted that reducing the local frequency content in the interferogram [22, 41] reduces the difficulties with phase unwrapping by reducing the number of residues (see Appendix C). In addition, the step of pre-processing with an existing D E M is also seen as one way to reduce phase unwrapping difficulties [22, 3]. 2.6 Baseline Estimation There are a number of different strategies for estimating the interferometric SAR baseline parameters. The strategy used depends on the data available (i.e. D E M or height ground control points), application and required accuracy (i.e. differential InSAR or topographic InSAR), and variability of the terrain. In some implementations, the baseline estimation procedure is a part of the overall processing while in others it is distinct stage of processing. In the following, we review the main implementations of baseline parameter estimation methods. After reviewing the methods, we discuss them in the context of low-quality D E M data and poor orbit data. 34 2.6.1 Baseline Estimation By Identification of Flat Areas For a flat area of terrain, the local interferogram fringe frequency rate may be calculated as [42]: _ 2Bsin<?!> A r ^ ~ A#tan6> ' where: f r Bsmcf) 0 = the range fringe frequency term, = the normal baseline, = the ofT-nadir angle, A = the wavelength of the center frequency of the SAR, H = the satellite altitude above a reference plane, Ar = the slant range sampling interval. This method relies on the assumption that areas of flat topography exist and that they are identifiable in the image. In general, these are not good assumptions to make so this particular method is an approximation that is usually made for the initial flat Earth correction procedure. 2.6.2 Calibration Using Registration Parameters [1] At Politecnico di Milano (POLIMI), the registration parameters of the scene along with a local spherical model of the earth are used to estimate the baseline orientation and magnitude by considering the difference in ofT-nadir angle between the two data acquisitions. The look angle difference is a function of normal baseline and slant range distance. Because this method's baseline parameter accuracy is dependent on the accuracy of the registration, orbit accuracy and knowledge of the local topography, it is useful only for preliminary baseline parameter estimates that must be refined by other means. 35 2.6.3 DIAPASON [2, 3] CNES's DIAPASON InSAR processor is designed to perform differential InSAR images automatically but does not perform phase unwrapping. It implements the following procedure: 1. The external D E M and the SAR image orbit data are used along with interimage offsets estimated by image chip correlation to define the resampling grid required to make the initial interferogram. This provides an initial estimate of the baseline parameters for the scene through modeling the orbits of the satellites. 2. One of the radar images is registered to an absolute geographic reference by simulating a radar image from the input D E M as a function of the D E M slope. The simulated image is then correlated with the SAR image to register the SAR data with a known output grid. 3. Spectral shift filtering (see Appendix A). 4. Subtraction of the D E M based model of the interferogram phase. 5. Resampling of the interferogram into a usual geographic coordinate system. Running automatically, the DIAPASON'S output D E M will have residual fringes due to the accuracy limits of the initial orbit modeling attempt [2]. A post-processing step is then required to complete the modeling of the baseline geometry. However, CNES has developed a tool external to DIAPASON that positions the satellite tracks precisely using four ground control points and a "count" of the number of fringes between the pairs of points [2]. The focus of CNES's work has been on differential InSAR so their processing algorithm naturally includes an accurate D E M and accurate orbit data as part of the 36 processing stream. It is not clear what topographic accuracies are achievable with this processing method and how well the method works when the input D E M is of poor quality. In particular, there is the implied requirement for phase unwrapping in the fringe counting procedure and manual intervention in the choice of control points. Note that CNES did publish a method of simplifying phase unwrapping that used the DIAPASON processor [40]. The method uses integer interferometric combination (IIC) of interferograms and involves differencing flattened interferograms that have had their phase values multiplied by integers. The net effect is to create an interferogram with a large ambiguity height that is (subject to noise considerations) easier to unwrap. 2.6.4 JPL's Baseline Estimation Method [4] JPL's approach to baseline estimation algorithm is an empirical method that requires ground control points and unwrapped interferometric phase. They start from the expression relating unwrapped phase to topography for parallel orbits. However, in general satellite orbits converge or diverge slightly so that the interferogram phase appears to have a component that varies with azimuth position. They then modeled the slant range difference due to topography as 6 = ^[y-m^ + Ve]] (2.28) where m n = orbit convergence rate. (2.29) They then fixed the vertical component of the baseline based on knowledge of ERSl's vertical decay rate and performed an optimization using height ground control 37 points for the cross-track component of the baseline, the baseline orientation, and the constant phase offset. Using this technique, they were able to generate fairly accurate spot heights (four errors in the region of 1 to 5 m and one large outlier at 31 m). Small [43] reports that with 18 ground control points, an InSAR D E M with RMS error of 8.2 meters was attained for a standard ERS quarter scene. A similar procedure was reported for a SIR-C differential InSAR experiment in [23]. In this experiment, an orbit model was used to firstly to flatten the interferogram phase using an external D E M . The residual interferogram phase was filtered, downsampled, unwrapped and used to refine the orbital baseline estimate to estimate the local scene displacement velocity. 2.6.5 Image Modeling Approach [5] In the image modeling technique, the conversion from unwrapped phase to height is assumed to be a quadratic function: h = k + k {V + V ) + k {V + * ) . (2.30) 2 1 2 0 3 0 The coefficients for the quadratic function are estimated by fitting the quadratic relation to a calculation of the unwrapped phase based on orbit models alone. In principle, this method can provide reasonably accurate height estimates with only one ground control point that establishes the unwrapped phase offset (^ )0 In practice, more ground control points are needed to compensate for the orbit data errors. About 30 ground control points per full frame were required for a mapping exercise performed for the Czech Republic [44]. 2.6.6 Review of Baseline Estimation The methods for baseline estimation reviewed above vary in their complexity and how exactly they are applied. One common element is that the methods for gen38 erating accurate height estimates require using unwrapped phase and some sort of explicit manual interaction to generate accurate baseline parameters. If the only information that is available is a coarse low-quality D E M , identification of individual ground control points is unlikely. In addition, a small number of low-quality ground control points are unlikely to provide accurate topography estimates. Another consideration is the need for unwrapped phase in the estimation problem. Phase unwrapping errors [23] also have the potential to cause biases in the output height estimates. It is therefore useful to consider estimating baseline parameters by other means not using unwrapped phase. 2.7 Summary We have reviewed an InSAR data model and data processing at a high level to describe the technology and current state of the art. The baseline magnitude and orientation parameters are important for InSAR processing because they relate the measured interferogram phase to topography. They must be accurately estimated to provide reliable accurate topographic height estimates. To date, the only thoroughly tested method to obtain accurate baseline estimates for topographic estimation is the method using existing height control points. Existing height control points must be identified in the reference SAR image to calculate the baseline [4, 44, 3] using the unwrapped phase of the interferogram. Orbit data provides baselines which are generally not accurate enough for topographic height calculation [17, 18] for both RADARSAT-1 and ERS data. The CNES and POLIMI methods yield reasonable estimates of the interferometer geometry which give (with some manual intervention) interferograms with no residual phase trends. In the next chapter we present algorithms for flattening and baseline estimation that do not depend on accurate orbit data or manually choosing height ground control points for baseline estimation. 39 Chapter DEM 3 Refinement A l g o r i t h m This chapter describes the D E M refinement algorithm and analyzes the algorithm performance. We start by giving an overview of the method in Section 3.1. The method has two novel sub-algorithms: D E M flattening and D E M updating. The D E M flattening algorithm is developed in Section 3.2 followed by the D E M updating algorithm in Section 3.3. Some implementation details of the algorithm are discussed in Section 3.4 (a concise and complete algorithm description is given in Appendix E). Section 3.5 discusses issues affecting the algorithm's performance. We predict the accuracy limits of the topographic estimate in Section 3.5.4. Some concluding remarks are given in Section 3.6. 3.1 Overview of the D E M Refinement algorithm The algorithm for D E M improvement using InSAR techniques consists of three parts (see Figure (3.1)): 1. D E M "flattening" using a coarse, low-quality DEM; 2. Phase unwrapping of the flattened phase signal; 40 3. D E M "updating" using the unwrapped flattened phase. The main focus of this research is the first and third stages of the D E M refinement algorithm that are described in detail in Section 3.2 and Section 3.3 respectively. In the flattening algorithm, we attempt to generate an interferogram model from the input coarse D E M that is as close as possible to the input interferogram (see Figure (3.2)). The residual interferogram formed by the difference between the input interferogram and interferogram model can be filtered and downsampled to facilitate phase unwrapping processing. Finally, the unwrapped residual interferogram phase is then used by the D E M updating algorithm to find the best fit InSAR D E M with respect to the input D E M (see Figure (3.3)). Both the flattening and updating algorithms rely on estimation of the baseline parameters as these parameters relate topography and interferogram phase. In the flattening algorithm, the interferogram phase is generated from the input D E M using baseline parameters. In the updating algorithm, the InSAR D E M is generated from the unwrapped phase using the baseline parameters. In both algorithms, the baseline parameters are estimated using non-linear optimization. In the following two sections, we develop the two algorithms in detail. 41 Interferogram DEM Flattening Residual Interferogram Input Phase Unwrapping Coarse DEM Unwrapped Phase Refined DEM Figure 3.1: Algorithm for updating DEMs. 42 Low-quality DEM Interferogram Figure 3.2: Overview of the D E M "flattening" algorithm. 43 Unwrapped Phase Low-quality DEM Figure 3.3: Overview of the D E M "updating" algorithm, following phase unwrapping. 44 3.2 Development of the D E M Flattening Algorithm The D E M flattening algorithm tries to create an accurate model of the interferogram phase using the input D E M . Calculation of the interferogram phase model using the input D E M requires several parameters (see eqn. (2.6)): • slant range (r), • baseline orientation ( 0 ) , • baseline magnitude (B), and • reference sensor altitude (a). The reference slant range is given by the reference SAR sensor sampling times and for the moment, we assume that the sensor altitude is accurate and available. Therefore generating the interferogram phase model requires estimation of the baseline orientation ( 0 ) and baseline magnitude (B). The baseline parameter estimation algorithm uses a non-linear least-squares fit between a geometric model of the interferogram phase and a parametric approximation to the interferogram phase itself. This fit proceeds without phase unwrapping and therefore avoids biases due to possible phase unwrapping errors [23] and facilitates further processing of the interferogram such as straight forward lowpass filtering. It is an iterative procedure whereby erroneous phase trends in the residual interferogram are used to update the interferogram model of phase. The general idea of the baseline estimation algorithm is to fit S((B, 0 , a, h, r) = r2(J3, 0 , a, h,r) — r calculated from the input D E M , satellite altitude and the baseline parameters to an estimate of 6(r) derived from the interferogram: min(B, 0 ) £ > ( r ) k f c - (r2 {B, 0 , a, h, r) - r ))\ k 45 k (3.1) where k is an index running over all available samples. By applying this minimization equation, we are implicitly assuming Gaussian independent statistics for the differences in the fitted terms. This approach tries to capture the global information of the D E M rather than looking at local values of the unwrapped interferogram phase as with most of the baseline parameter estimation methods. The values of B and 0 which minimize eqn. (3.1) are found using the Levenburg-Marquardt algorithm [45] (see Appendix D.l). For the case of topographic estimation, we do not know the value of 5(r) a-priori. We are therefore reduced to using an approximation to S(r) which is updated as the algorithm progresses. As a starting point, 8{r) is estimated from the registration relation between the two SAR images used to form the interferogram. One can usually register the images assuming that there is slight rotation of the image (a range dependent azimuth shift) as well as a linear or quadratic range stretch. The range stretch component gives an approximation of S(r) which can be used in eqn. (3.1) to give initial estimates of the baseline parameters: B and 0 . The baseline parameters estimated from the minimization algorithm are used to form a residual interferogram (ig id al) res U which is (neglecting the interferogram amplitude and interferogram noise): igreMval = e * * *> e ~* * ^(B,%^)-r), ( 3 . 2 ) where: B,Q = values from the optimization algorithm. Assuming correct estimation of baseline parameters and a reasonable quality input D E M , the residual interferogram spectra will be centered at zero frequency as the constant range frequency term of the original interferogram is removed (see page 22). In addition to shifting the spectrum of the interferogram to lower frequencies as conventional flattening does, removing the geometric phase model also makes 46 the residual interferogram's bandwidth smaller than the original interferogram as a function of the input DEM's quality. If the input D E M had no errors, all of the energy in the bandwidth of the original interferogram would be concentrated at 0 frequency. However usually, there are both baseline parameter errors and D E M errors to contend with in the residual interferogram. The ideal noise free residual interferogram consists of two phase terms: igresidual = J*"™ J**™; (3.3) where: ^dem comes from height errors in the input D E M (Ah) and is approximately given by: An B Ah 1 ^dem « 7 ,„ . (3.4) A , x A where J3 is the normal baseline and x is the ground range distance. X In general, ^dem are correlated variables due to interpolation of the input coarse D E M and correlation in the topographic samples under scrutiny. tygeom is the residual phase term due to errors in the baseline parameters and is approximately 47T Vgeom « — (-cos <f>(r) AB + B sin <f>(r) A G ) , « y (~~p~^B + B AQJ . X (3.5) Since the off-nadir angle subtended by the range swath is quite small for satellite SAR systems ( only « 6° for both RADARSAT and ERS standard beam images), one can model geometric phase error term as the sum of two distinct linear trends; one from the change in cos<f>(r) across range, and the other from the change in sin^>(r). Therefore, the geometric phase error term may be approximated by: ^geom ~ Wr,a(V ~ Vo) {r ~ T ) 0 + W (r r 47 - T) 0 + W (t] a - T} ) 0 + Vo, (3-6) where: w, r w ,w a range, azimuth, and cross fringe frequency terms in radians/s, = rta T — T = 0 7] — T] 0 ip 0 = — range pixel time centered at 0, azimuth pixel time centered at 0, constant phase value, (—7r,7r]. An example of residual phase after estimation of geometry using only the registration relation is shown in Figure (3.4) (a constant baseline was used in this simulation). Note that for this example, the phase variations are quite small in large parts of the residual interferogram - corresponding to areas where the D E M fits the gentle topography well. Areas where the topography has higher frequency or steep components are not fitted as well by the D E M (or layover is present), and large phase fluctuations are noted. The residual phase trends and phase offset of the residual interferogram therefore contain information about the error in the estimate of <5(r) derived from the interferogram registration relations. The slope of the ramp or frequency of the residual interferogram may be efficiently estimated using a D F T followed by a golden section search to accurately estimate the location of the peak. Note that the baseline geometry usually varies in azimuth so care must be taken to estimate the trend in range or azimuth frequency as a function of range or azimuth to estimate residual phase. The geometry phase term can then be added to the initial estimate of the geometry to provide a more accurate 8{r) term for the fitting procedure: S{r~) = (r2(B,Q,a,h,r)-r) + T- (w ,a(v ~ Vo)(r - r ) r 0 + w (r r - T ) + w (n - Vo) + i>o)0 a (3.7) The next estimate of the geometry is made by solving the minimization problem eqn. (3.1) using the updated estimate of 6(r). The process is iterated until the frequency estimates converge to values of 0 or a maximum number of iterations is reached. 48 In practice, we have found that five or six iterations provides good estimates of the baseline geometry for reasonable quality (i.e. DTED-1) DEMs. Range Index Figure 3.4: Example of residual phase term after initial estimate of geometry for 25m data set. 49 3.3 Development of the D E M Updating Algorithm After filtering and unwrapping the residual phase from the D E M flattening algorithm, a refined estimate of S(r) can be made by adding the scaled unwrapped phase to the estimate of S(r) derived in the D E M flattening portion of the algorithm S(r): 6$) = 6$) + (3.8) where: * = the unwrapped phase . The updated value of S(r) can be used with the baseline parameters derived during the D E M flattening process (5,0) and eqn. (2.14) to derive refined topographic heights. However, the baseline parameters derived from the D E Mflatteningalgorithm most likely will need to be re-estimated . One can re-estimate the baseline param1 eters by minimizing the difference between the input coarse D E M and the D E M computed from the unwrapped phase as a function of the baseline parameters [46]: mintB, O) ^ ( / j ( O E M ) - h(B, 0 | 6(7), a, r)) . 2 (3.9) where: h(DEM) h(B,Q | S(r),a,r) = the input coarse DEM, — the D E M calculated from the unwrapped phase. This function also assumes that the errors between the estimated and D E M height are also independent and Gaussian as with eqn. (3.1). The nonlinear minimization in this baseline estimation problem is handled using the Levenburg-Marquardt algorithm as with D E M flattening procedure. Here, we take advantage of the dense 'See Section 3.5.1 for details. 50 network of height ground control points afforded by the input coarse D E M to constrain the InSAR topographic height estimate. An analysis of the performance of this iteration is given in Section 3.5.2. 3.4 D E M Refinement Algorithm Implementation The previous section's development of the D E M refinement algorithm omitted some implementation details for clarity. We review implementation details specific to the D E M Flattening algorithm, the phase unwrapping processing, and the D E M updating algorithm in Section 3.4.1 through Section 3.4.3 respectively. Implementation details including a summary of the computational requirements of the algorithm are given in Section 3.4.4. See Appendix E for a concise overview of the entire algorithm. 3.4.1 D E M Flattening Algorithm Implementation Issues There is a requirement to resample the existing coarse D E M into radar coordinates. If accurate orbit data were available, this procedure could be performed in a wellknown fashion [8]. In the absence of accurate orbit data, an initial registration step between the D E M and the SAR data must be completed. For RADARSAT-1 data, where the accuracy of the orbit data is not sufficient to register the input D E M data to SAR data, we have found that a registration step which maximizes the coherence magnitude of the interferogram as a function of baseline parameter estimate provides reasonable registration between the InSAR data and the D E M . 3.4.2 Phase Unwrapping Implementation Issues We use standard phase unwrapping algorithms (see Appendix C and Chapter 5). However, we also assign each valid unwrapped phase value to a group. Unwrapped 51 phase values with a common group number represent phase values which are members of the same continuous unwrapped phase surface. The phase unwrapping groups are useful when processing data that cannot be unwrapped into one contiguous phase surface as they allow one to process each group independently. 3.4.3 D E M Updating Implementation Issues Each unwrapped phase group is checked to ensure that it has the correct absolute phase ambiguity. The ambiguity is checked by comparing the input D E M mean altitude and the output D E M mean altitude for each group. Mean D E M values which are inconsistent are a sign that the unwrapped phase is in error by at least one half-wavelength. The unwrapped phase can be re-corrected and the baseline estimation procedure can be redone with the input coarse D E M to yield a refined height estimate. 3.4.4 Computational Resource Requirements We implemented the D E M refinement algorithm in Matlab [47]. Our system was constrained by • system memory, • disk space, and • amount of computation. In the following paragraphs, we discuss these resource requirements in the context of the D E M refinement algorithm. The largest memory requirement occurs during the phase unwrapping processing. The weighted least squares phase unwrapping algorithm we implemented 52 Data store input coarse D E M input interferogram residual interferogram second slant range data (r2) filtered, residual interferogram unwrapped phase phase unwrapping mask phase unwrapping groups output refined D E M Bytes per sample 8 8 8 8 8 8 1 4 8 Number of range samples Number of azimuth samples Nr Nr Nr Nr Nr/3 Nr/3 Nr/3 NaNa Na Na Na/3 Na/3 Na/3 Na/3 Na/3 ' iVr/3 Nr/3 Table 3.1: Data stores required for the D E M refinement algorithm. See Section 3.4.4 for discussion of the table values. requires all the wrapped phase data to be in memory. Because of the typically large size of input interferogram, Matlab's choice to represent all numbers internally as doubles, and a system memory limit (including swap space) of 1 Gigabyte; we were required to subsample the residual interferogram before phase unwrapping. However, the subsampling was easily implemented on the residual interferogram because of the filtering facilitated by the D E M flattening. Note that the initial generation and filtering of the interferogram are performed on a small number of contiguous range lines (typically less than 25) out of the whole whole data set so that very large interferograms can conceivably be generated. A large amount of disk space (on the order of 1/2 gigabyte per input interferogram for the data processed in this thesis) was required as different data components were all written to disk for each input D E M during processing. The data stores of significant size are described in Table 3.1. Assuming the original interferometric SAR data size was Nr range samples by Na azimuth samples, and a nominal subsampling ratio of 3 in both range and azimuth was applied, the net amount of disk space required is on the order of 35 Na Nr bytes. The computational complexity (sometimes called the "work") of the D E M refinement algorithm was dominated in our implementation by the computation of 53 the range-azimuth frequencies estimated during the D E M flattening algorithm. This is because the range-azimuth spectrum was calculated using a DFT-like calculation that has 0(n ) computations. The number of operations was limited for this calcu4 lation by subsampling the residual interferogram aggressively in range and azimuth to give a square array of data with approximately 512 samples to a side and limiting the number of sample frequencies investigated in the spectrum. The peak frequency of the residual interferogram spectra was estimated by fitting the output spectra parametrically to determine the peak location. A similar approach using an F F T for coarse frequency estimation followed by a fitting procedure using a small number of D F T coefficients was used for the range and azimuth frequencies. Because of the 0(n log n) operations requirement for the FFTs, these calculations were performed far quicker than the range-azimuth frequency component calculations. There are two other stages in the algorithm that require more than linear amounts of calculations: • Producing the input coarse D E M in slant range coordinates requires an initial triangulation phase which is 0(npts log npts) work where npts is the number of samples in the input coarse D E M . These points were gridded as a function of the triangulation at a grid spacing commensurate with the input D E M data and then oversampled to the slant range grid spacing of the interferogram data. The gridding and oversampling functions require approximately linear number of operations that varies with the required number of output points. • The D E M flattening algorithm and D E M updating algorithm both minimize functions using the Levenburg-Marquardt algorithm (see Section D.l) to determine the optimal baseline parameters. The Levenburg-Marquardt algorithm requires an approximation to the Hessian matrix calculated using the Jacobian at each iteration. The Jacobian matrix has size Ndat by Nvars where Ndat is the number of sample points of the problem and Nvars are the number 54 of parameters. The Hessian matrix approximation requires 0(Nvars 2 Ndat) computations. Because low order models of the baseline parameters generally suffice to characterize the baseline, this computation is dominated by the number of D E M data points investigated. For example, for a 100 square kilometer scene, a linear model suffices for each baseline parameter (i.e. Nvars = 4) and there are approximately 100 • 100 independent GTOPO30 data points that could be used for the minimization problem. In the processing examples, we subsampled the input D E M in range and azimuth to decorrelate errors in the original input coarse D E M data and to decrease the amount of computations required for these stages of the processing. Generally, we attempted to maintain the data at a sample spacing consistent with the input coarse D E M sample spacing. 3.5 A l g o r i t h m P e r f o r m a n c e Issues In the following subsections, we analyze the performance of the two novel stages of the D E M refinement algorithm. In Section 3.5.1 we discuss the factors which impact the flattening procedure. In Section 3.5.2 we discuss the issues affecting the D E M updating procedure. Both algorithms are sensitive to mean offsets, linear trends in range and azimuth, and bi-linear trends in the input coarse D E M as discussed in Section 3.5.3. A performance bound for the output InSAR D E M is derived in Section 3.5.4. 3.5.1 D E M Flattening Algorithm Performance Issues The effect of error in the initial estimate of S(r) is discussed in Section 3.5.1.1. We discuss the residual phase ramp estimate in Section 3.5.1.2. We analyze the spectrum of the residual interferogram in Section 3.5.1.3. 55 3.5.1.1 Errors in the Initial S(r) Estimate The function of the initial estimate of 8(r) used in the D E M flattening algorithm is to provide an accurate estimate of the mean value of S(r) that is to be fitted . 2 If the mean value of 5(r) is specified accurately, the requirement for absolute phase unwrapping in the algorithm is eliminated. In addition, least squares fits preserve the mean value of the target function if the noise errors are well-behaved (see Appendix D.4). Therefore one must assess the accuracy of the mean from the initial estimates of 6{r). One procedure for registering satellite SAR images for interferometry is to estimate the offset between image chips of the two images used to form the interferogram. These offsets are then modeled and a fit is made to determine the resampling between the reference and candidate images. A tradeoff must be made between the accuracy of the estimated offsets and the number of chips that are processed. We would like to calculate the number of chips required so that most of the time, the error in the mean of the calculated S(r) corresponds to less than ir error in the output phase: 2^(r) < ^ , where G$(r) 1 S (3-10) the standard deviation of the <5(r) estimates. Assuming that the errors in the offsets estimates are distributed as independent zero mean Gaussian random variables, the factor of 2 limits one to achieving accurate estimates of the mean of S(r) 95% of the time. One can calculate a numerical approximation for the mean of the average assuming independent samples with constant standard deviation as "M\JN < reg where N reg is the number of image chips examined. From [1] 2 _ *'<••> 2 (3.1D u 1-M — 3p 2 ^AT' The slopes of the fitted function are calculated using the frequency estimates. 56 . ( 3 . - 1 2 ) where H — coherence magnitude, p r — the slant range sampling interval, N s = the number of samples in the image chips. The condition for registration becomes p nN 2 s N reg < 64- [ 6 A 6 ) Rearranging eqn. (3.13) gives the required number of chips on the basis of coherence magnitude (fi) and samples in the image chips as This function is plotted for different chip sizes in Figure (3.5) for ERS standard beam data. In general in satellite InSAR, one is interested in processing interferograms of reasonable coherence (fi > 0.3) which implies that a tradeoff between number of chips and chip sizes must be made to keep N reg to a reasonable number. See Section 3.5.3 for more discussion of the effect of mean S(r) errors. 3.5.1.2 Estimation of Residual Phase Ramp The maximum likelihood (ML) estimates of the residual interferogram phase amount to maximizing the D F T magnitude in the range and azimuth directions independently and then together and and calculating the phase of the sum of the residuals. The asymptotically best performance for the estimation process is defined by the Cramer-Rao (CR) bounds for the process (see Appendix D.3). Referring to the variables of eqn. (3.6),the CR bounds assuming independent samples are > 72(l-/x ) p Nr Na (Na - 1) (Nr - 1)' - n NaNr(Na 2 (w rA - E \(w a W ) r<a 2 2 6(1-M ) 2 w) 2 a 2 2 2 2 57 - 1) (3.15) (3.16) 0.3 0.5 0.4 0.6 0.7 Coherence Magnitude Figure 3.5: Number of chips required to ensure that no absolute phase unwrapping is required versus coherence magnitude and chip sizes of 16, 32, 64, 128, and 256 for ERS standard beam data. 6 (1-M ) 2 (w r [(1>o WrY > H NaNr(Nr' 2 - &y> 2 - 1)' (1 ~ M ) 2iJ?NaNr' 2 (3.17) (3.18) where Nr the number of range samples, Na the number of azimuth samples, EN the expectation operation, the mean value of the estimate. These functions are plotted in Figure (3.6). It is clear that for interferometric applications where one assumes that the coherence magnitude (fx) is greater than about 0.3, one should in theory always be able to estimate the residual phase term 58 accurately. These values are analogous to those found for the single tone in additive white noise case [48, 41]. However, CR bounds are derived using asymptotic arguments. In general, the performance of maximum-likelihood estimates will undergo a sudden sharp decrease in performance from the CR bound at some specific SNR, called the threshold. This effect has been well documented for the single tone frequency estimation problem [48, 49, 50] and could play a role in the residual interferogram phase estimation problem also. James [50] proposed a semi-empirical method to determine the threshold SNR based on an indicator function that is a function of the mean-squared phase error. The onset of threshold can be approximately estimated for the one-dimensional sinusoid as occurring when the phase variance of the ML estimate is 0.0625 rads . 2 We can use this indicator function for the residual phase estimation problem as well since essentially the same process (a DFT) is being used to determine the frequencies. For our use then, the threshold is defined by ^o'J^lr 2/j, NaNr = 0.0625. 2 (3.19) ' v The threshold coherence magnitude can be derived from this relation if Na and Nr are defined. A plot of the threshold value of fi assuming Na = Nr is shown in Figure (3.7). The receiver noise should not be a limiting factor because of the large number of samples in the interferogram which may be used to make the frequency estimate [48, 49, 50, 41]. In addition there is an implicit assumption that the SNR of the interferogram is sufficient to exhibit phase fringes, which implies that some useful frequency estimates may be extracted from the interferogram. 59 Cramer-Rao Bound lor wxy in rada^sample 0.1 0.2 0.3 0.4 0.5 0.6 Coherence Magnitude 0.7 0.8 0.9 (a) C R Bound for w ,a r Cramer-Roo Bound lor wx, wy in rad^/aample 0.1 0.2 0.3 0.4 0.5 0.6 Coherence Magnitude 0.7 0.8 0.9 (b) C R Bound for w and w a Cramer-Rao Bound forty in rads 0 1 02 03 04 OS 06 Coherence Magnitude a 2 07 08 09 (c) C R Bound for ip 0 Figure 3.6: Cramer-Rao bounds calculated for Nr = Na for the model paramet of the residual interferogram phase. 60 61 3.5.1.3 Spectrum of the Residual Interferogram A critical step in the flattening algorithm is the estimation of the residual phase ramps from the maximum spectral peaks of the residual interferogram. In the following we discuss the possibility of error in this process. The simplified model for the interferogram used to develop the D E M flattening algorithm in Section 3.2 neglected one important component of InSAR data: the noise signal. Noise in the interferogram is generated by system noise in the two SAR images and by the uncorrelated portions of the SAR signal data spectrum (see Appendix A). Assuming that the noise level in both images is constant and that no weighting is applied during SAR processing, the spectrum of the resulting noise in the interferogram will be a triangle function in range and azimuth because of the multiplication of the oversampled SAR images used to form the interferogram. When forming the residual interferogram, the constant (flat-earth) part of the geometrical phase model tends to shift the peak of the noise spectrum away from zero in the same manner as the interferogram signal data is shifted towards zero. Normally, one minimizes the contribution of interferogram noise to errors by halfband filtering the data (see Appendix B). Atmospheric disturbances are another source of noise in the interferogram phase. These disturbances tend to be low-frequency, large area disturbances. They are more likely to distort the peak location but not the existence of a peak. Errors of this sort are discussed in Section 3.5.3. One can view the spectrum of the residual interferogram as the convolution of the spectra of the geometry error phase term (approximately a sinusoid) with the D E M error term. A necessary condition for the residual interferogram spectral peak to correspond to the geometric phase ramp is that the spectrum of the D E M error term / e * J dem have a significant maximum at zero frequency. Statistical character- ization of the spectrum of / e * J d e m is difficult because of the unknown topography 62 and the impact of interpolation of noisy D E M data from the D E M sample positions ( e.g. a U T M grid ) to the radar data grid. However we can make some observations about the impact of this noise source on the algorithm's performance. The deviation of ^dem from zero depends on the normal baseline of the data in concert with the distribution of D E M errors (Ah) and the ground range (x) (see eqn. (3.4)). In general, one expects that the mean of the D E M error will be at or close to 0 so the distribution of ^dem will have a local peak at zero frequency. If this peak at zero frequency is the only peak in the spectrum due to D E M error, then the spectrum of the residual interferogram will be directly related to the geometry error. If there is a larger peak at non-zero frequency in the spectrum of e * J dem , the phase ramp estimated as a correction to the second slant range will be in error and the flattening algorithm's results will be unreliable. For a given set of D E M errors, the normal baseline essentially determines how broad the spectrum of the D E M error signal becomes. As the normal baseline becomes larger and hence the altitude of ambiguity becomes smaller, the energy in the spectrum of the residual interferogram spreads and the possibility of having multiple peaks near zero frequency increases. We illustrate this effect diagrammatically in Figure (3.8). For the D E M errors shown in Figure (3.8) (a), the spectrum of the complex valued D E M error phase term was plotted for altitudes of ambiguity of 500 m, 250 m, and 50 m respectively as shown in Figure (3.8) (b). For the 500 m case, the D E M phase error spectrum has a sharp peak at 0 frequency. A geometry error could therefore be recognized easily from the residual spectrum of the interferogram during the unwrapping processing with this data. When the altitude of ambiguity is 250 m, the magnitude of the peak at 0 frequency is less prominent and one of the sidelobes has become significantly larger. This is likely a marginal case for recognition of the geometrical error frequency component. At an altitude of ambiguity of 50 m, 63 the spectrum of the residual interferogram is very broad with no prominent peaks and no chance of recognizing a particular frequency component as the frequency associated with the geometry error. The spectrum of the residual interferogram indicates how well the algorithm has performed. If there are a large number of discrete peaks near zero frequency, then it is likely that an error exists in the geometrical phase, and the baseline parameters will also likely be in error. The baseline parameters can be re-estimated after phase unwrapping to ensure accurate height estimates are generated. However, if the spectrum contains only a single peak at zero frequency, then we know that the input D E M matches the structural information in the interferogram to sufficient detail so that the geometrical phase errors can be correctly recognized. The baseline parameters generated from such a process will likely be reasonably accurate subject to trends and biases in the input D E M (see Section 3.5.3). Note that in general, it is always useful to refine the baseline estimates by fitting the InSAR topography to the input coarse D E M (see Section 4.1.2). 64 -0.015 -0.01 -0.005 0 frequency 0.005 0.01 (b) Spectrum of D E M error term. Figure 3.8: Sample spectra for D E M errors. 65 0.015 3.5.2 D E M Updating Performance Issues In the final stages of processing, the baseline parameters may have to be re-estimated using the input coarse D E M as a template function. The errors that can specifically affect this calculation are • Random phase noise. Random phase noise lies between — ir and TT. It contributes to the variance of the baseline parameter estimates and therefore also the output terrain height estimates. Usually there is an abundance of InSAR data points so that averaging may be taken advantage of to ensure good quality baseline parameter estimates and hence height estimates are achieved. • Phase unwrapping errors. Phase unwrapping errors can either be local (an isolated pixel or small number of pixels) or global in the sense of a group of pixels unwrapped to the wrong ambiguity. Unwrapping errors tend to occur at multiples of the ambiguity height of the interferometer and tend to create long "tails" in the error distribution of the topographic height estimates. Both local and global unwrapping errors tend to bias the estimate of the baseline parameters with the global errors causing severe errors in the output height estimates. By checking that the average height of contiguously unwrapped groups of pixels is consistent with the input coarse D E M mean, one can eliminate the unwrapping errors for large groups. • Random input D E M errors. Input D E M errors cause similar problems to the random phase noise. Again, we can take advantage of the large number of data points available to us to minimize the impact of D E M errors on the output topographic heights. 66 One other source of noise is trend and offset errors in the D E M . These errors are discussed in the next section as they affect both the D E M flattening and D E M updating algorithms. 3.5.3 Effects of D E M Trend Errors Considering the bi-linear nature of an error in h (input D E M height) or r2 — r (slant range difference) due to baseline parameter errors, the D E M refinement algorithm is sensitive to any error sources which have a similar linear character. Possible sources of this type of error include: • a, the satellite altitude, • interferogram phase, • mean value of S(r) estimated from the SAR image registration parameters, • coarse D E M model. The satellite altitude is used as a parameter for the geometry calculation. Errors in a result in bias and some slight error trends in the output D E M . Errors in the interferogram phase are caused by atmospheric artifacts [51, 15] for example. Linear errors in the interferogram phase are wrongly identified as geometry errors by both baseline estimation algorithms, leading to maladjustment of the estimated interferometer geometry. Finally, linear trends and offsets in the D E M used for the geometry calculation also lead to errors in the estimated geometry. In general, for the algorithm to arrive at the correct value of B and 0, there should be no linear or constant trends either in the D E M or <5(r) derived from the interferogram or from an error in the satellite altitude. However, if the eventual goal of the baseline parameter estimation exercise is to estimate topography from the unwrapped interferogram phase, the error in 67 the baseline estimate can be beneficial. In particular, trends in the interferogram phase, error in the estimate of the mean value of S(r), and altitude errors can all be compensated by estimating the "incorrect" baseline geometry. Because the height error terms due to baseline parameter errors are linear functions they "absorb" systematic errors such mean offsets and error trends in the satellite altitude and S(r). Note that trends in the input D E M used to calculate 6(r) are not compensated by improper baseline parameter estimates in the height calculation. The trends and mean values (see Appendix D.4) of the input D E M are conserved in any D E M calculated from the unwrapped phase. This is important because it means that the proposed method of D E M refinement is not suited to removing trend and mean errors in existing DEMs. However, fine details of the topography which are present in the interferogram may be used to update the existing biased coarse D E M even though it may have linear bias terms. These effects are demonstrated in Section 4.2. 3.5.4 E x p e c t e d Baseline and Height Accuracy- Prediction of the accuracy of the baseline parameter estimates and height estimates for the D E M refinement algorithm is hampered by the nonlinear nature of the minimization problem, and the interpolation of the D E M values to the slant range radar grid, and in the case of the flattening algorithm, the dependence of algorithm performance on estimation of peak frequencies in the residual interferogram. Following the discussion of Section 3.5.1.3, we do not expect the baseline parameters derived during the D E M flattening algorithm to be accurate because of the likelihood of extraneous peaks in the cases where the normal baseline is suitable for topographic estimation (i.e. the altitude of ambiguity is small). For the D E M updating algorithm, the net effect of errors in B and 0 is to create an approximate bilinear trend and possibly an offset in the InSAR output D E M . For the output InSAR D E M to match the input coarse D E M in the mean- 68 squared error sense, the output InSAR D E M must have the same mean and slope characteristics as the input D E M . Therefore, analyzing the mean and slope errors of a bilinear function fit to the input D E M errors will give some indication of the best possible accuracy of the output InSAR D E M estimates. It can be shown under the assumption of zero mean, Gaussian distributed, uncorrelated errors in the N by N element input coarse D E M (see Appendix D.5), that the standard deviation of the mean error of the output plane fit is _ °~DEM Pmean — — > (o 9 f 1 N the standard deviation of the output plane fit is 0~out = , (6.2L) the standard deviation of the maximum error across the swath is Ah th swa = ^ JV ( i v V l ) ^ ° D E M - (3.22) These functions establish lower bounds on predicted output D E M accuracies as a function of input D E M errors. They are tabulated for GTOPO30, D T E D , and TRIM DEMs in Table 3.2 for N = N = 100. From the table, it is clear that Inr a SAR DEMs derived from reasonable quality interferograms using D T E D and TRIM DEMs can potentially be very accurate. DTED-1 and DTED-2 quality DEMs can potentially be derived from a full 100 km by 100 km input GTOPO30 D E M . Further improvements in accuracy for GTOPO30 case could be achieved by processing larger amounts of data. Following Section 2.4, the error bound defined above can be used to define a range of baseline magnitude and orientation errors which yield InSAR DEMs consistent with the error properties of the fitted plane. Due to errors from additional 69 input Source D E M GTOPO30 DTED TRIM &DEM m 86 18 3.0 output std. of mean in m 0.9 0.2 0.03 output std. in m 2.6 0.5 0.09 output Ah f}i swa 5.9 1.2 0.2 Table 3.2: Estimated mean and across swath height errors in planar fit to D E M errors for N = N = 100. For the GTOPO30 D E M this corresponds to a standard ERS scene. For the D T E D and TRIM DEMs this corresponds to ground coverage of about 10 kms. r a sources such as independent phase noise in the unwrapped phase, phase unwrapping errors, and the assumption of uniform independent error statistics, the InSAR D E M error must be greater than or equal to the possible error of the plane fit to the D E M errors: Ah (InSAR) > Ah swath swath (DEM). (3.23) The mean height error in the output D E M is difficult to derive analytically but one can use the proxy estimate of the average of the height errors across the swath: Ak I«SAR) mm{ = ^ ( A B . A e ) + A/w(Afi.A9) ( J J J 4 ) to set a lower bound for the mean InSAR D E M height error which must be consistent with the standard deviation of the mean of the input D E M : 6h (InSAR) > mean ODEMmean- (3.25) Using the above equations and Appendix D.2, the range of possible values of baseline parameter errors can be found which produce InSAR D E M errors within the same range as the possible errors in the input D E M . An example of a plot of the range of valid errors for the GTOPO30 parameters for a standard ERS scene are shown Figure (3.9). Finally, one can write an expression for the output height accuracy of the InSAR D E M including parameter errors and random independent phase noise but 70 -0.08 h -0.06 -0.04 -0.02 h 0.02 h 0.04 0.06 0.08 b A@ in radians x 10 Figure 3.9: Range of baseline parameter errors for GTOPO30 input D E M for a standard ERS scene with B fa B fa 250 meters. The mean and swath errors are shown in Table 3.2. L excluding the effect of phase unwrapping errors as 2 _ 2 °~h — pars a , 2 + Otf, (3.26) where: a\ J 2 <7-„„ pars c _2 0$ = the output height standard deviation, = the standard deviation due to the cumulative effect of parameter errors, the standard deviation of the phase noise. All of the parameters for the InSAR topography problem (a, mean S(r), B, and 0 ) combine to give biases and linear trend errors in the output D E M . Using our approach, these errors can be reduced by adjusting the values of B and 0 to provide a minimum least squares fit to the input coarse D E M . The second independent error 71 term depends on the amount of phase noise which can be considered to be an additive independent error term. Note that phase unwrapping offset errors affecting a single large unwrapped phase group will cause significant errors in the baseline parameter estimates and hence the output height estimates. 3.6 Summary We have presented and analyzed an algorithm for estimating topography using coarse DEMs. Key to the algorithm's formulation is estimating the baseline both before the interferogram phase is unwrapped and after the interferogram phase is unwrapped. Estimating the baseline before the phase is unwrapped allows us to flatten the interferogram by removing the topographic model of the interferogram phase. Estimating the baseline parameters after the phase is unwrapped allows the baseline parameters to be refined for topographic height estimation. Lower bounds on the algorithm performance were derived using a hypothetical experiment involving plane fits to random noise. Flattening with an existing D E M facilitates phase unwrapping and can provide an accurate estimate of the interferometer geometry. One can check whether the algorithm has succeeded or not by examining the spectra of the residual interferogram. Multiple peaks near zero frequency suggest that the baseline parameter estimates will be in error. The performance of the D E M flattening algorithm is determined by the quality of the input D E M relative to the variability of the terrain and the sensitivity of InSAR geometry to height errors. Large height errors coupled with large normal baseline will not likely give reliable baseline parameter estimates. Even if the baseline parameter estimates are poor, the resulting residual interferogram can be post-processed (filtered and downsampled) to reduce the possibilities of problems in phase unwrapping processing. 72 In the D E M updating algorithm, baseline parameters can be mis-estimated to compensate for errors in other parameters such as satellite altitude errors and absolute phase unwrapping errors. The baseline estimation algorithm is sensitive to linear trends and biases in the input coarse DEM. Processing as large an area as possible will minimize the likelihood of the coarse D E M having spurious bias and trend errors. It is tempting to suggest an iterative approach to processing whereby the output InSAR height estimates are used as the input "coarse" D E M for another application of the D E M refinement algorithm. However, since the updating algorithm conserves trends in the input DEMs, application of the algorithm a second time will not lead to significant improvement in the output InSAR heights. In addition, the noise in the output InSAR D E M will be correlated with the phase noise in the interferogram so the improvement through filtering the residual interferogram will be minimal. 73 Chapter 4 Simulations of D E M Refinement Algorithm Simulated datasets were generated and processed to investigate the D E M refinement algorithm's performance. Two different scenarios were investigated: random D E M errors in Section 4.1 and low-frequency trend errors in Section 4.2. Summary remarks are given in Section 4.3. Experimental results for an ERS interferometric pair and a RADARSAT interferometric pair are presented in Chapter 5. 4.1 Random D E M Error Simulations For the random D E M simulations, we synthesized a noise free interferogram from an approximately 10 by 40 km sample of TRIM [52] data from Canadian map sheet 92O (Taseko Lakes). Using a noise free interferogram allows us to isolate the efficiency of the flattening and updating algorithms in estimating the baseline parameters. Nominal ERS SAR parameters were used to generate an interferogram with 2000 range by 350 azimuth samples for a succession of normal baseline values. 74 The TRIM data lies on an irregular grid with nominal sample spacing of 100 m in relatively flat terrain and 75 m in areas with significant slope [52]. The synthetic interferogram phase was obtained by interpolating the scattered TRIM height points to a slant-range/azimuth SAR grid. To test the effects of different D E M accuracies, lower resolution and quality DEMs were created by deleting selected points from the TRIM data [53]. When the pruned DEMs are interpolated to the slant range grid of the interferogram, random (but locally correlated) errors result from the missing data points because of the undersampling in the pruned D E M . Reduced-posting DEMs were generated with nominal 2a interpolated height errors of 15, 25, 45 and 80 meters. The error statistics of the D E M after interpolating to slant range are shown in Table 4.1. GTOPO30 data points [21] were also interpolated to generate a coarse D E M for processing. The GTOPO30 slant range interpolated D E M had errors roughly three times that of the 80 m simulated D E M . In addition, Table 4.1 also includes a prediction of the minimum output InSAR D E M standard deviation in the accuracy column using the analysis of Section 3.5.4. A graphical example of the D E M pruning algorithm is shown in Figure (4.1) for the 25 m and 80 m DEMs. The top row of sub-figures in Figure (4.1) applies to the 25 m dataset while the bottom row of sub-figures applies to the 80 m data set. Each row shows the input D E M sample locations, the interpolated slant range D E M , and the interpolated DEMs' errors with respect to a dataset generated using all the available TRIM samples. The 80 m dataset was generated from approximately 0.5% of the TRIM DEMs' samples while the 25 m data set required approximately 5% of the samples. The difference in the number of input data samples is reflected visually in the lack of detail in the 80 m interpolated dataset with respect to the 25 m dataset and quantitatively in the increased RMS height error of the 80 m data set compared with the 25 m data set (40 m versus 12 m, see Table 4.1). From the analysis of Section 3.5 one expects that the error performance of 75 the algorithm depends both on the baseline magnitude and orientation. We simulated data with a wide variety of baseline magnitudes and orientations and found roughly the same performance in terms of number of iterations to convergence and topographic accuracy of the final estimates in all cases. To illustrate the dominant effect of the perpendicular component of baseline magnitude, we chose one representative case of orientation with mean (f> « TT/4, which tended to have the largest number of iterations to convergence. The simulations were performed with mean baseline magnitudes of 70, 170, 270, 370 and 470 meters at this orientation. For a more realistic simulation, the baseline magnitude was allowed to vary by 10 cm and the baseline orientation was allowed to vary by 100 /xrads in the azimuth direction. Lower bounds on the accuracy predictions for the output InSAR DEMs are given in Table 4.1 also. DEM 15 m 25 m 45 m 80 m GTOPO30 Mean 0.1 0.8 1.2 7.6 2.9 Std. 7.6 12.0 22.5 40.1 116.5 90% 8.2 16.0 32.0 65.0 194.9 Max 246.3 231.9 235.3 367.2 847.5 Min -305.1 -344.1 -405.3 -409.7 -375.4 Output mean 0.1 0.2 0.3 0.5 1.5 Output std. 0.2 0.3 0.5 0.9 2.7 Table 4.1: Statistics and lower bound on output D E M accuracy for degraded D E M height errors in meters after interpolating the D E M samples to the simulated slant range grid. 76 400 0 5 10 Along Track in km Input D E M samples. Interpolated D E M . a) 25 m simulated dataset. D E M errors. • V*. •i -—i w • 40 1000 35 t " ,•.'•* 1 30 | • .'!'::V\s 1 25 fr 1 ?20 o 0 15 K*& • • . 10 . •*« .." .* 5 0 0 5 10 Along Track in km Input D E M samples. Interpolated D E M . b) 80 m simulated dataset. D E M errors. Figure 4.1: Example of interpolated DEMs in slant range for the 25 m and 80 m data set. See page 75 for more information. 77 4.1.1 D E M Flattening Simulations The results of the D E M flattening algorithm are summarized in Table 4.2. The following factors are reported in the table for each combination of input coarse D E M and baseline magnitude: • mean error of the baseline magnitude, B; • mean error of the baseline orientation, 0; • the state of spectrum of the residual interferogram, either single peaks (y/) in all the residual interferogram spectra or multiple peaks (x) in one or more of the residual interferogram spectra; and • the error of the height estimate using the algorithm baseline parameters and ideal unwrapped interferogram phase. The D E M height errors give concrete values for the D E M error caused by the given baseline parameter error and calculated with respect to the accurate D E M used to create the interferogram. From Table 4.2, we see that the algorithm works well with the higher quality input DEMs (the 15 m and 25 m data sets), where the refined baseline estimates are accurate to within a few millimeters. The flattening algorithm produces baseline parameter estimates yielding relatively accurate DEMs for the 15 m, 25 m, and 45 m case. These output InSAR DEMs have accuracy statistics on the same order as the predictions made in Table 4.1 but in general the errors in output InSAR DEMs are larger than predicted using the input D E M error statistics. For the 80 m simulated D E M (when the input D E M has 40 m standard deviation errors), the baseline parameters are not estimated correctly when the baseline magnitude becomes larger than 70 meters. As the baseline magnitude and hence the normal baseline increases, the residual interferogram's phase errors due to D E M errors become large enough 78 to make an accurate estimate of the required frequencies due to geometry errors unlikely. The far coarser and higher error GTOPO30 data set did not provide the correct baseline estimates for any baseline length. In further experiments we found that the normal baseline had to be less than 20 m for the algorithm to work with the GTOPO30 data. This is likely due to one large outlier in the original D E M which distorts the geometrically estimated S(r) significantly from the true value Baseline in m 70 170 270 370 470 Dataset 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 B Error in m 0.001 -0.003 0.009 0.009 -2.171 0.002 -0.006 0.010 0.508 -0.653 0.003 -0.008 -0.007 0.418 -0.856 0.004 -0.007 0.032 -0.260 -1.318 0.005 -0.002 0.101 -0.642 -2.382 Q Error in mrads 0.009 -0.036 0.130 0.144 -32.776 0.010 -0.032 0.060 2.932 -3.932 0.010 -0.026 -0.024 1.493 -3.155 0.011 -0.016 0.090 -0.702 -3.554 0.012 -0.002 0.214 -1.307 -5.056 Spectrum check V V V V X V V V X X V V V X X V V V X X V V X X Mean h in m 0.2 0.9 1.3 6.0 -62.0 0.3 1.0 1.3 2.6 -10.9 0.3 1.0 1.1 -18.0 12.9 0.3 0.9 1.1 0.9 12.0 0.2 0.7 -0.1 21.8 19.0 Std. h in m 0.2 0.5 1.6 2.4 299.7 0.2 0.4 1.1 66.3 81.2 0.2 0.3 0.8 28.3 38.4 0.2 0.3 1.5 9.0 41.9 0.2 0.3 2.6 15.4 47.9 Table 4.2: D E M flattening simulation results. The spectral check column refers to the condition of the residual interferogram spectrum. Residual interferograms having strong single peaks at 0 frequency are denoted by sj while interferograms with multiple peaks near zero frequency are denoted by X . 79 (see Figure (5.4)). From our analysis in Section 3.5 and more particularly Section 3.5.1.3, we expect that the spectrum of the residual interferogram in the range, azimuth, and range/azimuth directions will be useful quality measures for assessing the accuracy of the baseline parameter estimates. To demonstrate the utility of examining the spectra of the residual interferogram, we examine the 470 m baseline case in Figure (4.2). The left hand sub-figures in Figure (4.2) show the spectra of the residual interferograms generated from the flattening algorithm's estimate of the baseline parameters and the associated input D E M . The right hand sub-figures show the spectra of the residual interferogram generated with the true baseline parameters and the same DEMs as used in the flattening algorithm. Each one of the sub-figures shows the spectrum in the sub-figure's labeled direction (i.e. range, azimuth or range/azimuth) for all input DEMs and the ground truth (TRIM) D E M . The spectral estimates were made for each input D E M in the desired direction independently and then the spectra were plotted with the same vertical scale in a stack-like configuration to facilitate comparison of the different results. In each sub-figure, input D E M quality decreases (i.e. D E M errors increase) as one goes from the bottom plot to the top plot. The input D E M associated with each plot within each sub-figure is shown as a label to the left of the plots. Firstly considering the spectra generated using the true baseline parameters (i.e. the right hand plots of Figure (4.2)), we note that as the D E M quality decreases the spectra of these residual interferogram tends to "spread" across a larger range of frequencies. In the "0 m" data set (the lowest rows of each plot), the dominance of the DC term is apparent. At the correct baseline and with no input D E M errors, the "flattened" interferogram is a constant real value and the spectrum are therefore sine functions. As D E M errors are progressively introduced, the dominance of the DC term vanishes, and significant distortion enters the spectra of the residual interferogram. As the energy in the residual interferograms spectra spreads due to 80 the decrease in D E M quality, an increasing number of local peaks become apparent near zero frequency. In the development of the D E M flattening algorithm (see Section 3.2), it was assumed that the D E M errors contributed only one significant peak to the residual interferogram spectra that was convolved with the nearly linear baseline parameter phase error terms. As we see in the right-hand plots of Figure (4.2), this assumption is not valid when the D E M errors start to become large. Predictably, this results in poor algorithm performance for large D E M errors as shown in Table 4.2. Transferring our attention to the left-hand side plots of Figure (4.2) with the results of Table 4.2 in mind, the D E M flattening's algorithm goal of driving a significant spectral peak to zero frequency is realized for most input D E M data sets except for the GTOPO30 data set where the iterations were in general stopped after reaching a maximum number of iterations. It is evident that the flattening algorithm performs well when the spectra of the residual interferogram is dominated by its DC term but fails to yield correct baseline parameters when the spectra shows significant peaks at non-zero frequencies. The extra peaks in the residual interferogram spectra correspond to the input D E M height errors as shown in the right-hand plots of Figure (4.2) and are not introduced by the D E M flattening algorithm. These extra peaks confound the D E M flattening algorithm, which uses the frequency associated with the largest peak to calculate the refinement to 5(r) in the fitting algorithm. 81 (a) Azimuth: D E M flattening results. (b) Azimuth: true baseline parameters. (c) Range: D E M flattening results. (d) Range: true baseline parameters. (e) Range/Azimuth: D E M flattening results. (f) Range/Azimuth: true baseline parameters. Figure 4.2: Comparison of residual interferogram spectra plotted in dB for azimuth,range, and azimuth range frequency components. See page 80. 82 4.1.2 D E M Updating Simulations The D E M updating simulations processed the output of the D E M flattening algorithm to generate refined D E M estimates. The D E M updating simulations encompass both both the phase unwrapping processing and re-estimation of the baseline parameters using the procedure of Section 3.3. The phase unwrapping processing was performed in two stages. In the first stage, weighted-least squares phase unwrapping was performed using binary weighting. An initial weight mask was computed for all data sets based on the small number of zones of radar shadow throughout the image. This initial weighting was augmented by a mask that was set to zero when the magnitude of the range phase derivative of the 470 meter dataset exceeded n rads. The second mask was applied as an approximation for zones of low coherence magnitude associated with terrain effects. The 470 meter data set was used to calculate the second mask as it is the worst case for all the simulations. The same mask was used for all data sets to get a uniform basis of comparison. For each data set, this initial weighting was then augmented by the residues produced by the current residual interferogram phase. After the initial weighted least squares unwrapped phase estimate, the phase residual was unwrapped using Flynn's Minimum Discontinuity method [54, 55]. The unwrapped phase was then used to update the D E M flattening algorithm's estimate of <5(r) for estimation of the final baseline parameters and refined D E M estimate. Results of the experiment are shown in Table 4.3, Table 4.4, and Table 4.5. Table 4.3 reports the error statistics of the input coarse D E M pixels corresponding to the output refined D E M pixels. This ensures that we compare input and output D E M quality over the same pixels. These values are different from the values reported in Table 4.1 because the updating simulations only produces results for a subset (about 90 %) of the input coarse D E M values. The binary mask tends to discriminate against pixels in steeply varying terrain where height errors are likely 83 larger. The D E M error statistics are consistent with the input D E M error statistics in Table 4.1 but with smaller standard deviation in general. Table 4.4 reports the error statistics of the refined InSAR D E M along with the errors of the baseline parameters, the number of phase unwrapping errors, and improvement in the D E M quality with respect to the input coarse D E M . The improvement in the D E M quality is expressed as the ratio of the input coarse D E M standard deviation to the output InSAR D E M standard deviation. The performance statistic is greater than 1.0 for all input DEMs except for the GTOPO30 D E M which generally demonstrates the improvement due to phase unwrapping and re-estimation of the baseline parameters. Note that when checking for offset phase errors of unwrapped phase groups, the GTOPO30 dataset often failed because the mean topographic height calculation relied on too few input D E M points to establish an accurate estimate of the mean height of the group. The percentage of phase unwrapping errors is generally correlated with the decrease in output D E M quality. Unwrapping errors bias the baseline parameter estimates as unwrapping errors produce topographic heights which are centered at height error values corresponding to the 2n ambiguity height of the interferogram. The process of least squares fitting generally tries to drive the mean error of the fitted parameter towards 0 (see Appendix D.4). Phase unwrapping errors therefore tend to bias and inject error trends into the output height estimates. The mean values of the output InSAR D E M are approximately the same as the input mean D E M values shown in Table 4.3. This demonstrates the D E M updating algorithm's tendency to conserve the mean of the input data. In general, the baseline errors from the D E M updating procedure are either approximately equal or larger than the same errors for the D E M flattening procedure as the baseline parameters adjust to compensate for the mean value of the input coarse D E M . To assess whether or not the D E M updating procedure increased the output 84 InSAR D E M product quality, we calculated the error statistics of the D E M generated using the D E M flattening baseline parameters and the output unwrapped phase and compared them with the height errors of the input DEMs also. The statistics are reported in Table 4.5. Again, the ratio of input to output D E M errors using the D E M flattening estimate of the baseline parameters are reported using the unwrapped phase estimates. The performance ratio is smaller than the same ratio for the D E M updating algorithm in general although the baseline errors are in general equal or smaller than the baseline parameter errors from the D E M updating algorithm. Note also the effect of wrong baseline parameters because of multiple peaks in the residual interferogram spectra that generates output DEMs with large output errors. The improvement in output InSAR D E M quality with the D E M updating baseline parameters is unsurprising because the updating algorithm minimizes the output D E M variance with respect to the input coarse D E M . In contrast, the D E M flattening algorithm minimizes a fit of an estimate of the slant range difference between the two images. The two problems are not identical but are close enough to yield similar results. 85 Baseline in m 70 170 270 370 470 Dataset 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 Mean h in m 0.3 1.2 1.5 8.2 -2.8 0.3 1.3 1.6 8.3 -4.2 0.3 1.3 1.6 8.5 -2.8 0.3 1.3 1.6 8.8 -2.7 0.3 1.4 1.7 9.3 -2.5 Std. h in m 3.8 8.5 17.5 36.1 108.3 3.7 8.4 17.5 35.9 105.7 3.7 8.3 17.4 35.9 104.3 3.7 8.3 17.4 35.7 101.3 3.6 8.3 17.5 35.4 98.1 Table 4.3: Input D E M error statistics for D E M updating experiments. 86 Baseline in m 70 170 270 370 470 Dataset 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 B Error in m 0.001 -0.003 0.011 0.017 -0.265 0.002 -0.007 0.027 0.041 -2.247 0.003 -0.010 0.044 -0.347 -2.090 0.004 0.186 -0.212 -0.888 -1.198 0.005 0.166 0.237 -0.515 -3.260 0 Error in mrads 0.010 -0.035 0.161 0.262 -4.439 0.010 -0.035 0.165 0.267 -13.492 0.011 -0.034 0.167 -1.305 -7.646 0.012 0.518 -0.584 -2.338 -3.035 0.012 0.364 0.522 -1.019 -6.848 Mean h in m 0.3 1.3 1.5 8.1 -2.9 0.3 1.3 1.6 8.3 -3.9 0.3 1.3 1.6 8.5 -2.9 0.3 1.3 1.6 8.8 -2.8 0.3 1.4 1.7 9.3 -2.4 Std. h in m 0.3 0.5 1.8 4.3 47.8 0.4 0.5 1.9 5.0 112.6 0.4 0.6 6.3 12.6 80.5 0.5 7.7 9.9 18.7 68.2 0.9 5.8 7.9 14.3 67.2 % o(in) Unwrapping Errors a(out) 12.3 0 15.7 0 9.9 0 8.4 0 2.3 49 0 8.7 17.2 0 9.2 0 0 7.3 0.9 45 9.2 0 14.3 0 1 2.8 11 2.8 1.3 79 7.4 0 1.1 4 5 1.8 1.9 40 1.5 75 4.3 0 1.4 4 4 2.2 34 2.5 1.5 46 Table 4.4: Output D E M error statistics for D E M updating experiments. 87 Baseline in m 70 170 270 370 470 Dataset 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 15 m 25 m 45 m 80 m GTOPO30 Mean h in m 0.2 0.9 1.2 5.8 161.8 0.3 1.0 1.1 2.8 -21.3 0.2 0.9 1.8 -4.1 -58.0 0.2 -2.4 5.5 -16.5 -66.6 0.2 -1.5 -3.7 3.4 -39.2 Std. h in m 1.3 1.6 3.0 5.5 302.2 1.6 1.9 3.8 67.3 136.9 1.9 2.2 7.4 34.3 99.6 1.6 10.7 13.3 29.9 101.5 1.7 8.1 10.9 22.0 91.3 a(in) a (out) 2.9 5.3 5.9 6.6 0.4 2.3 4.5 4.6 0.5 0.8 2.0 3.7 2.3 1.0 1.0 2.4 0.8 1.3 1.2 1.0 2.1 1.0 1.6 1.6 1.1 Table 4.5: Output D E M error statistics using unwrapped phase with D E M flattening baseline parameters. 88 4.2 Effect o f T r e n d Errors As noted in Section 3.5.3, errors in certain parameters result in low frequency error trends in 8(r) which may cause errors in the baseline estimate. However, in some cases the algorithm can converge upon an incorrect baseline which counterbalances errors in the parameters, resulting in reasonable quality InSAR terrain height estimates. To demonstrate this feature of the algorithm, we consider four types of errors: • a satellite altitude error of 100 meters, as might be experienced when ERS predicted orbit data are used for estimating the satellite position, • (5(r) errors consisting of a linear trend of two phase cycles in the interferogram (modeling an atmospheric artifact) and an offset of 2 A (modeling an error in estimating the constant offset between the two SAR image slant ranges), • linear errors on the input D E M consisting of an offset of 15 meters with a 5 meter slope in range and a 2 meter slope in azimuth consistent with errors observed between two different DEMs we have processed, • zero mean quadratic errors on the input D E M consisting of quadratic functions in range and azimuth directions with maximum amplitude of 5 meters corresponding to a topographic error with standard deviation of « 2 meters. The baseline estimation algorithm was applied to a noise-free synthetic interferogram with baseline of B = 270 m and 0 = 1.15390 f» 37r/8. The results from the D E M flattening algorithm and the D E M updating algorithm are reported below. Terrain heights were generated with the true unwrapped phase using the baseline parameters generated by the D E M flattening and D E M updating algorithms. The results are summarized in Table 4.6 and Table 4.7. Height errors are 89 reported for both the algorithms' baseline parameter estimates (denoted "Algorithm Baseline" in the table) and the true baseline parameters (denote "True Baseline" in the table) for both algorithms. The D E M error statistics are the same for both algorithms to within the number of significant figures reported in the table. On comparing the height error statistics of the DEMs derived using the algorithm parameter estimates and the true baseline parameters, we see that the baseline estimates for the altitude and slant range errors, despite being in error, can eliminate some of these errors in the terrain height estimate. For example, the altitude error and <5(r) error simulations show an improvement in the quality of the output D E M by using the algorithm estimate of the baseline parameters instead of the true baseline. In these cases the flattening and updating algorithms' baseline parameter estimates are in error, but the baseline error compensates for the error of the other parameters yielding an improved quality output D E M . The quadratic D E M error does not effect the baseline estimate of the algorithm - the correct baseline estimate is produced by the flattening algorithm and almost exactly by updating algorithm. Note, however, that the input D E M error standard deviation was approximately 2 m while the output D E M error using the unwrapped phase from the residual interferogram have standard deviation of about 0.2 meters — a ten-fold improvement. The lack of improvement in the linear D E M error cases demonstrates the both theflatteningand updating algorithm's sensitivity to biases and linear trends in the D E M errors. 90 Estimated Geometry Type of error B Height Errors Algorithm Baseline e mean stdev True Baseline mean stdev Act = 100 meters 269.88 1.15317 0.0 0.2 99.9 0.0 f5(r)errors 270.73 1.15776 0.0 1.2 -410.5 8.2 Linear D E M errors 269.97 1.15384 15.0 1.6 0.0 0.0 Quadratic D E M errors 270.00 1.15390 0.0 0.2 0.0 0.2 Table 4.6: Baseline estimates and height errors in the trend error experiments. Estimated Geometry Type of error B e Height Errors Algorithm Baseline mean stdev True Baseline mean stdev Aa = 100 meters 269.88 1.15317 0.0 0.2 99.9 0.0 <5(r)errors 270.73 1.15776 0.0 1.2 -410.5 8.2 Linear D E M errors 269.97 1.15384 15.0 1.5 0.0 0.0 Quadratic D E M errors 270.00 1.15389 -0.0 0.2 0.0 0.2 Table 4.7: Baseline estimates and height errors in the trend error experiments. 91 4.3 Summary We have simulated an algorithm for refining DEMs using SAR interferometry. Substantial improvement in terms of D E M errors was made using the baseline parameters estimated during the D E M flattening and D E M updating portion of the algorithm. The D E M quality improvement was always greater when the baseline parameters were updated by fitting the InSAR D E M to the input D E M (subject to lack of phase unwrapping errors). This procedure should always be applied in practice to ensure the best possible output D E M quality. The D E M flattening procedure facilitates phase unwrapping by reducing the phase variability and can provide a reasonably accurate estimate of the interferometer geometry. One can check whether the algorithm has succeeded or not by examining the spectra of the residual interferogram. Multiple spectral peaks near zero frequency in the D E M flattened interferogram suggest that the baseline parameter estimates will be in error. Note also that for variable topography and low-quality DEMs, the D E M flattening algorithm will likely not produce accurate DEMs. In addition, the coarse D E M can be used to equalize the phase offsets between different unwrapped phase groups. The accuracy of this procedure is dependent on the accuracy of the mean estimated from coarse input D E M samples coming from the same pixels as the unwrapped phase group. Baseline parameters can be mis-estimated to compensate for errors in other parameters such as satellite altitude errors and absolute phase unwrapping errors. Precision orbit data or manual ground control point selection are therefore not required to perform accurate InSAR processing if a reasonably accurate coarse D E M is available to calibrate the interferometer. The baseline estimation algorithm is sensitive to linear trends and biases in the input coarse D E M . Processing as large an area as possible will minimize the 92 likelihood of the coarse D E M having spurious bias and trend errors. 93 Chapter 5 D E M Refinement Processing Experiments We report on experiments performed with two interferograms using the proposed technique: one generated from ERS Tandem Mission data and the other from a RADARSAT interferometric pair . The ERS Tandem Mission data has high coher1 ence with a mean coherence magnitude of approximately 0.7 and a (relatively) small normal baseline that makes this data set marginal for topographic height estimation and more suitable for differential interferometry. The RADARSAT interferogram has low coherence magnitude with mean value of approximately 0.35 and a larger normal baseline. The ERS data is in some respects an easier test case than the RADARSAT data because there are likely to be less phase unwrapping problems due to low coherence areas in the image. In addition, the RADARSAT interferogram has a much smaller (about five times) ambiguity height than the ERS data. When one models the interferogram using the coarse input D E M , D E M errors impact the D E M flattening algorithm more in the RADARSAT data case than in the ERS data 'We have also generated elevation data for another ERS test site, Sardegna, [46] characterized by high coherence magnitude but very difficult terrain. Lack of a complete reference DEM prevented full analysis of the results. 94 case. We therefore expect that the ERS dataset will respond better to the flattening algorithm than the RADARSAT data. However, note that the ERS InSAR heights will be more sensitive to phase noise because of the larger ambiguity height of the ERS dataset. Section 5.1 contains a description of the test site while Section 5.2 reviews the D E M data used during the processing. The experimental procedure is detailed in Section 5.3. The ERS results are reported in Section 5.4 and results for the RADARSAT-1 data are discussed in Section 5.5. Finally, conclusions drawn from the experimental results are given in Section 5.6. 5.1 Test Site Overview The Chilcotin area of British Columbia was chosen as a test-site because DEMs of 3 different qualities were available (TRIM, DTED-1, and GTOPO30 ) and because the topography is very challenging with large height variations and some layover in the ERS case. A collage of air-photos over the test site is shown in Figure 2 (5.1). The white rectangles on the collage show the approximate location of the SAR data investigated. Despite the rather low contrast of the scanned images, several features of the scene are clearly visible. The data covers a portion of the Fraser Canyon with elevation variation of about 1400 meters between the top of the canyon and the Fraser River below. The Chilcotin River flows in from the west to meet the Fraser River which flows from north to south. The darker areas of the figure likely correspond to forest or bush. This area is usually dry and there is limited vegetation so some portions of the scene can provide good repeat-pass interferometric coherence. 2 © 1 9 6 6 , 1 9 6 7 British Columbia Provincial Government 95 Figure 5.1: Aerial photos of the Chilcotin test site. The larger polygon shows the approximate boundaries of the ERS data set while the dashed lines show the approximate boundaries of the RADARSAT-1 data set. 96 5.2 D E M Data Overview The three DEMs used in this experiment were: T R I M D E M [52] The TRIM digital dataset is a recent product (1997) of the government of the province of British Columbia that consists of a D E M and possibly other cartographic information such as planimetric features and landcover. The TRIM D E M data were generated using photogrammetric techniques. Elevation data were captured manually in an approximate grid pattern with nominal data spacing of 100 m in flat areas and 75 meters in areas of higher relief. The height accuracy of the sampled data points is 5 m linear error at 90 % confidence (LE90) and 12 m circular error at 90 % confidence (CE90). DTED-1 DEM[56] The DTED-1 D E M is a discontinued product from the Government of Canada which was generated using 1950's era mapping technology. It is a gridded data set with sample spacing of 3 arcseconds or approximately 90 m. The accuracy goal of DTED-1 data is 30 m LE90 and 50 m CE90. Note that the quality of these discontinued Canadian DTED-1 data is known to vary widely. GTOPO30 DEM[21] The GTOPO30 D E M is a project of the United States Geological Survey (USGS) to provide a publicly available (via ftp) coarse, 3 low-resolution data set suitable for modeling large-scale processes. A similar data set is available from the GLOBE project . This D E M was generated 4 by combining various sources of D E M data to generate a composite D E M with data postings of 30 arc-seconds (approximately 1 km). Where possible, the GTOPO30 data set was generated from higher resolution D E M data. In this case the data was carefully processed to ensure that drainage patterns 3 4 See See http://edcwww.cr.usgs.gov/landdaac/gtopo30/gtopo30.html http://www.ngdc.noaa.gov/seg/topo/globe.shtml 97 in the GTOPO30 D E M were consistent with drainage patterns in the higher resolution data. The nominal accuracy according to the GTOPO30 D E M documentation is • ± 650 meters absolute linear error (vertical) at 90% confidence, • 2000 meters circular error (horizontal), which is quite poor. However, recent investigations to characterize the accuracy of the data have found that the average vertical accuracy of the whole data set is much better with an RMSE of approximately 70 meters [21]. The accuracy of the GTOPO30 varies with the underlying data sources and has an estimated RMSE of approximately 86 meters for the Chilcotin area[21]. The resolution and accuracy properties of the DEMs are summarized in Table 5.1. The three D E M data sets are plotted in Figure (5.2), Figure (5.3), and Figure (5.4). The DTED-1 and GTOPO30 datasets are plotted at their respective full resolutions while the scattered data points of the TRIM data were interpolated to a regular grid of approximately 75 m. The two rectangles on the figures denote the approximate coverage of the RADARSAT-1 and ERS 1/2 interferograms investigated. The larger rectangle shows the approximate location of ERS data processed while the smaller dashed line rectangle shows the approximate location of the RADARSAT data. Note the topographic variation of the terrain, especially in the main canyon of the Fraser River and the smaller canyon of the Chilcotin River. The difference in resolution and hence detail between the DEMs is clearly shown by these figures. We could only procure a small number of TRIM datasets so the TRIM data has limited coverage compared with the other two DEMs. Closer examination of a small area of the DEMs gives an idea of the relative magnitude of error characteristics of the data. For example, the region centered at 51.75° N and - 1 2 2 . 3 ° W has three different characteristics depending on which D E M 98 DEM TRIM [52] DTED[56] GTOPO30[21] Sample Spacing « 100 meters in flat terrain,m 75 meters in rough terrain 3 arc-seconds gridded ( « 90 m) 30 arc-seconds gridded ( « 1 km) Vertical Accuracy (meters) 5 (LE 90%) Horizontal Accuracy (meters) 12 (CE 90%) 30 (LE 90%) 50 (CE 90%) 86 (RMSE) N/A Table 5.1: Summary of D E M characteristics for the Chilcotin test site. is examined (see Figure (5.5) through Figure (5.7)). In the TRIM data, it is evident that there is a small "hill" centered at this location. In the D T E D data, the hill appears to be partially missing. Errors of this type can occur at the edges of areas processed in separate stages of the map production. In the GTOPO30 data, there is a large spike in the D E M values. It has a magnitude nearly 650 meters larger than the TRIM data. Since GTOPO30 are derived in part from D T E D data, it is possible that the missing portion of the hill in the D T E D data forces the GTOPO30 D E M to have this large value to keep the hydrology consistent with the higher resolution data. 99 -122.5 -122.4 -122.3 Longitude (W) in degrees -122.2 -122.1 Figure 5.2: TRIM D E M of the Chilcotin test site. Longitude (W) in degrees Figure 5.3: DTED-1 D E M of the Chilcotin test site. 100 ii6oo 51.95 i 1400 51.9 h 1000 S 51.8 p 51.75 1800 51.7 1600 1400 1200 51.55 51.5 -122.4 -122.3 -122.2 Longitude (W) in degrees -122.1 Figure 5.4: GTOPO30 D E M of the Chilcotin test site. 1I6OO 1400 1200 1000 1800 1600 1400 200 51.7 -122.34 -122.32 -122.3 -122.28 Longitude (W) in degrees -122.26 Figure 5.5: Zoomed plot of TRIM data. 101 51.7 -122.34 -122.32 -122.3 -122.28 Longitude (W) in degrees -12226 Figure 5.6: Zoomed plot of D T E D data. 51.7 -122.34 -122.32 -122.3 -122.28 Longitude (W) in degrees -122.26 Figure 5.7: Zoomed plot of GTOPO30 data. 102 5.3 Experimental Procedure The experiments implemented the algorithm defined in Appendix E . The baseline magnitude and baseline orientation were estimated using a linear model in azimuth. The TRIM, D T E D , and GTOPO30 data sets were each used in turn as the reference D E M for the interferogram calibration. Baseline parameters were also calculated from the orbit data of the satellites to provide another reference for comparison. The D E M and radar data were registered coarsely by matching the river in the SAR images and the bottom of the river valley in the resampled D E M . The registration between the D E M and SAR imagery was optimized by maximizing the coherence magnitude as a function of location and estimated baseline magnitude and orientation. The D E M flattening algorithm was then used to produce the residual interferograms. A four stage procedure was used to unwrap the filtered and subsampled residual phase of the interferogram: 1. Weighted least-squares phase unwrapping [57], 2. Unwrapping the residual phase from the least-squares procedure using either region-growing[58] or Flynn's minimum discontinuity procedure[54], 3. Phase offset checking by checking mean height offsets of unwrapped phase groups. 4. Interpolation of small "holes" in the unwrapped phase surface. The first two stages of the phase unwrapping procedure are necessary to deal with the bias problems of weighted least squares phase unwrapping (see Appendix C). Depending on the imaged scene content, the phase may be unwrapped into disconnected groups of unwrapped phase and phase offset errors may be introduced during the previous two phase unwrapping steps. These phase offset errors must be 103 corrected to generate a contiguous output unwrapped phase surface. Finally, small residual holes in the unwrapped interferogram phase are interpolated to generate a continuous output height estimate. "Larger" holes in the unwrapped phase map are indicative of areas of invalid data which are not likely to be successfully interpolated; these are filled with "NaNs" in the output height mask. For the initial weighted least-squares phase unwrapping we use a binary mask to eliminate unreliable areas from the phase unwrapping processing. We generate the unwrapping mask by "or-ing" masks generated using the following criterion: Residues are generally indicative of local inhomogeneities in the interferogram phase which should be avoided during the phase unwrapping process. If there is an area where the residue density is locally very high, this mask serves to minimize the influence the area has on the weighted least squares unwrapped phase. Magnitude thresholds are used to eliminate areas of low magnitude (i.e. lowest 5 %) and high magnitude (i.e. highest 5%) backscatter from the unwrapping procedure. Areas of low backscatter correspond to areas with low coherence magnitude or radar shadow zones. Shadows are areas of invalid data that may be interpolated or calculated using SAR data collected from the opposite side of the target scene. Areas of high magnitude backscatter generally correspond to areas of steep slope or layover and should be avoided during phase unwrapping as the likelihood of phase unwrapping errors increases in these regions. Coherence magnitude is used to discriminate between low and high SNR data. A threshold is set to eliminate areas of low coherence from the data that are not already masked by the magnitude threshold mask. Low coherence in areas of reasonable backscatter generally corresponds to target areas that have changed significantly between the imaging passes. The interferogram 104 phase from such areas is not reliable and should be avoided during.phase unwrapping processing. For the ERS-1 data, the coherence magnitude is high enough that not many areas were masked out of the phase unwrapping processing with the exception of the rivers. The residual phase after weighted least-squares phase unwrapping for the ERS-1 data was almost always within an ambiguity of the data so that generation of the final unwrapped phase data was easily performed using a region-growing approach [58]. However, the RADARSAT-1 data has lower coherence magnitude in general due to the 24 day repeat orbit period. Note that there is an area in the lower right middle of the RADARSAT-1 data that has substantially lower coherence magnitude than the rest of the data. We found that the region growing approach was extremely slow so we applied Flynn's minimum cutline approach [54, 55] to unwrap the residual phase for this dataset. Finally, the D E M updating algorithm was performed including checks for missed ambiguities in the unwrapped phase. In the following subsections, we discuss the results of the D E M refinement algorithm applied to the ERS and RADARSAT data sets. 105 Incidence Angle: near swath far swath Slant Range: near range far range Swath Width: ground range slant range Wavelength: altitude: 20° 26° 833 km 873 km 100 km 40 km 0.0566 m 782 km Table 5.2: Nominal ERS SAR parameters [6]. 5.4 E R SD E M Refinement Experiments The reference SAR image, the interferogram phase, and the coherence magnitude of the ERS dataset are shown in Figure (5.8). All images are oriented with north approximately at the top of the page. The ERS Tandem mission data are an ERS-1 SAR image (orbit 21730) and a co-located ERS-2 SAR image (orbit 2057), collected on September 10, 1995 and September 11, 1995 respectively at a local time of approximately 10 A M PST. The ERS parameters are given in Table 5.2. The SAR data was collected.on track 199 at frame 2565 for both satellites. The ERS InSAR baseline listing [42] estimates the normal baseline to be -42 meters, well within the desirable range for SAR interferometry. The estimated contribution of the baseline parameter errors to the output InSAR D E M error are shown in Table 5.3 while the estimated baseline parameter accuracies are shown in Figure (5.9). 106 DEM TRIM DTED-1 GTOPO30 #of Azimuth Samples #of Range Samples 0~mean in m in m in m 330 367 34 147 233 23 0.01 0.06 3.3 0.02 0.1 5.7 0.1 0.43 21.2 ^trend Table 5.3: Summary of estimated baseline parameter errors contributions to output InSAR D E M accuracy for ERS data. The accuracy of an output InSAR D E M estimated using the input coarse D E M was calculated using statistics of the input DEMs. 107 (a) ERS SAR image (b) Raw ERS phase (c) ERS coherence magnitude (d) Conventional "flattened" ERS phase Figure 5.8: Summary of ERS Tandem mission processed data set. The "flattened" interferogram phase was generated by removing the dominant range frequency from the raw interferogram phase. 108 (a) TRIM D E M parameter errors. (b) DTED-1 D E M parameter errors. (c) GTOPO30 D E M parameter errors. £ure 5.9: Expected parameter errors for input DEMs for ERS data. 109 5.4.1 ERS D E M Flattening Experiments The slant range registered DEMs for the ERS data are shown in Figure (5.10). These DEMs were each used with the full interferogram shown in Figure (5.8) '(b), to generate the baseline parameters summarized in Figure (5.11). The interferogram models generated using the baseline parameters and the initial coarse DEMs are shown in Figure (5.12). The filtered residual interferograms representing the filtered differences between the interferogram models and the full interferogram phase are shown in Figure (5.14) . 5 The residual interferograms are consistent with the quality of the input DEMs. The higher quality TRIM and DTED-1 data sets have less residual phase signals than the GTOPO30 input D E M . We examined the spectra of the residual interferograms and found that the DTED-1 and TRIM data sets could provide accurate baseline estimates. However, the GTOPO30 dataset had several significant peaks near zero frequency so we did not expect the GTOPO30 dataset to yield accurate baseline parameter estimates. The TRIM residual interferogram has minimal residual phase compared with the DTED-1 and GTOPO30 residual interferograms. There is an obvious correlation between the features in the input TRIM D E M and the residual interferogram. The question is whether or not the features in the residual interferogram come from the D E M or from the interferometric SAR data? To investigate this, we plotted the full resolution DEMs and residual interferograms for the TRIM and DTED-1 data in Figure (5.16) and Figure (5.17) respectively. The pixel sample along the ground down the page (in the azimuth direction) is constant and approximately 4 meters. The pixel spacing across the page (the slant range direction) varies as a function of surface slope. Going down into the canyon (passing from the right middle to the T h e raw differences between the interferogram models and the interferometry data are shown in the appendix in Figure (5.13). 5 110 left hand side of the page), the pixel sample spacing is on the order of 10-12 meters while in the flatter regions it is on the order of 22 meters. From comparing the TRIM residual D E M and the DTED-1 D E M it is clear that the input TRIM D E M is having an effect on the residual interferogram; particularly going down into the canyon. Features are in general much broader in the TRIM residual interferogram than the DTED-1 residual interferogram. The input DTED-1 D E M lacks the detail of the input TRIM D E M so the DTED-1 residual interferogram is a better indicator of the information coming from the interferometric SAR data. The raw TRIM data sample spacing sloping regions is approximately 75 meters which is far larger than the sample spacing of the interferograms. Therefore, the correlation between the TRIM D E M features and the TRIM residual interferogram is likely caused by the smearing of topographic features sampled at a coarser sample spacing than the interferometric SAR data. There is also evidence of a slight mis-registration between the SAR data and the input TRIM D E M as the residual interferogram phase signatures which are correlated to the D E M features tend to negative to the top of the page and positive to the bottom. In any case, it is clear that the D E M features are much larger (have smaller bandwidth) than the residual interferogram phase features. Note the detail visible in the DTED-1 and GTOPO30 residual interferograms in Figure (5.14) which is not clear in the original DEMs. The InSAR data has roughly four times finer resolution than the DTED-1 D E M and far more improvement relative to the GTOPO30 D E M . Topographic features not present in the reference D E M are present in the residual interferogram. The residual interferogram is easy to unwrap, and can be used to improve the resolution and accuracy of the coarser DEMs as long as the phase noise is moderate. The baselines estimated by the algorithm are shown in Figure (5.11), as a function of azimuth sample number. The normal baseline is approximately 45 meters 111 which gives a 2ir ambiguity of approximately 200 m for the ERS Tandem data. The accuracy of the algorithms' baseline estimates were estimated using the predicted error of a bilinear fit to errors of the given data as shown in Table 5.3. These height estimation errors were converted to equivalent baseline errors as shown in Figure (5.9). The baseline error for the precision orbit data was estimated from the consistency of successive precision orbit estimates [17] as approximately 0.7 m [17]. The baselines estimated from each of the three DEMs differ from the precision orbit data by up to 1 meter and differ from each other by more than the predicted error shown in Table 5.3. The failure of the GTOPO30 result to match the other results comes from the existence of many peaks near zero frequency for the GTOPO30 residual interferogram. For the ERS test case, the disagreement between the DTED-1 and TRIM based results for the ERS Tandem mission data set was puzzling. Both the DTED-1 and TRIM data sets are fairly good quality DEMs (DTED-1 90 % significance level is 30 meters [56], TRIM 90 % significance level is 5 meters [52]). The disagreement is attributable to a significant trend in range in the difference between the DTED-1 and TRIM DEMs. Because the TRIM data were generated relatively recently (1997) under stringent accuracy constraints, the error trend is more likely in the DTED-1 data which was produced during cold-war era mapping. The average difference between the DEMs varied from -4 meters to -22 meters across the range swath. This "tilt" is likely due to errors in the DTED-1 data set as well as some misregistration in range between the DEMs. In addition, the difference between the input DTED-1 and TRIM data (after the error trend was removed) had standard deviation of approximately 34 meters. This is far larger than the approximately 18 meters standard deviation expected assuming that the DEMs are uncorrelated and can be modeled as independent Gaussian random variables. Further investigation of the difference in the results was performed on a patch of the interferogram with complete coverage in both DTED-1 and TRIM DEMs for 112 the ERS data. After correcting for the offsets, range trends, azimuth trends, and the range-azimuth trends in the difference between the TRIM and DTED-1 DEMs, the estimates were significantly closer as shown in Figure (5.15). The slight lingering difference in the results can be traced to the slight difference in sample mean of the subset of the D E M data actually processed. To double-check the modified baseline estimates, S(r) derived from the TRIM D E M data with the TRIM baseline parameters was used with the DTED-1 and corrected DTED-1 baseline parameters to calculate topographic heights. Using the modified D T E D baseline parameter estimates yielded a D E M which had a mean difference of 0.4 meters and a standard deviation of 0.8 meters compared to the TRIM D E M . Using the unmodified DTED-1 baseline parameters yielded a D E M with a mean difference of -15 meters and standard deviation of 12 meters compared to the TRIM D E M . 113 (a) Input TRIM D E M . (b) Input DTED-1 D E M . (c) Input GTOPO30 DEM. Figure 5.10: Input DEMs for ERS interferograms in SAR coordinates. The data are approximately 33 km in azimuth (down the page) and 21 km in range (across the page from right to left). 114 46.8 -Q TRIM DTED GT0P03C| Orbit S.6ir : 46.2 45.8 " v. ' "V ^ v.. 45.6 Is 45.4 45.2 0 200 400 600 800 1000 1200 Azimuth Index 1400 1600 1800 2000 (a) Estimate of B in m as a function of azimuth index. 1.88 1.86 <t—e—e—e—e—e—e- o o 6 c - -o b p o TRIM DTED GTOPO30 Orbit 1.84 1.82 1.8 1.78 1.76 0 200 400 600 800 1000 1200 Azimuth Index 1400 1600 1800 2000 (b) Estimate of 0 in rads as a function of azimuth index. Figure 5.11: Results of the baseline estimation algorithm applied with GTOPO30, DTED-1, and TRIM data over the Chilcotin test site for ERS Tandem mission data. 115 (a) TRIM data used for reference D E M . (b) DTED-1 data used for reference D E M . (c) GTOPO30 data used for reference D E M . Figure 5.12: Model interferograms using estimated baseline parameters from ERS processing. 116 (c) GTOPO30 data used for reference DEM. Figure 5.13: Raw residual interferograms from ERS processing. 117 (a) TRIM data used for reference D E M . (b) DTED-1 data used for reference D E M . (c) GTOPO30 data used for reference D E M . Figure 5.14: Filtered residual interferograms from ERS processing. 118 (a) Estimate of B in m as a function of azimuth index. 1.769 DTED TRIM Corrected DTED 1.768 1.767 £ 1766 » 1.765 1.764 1.763 1.762 200 400 600 800 xindx 1000 1200 1400 1600 1800 (b) Estimate of 0 in rads as a function of azimuth index. Figure 5.15: Comparison between DTED, trend corrected D T E D , and TRIM baseline estimates for the ERS dataset. 119 TRIM input D E M . Flattened residual interferogram Figure 5.16: Full resolution comparison of input TRIM D E M and filtered residual interferogram phase. 120 DTED-1 input D E M . Flattened residual interferogram Figure 5.17: Full resolution comparison of input DTED-1 D E M and filtered residual interferogram phase. 121 5.4.2 ERS D E M Updating Results The input and output refined DEMs for the ERS-1 data based on TRIM, DTED-1, modified DTED-1, GTOPO30, and modified GTOPO30 are shown in Figure (5.19) through Figure (5.23). The "modified" D E M trials were performed after removing biases and error trends with respect to the input TRIM D E M from the input D T E D 1 and GTOPO30 DEMs. In the case of the GTOPO30 input D E M , the modified D E M also includes a correction for a group phase unwrapping error. The gray scale for all the D E M images goes from -50 meters to 1350 meters. The light areas in the output DEMs are areas where valid height samples could not be calculated. The baseline parameter estimates for all input DEMs (including the modified ones) are shown in Figure (5.18). The baseline parameters are consistent with the analysis and results of previous chapters whereby baseline parameters depend on the bias and trend errors in the input D E M . The modified D E M baseline parameters are much closer to the TRIM baseline parameters except for the GTOPO30 baseline parameters. The residual differences of the GTOPO30 baseline parameters seem to be a result of a combination of the correlation of D E M errors with position and phase unwrapping errors. The GTOPO30 errors are not uniformly distributed with a large outlier at the center of the D E M and a bias to large negative errors in the lower left portion of the D E M . The height errors of each of the input and output InSAR DEMs are summarized in Table 5.4 . 6 The error statistics were calculated from the difference of the D E M of interest with the TRIM input D E M . The height error statistics were accumulated on all the valid pixels in the DEMs; including output D E M values that were interpolated. The 90 th percentile error is the mean corrected D E M error that was greater than 90% of the absolute mean error corrected D E M errors. All DEMs are labeled as being input or output. The input DEMs are the DEMs used in the 6 T h e histograms are in Appendix F from Figure (F.l) to Figure (F.9) 122 flattening processing of the data. The InSAR TRIM result of approximately 26 meters represents the noise threshold of the InSAR height estimation procedure as the baseline parameter errors for this input D E M should be negligible (see Table 5.3). The relatively high height error is because of the small baseline and hence large ambiguity height which makes the interferogram D E M sensitive to phase noise. The two modified output DEMs are close to the 26 meter result from the TRIM data set. The GTOPO30 results in Table 5.4 demonstrate that substantial improvement in D E M quality can be gained by InSAR processing subject to presence of biases in the input coarse D E M . Bias errors in the input coarse D E M are represented in the output InSAR D E M as expected. This is most clearly seen in the 90% error which tends to be roughly the same order of magnitude for the uncorrected input DEMs and the output InSAR DEMs derived from them. The algorithm for equalizing the constant phase offset between different unwrapped groups of pixels failed for the GTOPO30 dataset as shown in Figure (5.24) (b). The errors on the left hand side of the river all tend to be much lower than the rest of the D E M in the original input D E M shown in Figure (5.24) (a). When the equalization process takes place, the locally low bias for the errors on left hand side of the river provoke an error in the relative phase offset in the output D E M that is clearly shown in the error image . 7 The effect of the baseline parameter differences is shown graphically in Figure (5.25) where the differences between the InSAR-TRIM output D E M and the modified DTED-l-InSAR and GTOPO30-InSAR output DEMs is shown. The only differences here are differences due to the baseline parameters and also values of the unwrapped phase. The gray scale for these plots goes from -50 to 50 meters. The InSAR DTED-1 D E M is relatively close to the output TRIM D E M over most of See the double peak in the histogram in the appendixfigureFigure (F.7) also. 7 123 DEM InSAR TRIM DTED-1 InSAR DTED-1 Modified DTED-1 InSAR Modified D T E D GTOPO30 InSAR GTOPO30 Modified GTOPO30 InSAR Modified GTOPO30 Input/Output mean error (m) Output Input Output Input Output Input Output Input Output -0.83 -18.58 -16.67 -1.1 -1.25 9.72 3.85 -0.9 -0.9 a (m) 26.1 34.3 29.1 32.0 26.9 . 115.8 60.9 108.6 29.3 90*"% (m) 42.1 61.2 52.9 52.6 43.2 158.6 104.8 156.0 48. 3 Table 5.4: Summary of results for ERS topography estimation. the image with some small differences due to phase unwrapping errors (the darkest and brightest pixels). However, the output InSAR-GTOPO30 D E M shows signs of significant differences across the swath, on the order of 35 meters which is larger than the predictions made in Table 5.3. Even after correcting for the bias and trend errors in the whole GTOPO30 DEM, the errors are across the D E M are correlated by location with large errors along the vertical centerline of the D E M and large negative errors to the left hand (west) side of the river. 124 47.51 1 1 \ 1 1 ! 1 r X o-o-o-o -o &• ex 46 </) -e- e e e - 3 - G - 0 - <>o-o-o- x S X "X 15 X 'X E 45.5 - "X x 'X --B- o- 0 200 TRIM DTED-1 Modified DTED-1 GTOPO30 Modified GTOP03 ) 400 600 x "X 'X X 800 1000 1200 1400 1600 1800 azimuth index (a) Baseline magnitude (B) estimates. 1.82 e- o o - (> - o - o - o - €1 O - O - O -©' o - o - o - © " - & •0' ^ - e -©-43- c 1.81 1.79 TRIM DTED-1 Modified DTED-1 GTOPO30 Modified GTOP03 ) --Q- o- - 1.78 —& o • 1.76 0 200 •n i n 400 C L Q 600 Q UJ .0 800 • 1000 a n a n CJ UI 1 1200 _ LJ 1400 1600 azimuth index (b) Baseline orientation (0) estimates. Figure 5.18: Summary of baseline parameter estimates 125 1800 (a) Input TRIM D E M . (b) Output InSAR-TRIM D E M . Figure 5.19: Results of InSAR-TRIM D E M Experiment with ERS Tandem data. 126 (a) Input DTED-1 D E M . (b) Output InSAR-DTED-1 D E M . Figure 5.20: Results of InSAR-DTED-1 D E M experiment with ERS Tandem data. (a) Input Modified DTED-1 D E M . (b) Output InSAR-Modified DTED-1 D E M . Figure 5.21: Results of InSAR-Modified DTED-1 D E M experiment with ERS Tandem data. 127 (a) Input GTOPO30 D E M . Figure 5.22: data. (b) Output InSAR-GTOPO30 D E M . Results of InSAR-GTOPO30 D E M experiment with ERS Tandem (a) Input Modified GTOPO30 D E M . (b) Output InSAR-Modified GTOPO30 D E M . Figure 5.23: Results of InSAR-Modified Tandem data. GTOPO30 D E M experiment with ERS 128 (a) Input GTOPO30 D E M Errors, (b) Output InSAR GTOPO30 D E M Errors. Figure 5.24: Errors of input and output GTOPO30 DEMs Incidence Angle: near swath far swath Slant Range: near range far range Swath Width: ground range slant range Wavelength: altitude: 41.6° 44.2° 1020 km 1060 km 50 km 40 km 0.056 m 807 km Table 5.5: Nominal RADARSAT Fine Beam 3 SAR parameters. 5.5 R A D A R S A T D E M Refinement Experiments The reference SAR image, the interferogram phase, and the coherence magnitude of the RADARSAT dataset is shown in Figure (5.26). The RADARSAT data covers approximately 10 km in range and 30 km in azimuth (or about 12% of a RADARSAT-1 Fine beam dataset) All images are oriented with north approximately at the top of the page. The RADARSAT processing example was taken from data collected in "FINE 3" beam mode, with the parameters given in Table 5.5. The two passes were collected on April 24, 1997 (orbit 7675) and May 18, 1997 (orbit 8018) respectively at a local time of approximately 6 A M PST. The normal baseline was calculated geometrically from RADARSAT's orbit data to be approximately 200 meters which is well within the preferred range for RADARSAT fine beam data. The estimated contribution of the baseline parameter errors to the output InSAR D E M error are shown in Table 5.6 while the estimated baseline parameter accuracies are shown in Figure (5.27). 130 DEM TRIM DTED GTOPO30 #of Azimuth Samples #of Range Samples &mean in m in m in m 300 333 30 100 111 10 0.02 0.1 5.0 0.03 0.3 8.6 0.12 0.65 30.9 &trend Table 5.6: Summary of estimated baseline parameter errors' contributions to output InSAR D E M accuracy for RADARSAT data. The error estimates were calculated using the error statistics of Table 5.1. 131 (a) RADARSAT-1 SAR Image (b) Raw RADARSAT-1 phase. (c) RADARSAT-1 Coherence Magnitude (d) Conventional "flattened" RADARSAT-1 phase. Figure 5.26: Summary of RADARSAT processed data sets. The "flattened" interferogram phase was generated by removing the dominant range frequency from the raw interferogram phase. 132 x 10" 3 -1.5 -1 -0.5 0 A S in fads 0.5 1 1.5 (a) TRIM D E M parameter errors. -0.031 -0.02 x_ X. X. x_ X. X, X X, 0.02 -1 0 A 8 In rada 1 (b) D T E D D E M parameter errors. (c) GTOPO30 D E M parameter errors. Figure 5.27: Estimated parameter errors for input DEMs for RADARSAT data. 133 5.5.1 RADARSAT D E M Flattening Experiment The slant range registered DEMs for the RADARSAT data are shown in Figure (5.28). These DEMs were used with the raw interferogram phase shown in Figure (5.26) (c), to generate the baseline parameters summarized in Figure (5.29). The interferogram models generated using the baseline parameters and the initial coarse DEMs are shown in Figure (5.30). The raw differences between the interferogram models and the interferometry data are shown in Figure (5.31). Finally, the filtered residual interferograms are shown in Figure (5.32). There is a wide spread between the baseline parameters estimated using the different DEMs and the orbit data. Both the D T E D and GTOPO30 baseline parameter estimates can be considered invalid because there is no dominant peak at zero frequency in both the range and azimuth spectra (see Figure (5.33) for the azimuth spectrum). The TRIM D E M baseline parameter estimates are expected to be accurate based on the absence of any sizeable peaks at non-zero frequencies in the spectrum of the residual interferogram. Despite the failure of the algorithm to converge using the D T E D and GTOPO30 reference DEMs, these DEMs allow a substantial "flattening" of the interferogram which can benefit further InSAR processing. For example, because the D E M - flattened interferogram has a reduced bandwidth compared to the original interferogram, a greater degree of phase smoothing (band-pass filtering) can be applied to make phase unwrapping more reliable. This is illustrated in Figure (5.31) and Figure (5.32), where residual interferograms, before and after band-pass filtering, are shown. Compared with the flat-earth flattened interferogram in Figure (5.26), it can be seen how the interferograms in Figure (5.32) are much easier to unwrap. Note however, that bandpass filtering may also eliminate detail so the amount of filtering must be optimized to suppress noise while conserving detail. 134 (c) GTOPO30 D E M . Figure 5.28: Input DEMs for RADARSAT data in SAR coordinates. 135 3 j e 1PRIM [3TED — (3TOP03C (Drbit 0 x - T f " ^s r e e— w o 365 360 0 200 400 600 800 Azimuth Index 1000 1200 1400 (a) Estimate of B in m as a function of azimuth index. I o > e (j e e3 0 «3 ^ «j a c3 I o TRIM DTED GTOP03C Orbit C 1.24 1.22 •-• 1.2 '0 200 400 600 800 Azimuth Index _ 1000 1200 1400 (b) Estimate of 0 in rads as a function of azimuth index. Figure 5.29: Results of the flattening algorithm applied with GTOPO30, DTED-1, and TRIM data over the Chilcotin test site for RADARSAT data. 136 (a) TRIM data used for reference D E M . (b) D T E D data used for reference D E M . (c) GTOPO30 data used for reference D E M . Figure 5.30: Model interferograms using estimated baseline parameters from RADARSAT-1 processing. 137 (b) D T E D residual interferogram. (a) TRIM residual interferogram. (c) GTOPO30 residual interferogram Figure 5.31: Residual interferograms from processed RADARSAT data. 138 (a) Filtered TRIM residual interferogram. (b) Filtered D T E D residual interferogram. (c) Filtered GTOPO30 residual interferogram Figure 5.32: Filtered residual interferograms from processed RADARSAT data. 139 (a) TRIM reference D E M (b) D T E D reference D E M (c) GTOPO30 reference D E M Figure 5.33: Azimuth spectra of residual RADARSAT interferograms from each of the processed input D E M datasets. 140 5.5.2 RADARSAT D E M Updating Results The input and output refined DEMs for the RADARSAT-1 data based on TRIM, DTED-1, modified DTED-1, GTOPO30, and modified GTOPO30 are shown in Figure (5.35) through Figure (5.39). As with the ERS data, the "modified" D E M trials were performed after removing biases and error trends with respect to the input TRIM D E M from the input DTED-1 and GTOPO30 DEMs and also correcting for some phase unwrapping errors in the GTOPO30 case. The gray scale for all the D E M images goes from -50 meters to 1100 meters. The dark areas in the InSAR DEMs are areas for which valid height samples could not be calculated. The height errors are summarized in Table. 5.7 . The error statistics were 8 calculated from the difference of the D E M of interest with the TRIM input D E M . The D E M error statistics were accumulated on all the valid pixels in the DEMs; including output D E M values that were interpolated. The 90 th percentile error is the mean corrected D E M error that was greater than 90% of the absolute mean error corrected D E M errors. All DEMs are labeled as being input or output. The input DEMs are the DEMs used in the flattening processing of the data. The DTED-1 results for the RADARSAT data have similar characteristics to the DTED-1 results of the ERS Tandem Mission data set in that the majority of the data was unwrapped correctly and reasonable improvement to the initial input D E M was possible. After removing the biases and error trends in the results, it is apparent that the output D E M approaches the TRIM data quality and agrees with TRIM data to within DTED-2 standards. Note that the compared to the ERS results, there are more phase unwrapping errors in the RADARSAT-1 data as compared with the ERS data. In particular, there is one area in the lower left of the output InSAR D E M shown in the error plot of Figure (5.40) (b) that is not T h e histograms for each of the input and output InSAR DEMs are shown in Figure (G.l) through Figure (G.9) 8 141 unwrapped correctly although it is part of a larger contiguous group on the left hand (west) side of the river canyon. However, this area actually has a large error in the input D E M as shown in Figure (5.40) (a) also. Furthermore, the mask used to exclude noisy pixels from the phase unwrapping process actually isolates the error from the surrounding pixels on three sides so there is no path to take to unwrap this area correctly. The GTOPO30 dataset does not perform as well. There are a number of different phase groups which are not unwrapped contiguously. These errors are due to the coarse nature of the input D E M which generates correlated errors that are biased low on the left (west) side of the Fraser River. After removing the phase unwrapping errors and the biases and trend errors in the initial D E M , we get results which are reasonably close to those of the input TRIM D E M . The baseline parameter estimates for all input DEMs (including the modified ones) are shown in Figure (5.34). The baseline parameters are consistent with the results of the analysis of the previous chapters whereby baseline parameters depend on the bias and trend errors in the input D E M . The modified D E M baseline parameters are much closer to the TRIM baseline parameters. The effect of the baseline parameter differences is shown graphically in Figure (5.41). The differences between the InSAR-TRIM output D E M and the modified DTED-l-InSAR and modified GTOPO30-InSAR output DEMs is shown. The only differences here are differences due to the baseline parameters and also values of the unwrapped phase. The gray scale for these plots goes from -20 to 20 meters. The InSAR DTED-1 D E M is relatively close to the output TRIM D E M over most of the image with some small differences due phase unwrapping errors (the darkest and brightest pixels). However, the output InSAR-GTOPO30 D E M shows signs of significant differences across the swath, on the order of 20 meters which is within reasonable agreement with the predictions made in Table 5.6. As with the ERS case, the differences for the GTOPO30 D E M seem to come from the large outlier in the center of the D E M 142 and the area of negative large D E M errors on the left side (west) of the river. The distribution of the D E M errors in the GTOPO30 input D E M are not Gaussian and they appear to be correlated with position. We therefore do not expect to get reliable height estimates from the D E M updating algorithm. DEM InSAR TRIM DTED-1 InSAR DTED-1 Modified DTED-1 InSAR Modified DTED-1 GTOPO30 InSAR GTOPO30 Modified GTOPO30 InSAR Modified GTOPO30 Input/Output mean error (m) cr (m) 90 % (m) Output Input Output Input Output Input Output Input Output -0.4 -16.4 -17.0 -0.1 -0.6 7.83 6.8 0.1 -0.8 12.2 35.6 18.8 34.5 14.3 98.3 63.8 91.1 19.0 18.6 61.3 37.0 56.9 19.7 143.7 112.2 134.2 20.8 th Table 5.7: Input D E M and InSAR output D E M error statistics (m) derived from comparison of DEMs with the reference TRIM dataset. 143 384 383 382 381 £ 380 CD 379 378 TRIM DTED-1 - -a 377 376 Modified DTED-1 GTOPO30 • - o- • Modified GTOP03 1000 500 1500 azimuth index 2000 2500 3000 2500 3000 (a) Baseline magnitude (B) estimates. 1.295 1.285 1.28 1.275 1.27 TRIM 1.265 1.26 DTED-1 -Q — Modified DTED-1 GTOPO30 Modified GTOPO30 500 1000 1500 azimuth index 2000 (b) Baseline orientation (0) estimates, figure 5.34: Summary of baseline parameter estimates for RADARSAT-1 data. 144 (a) Input T R I M D E M . (b) Output I n S A R - T R I M D E M . Figure 5.35: Results of I n S A R - T R I M D E M Experiment with R A D A R S A T - 1 data. 145 (a) Input DTED-1 D E M . (b) Output InSAR-DTED-1 D E M . Figure 5.36: Results of InSAR-DTED-1 D E M experiment with RADARSAT-1 data. (a) Modified Input DTED-1 D E M . (b) Output InSAR-DTED-1 (Modified) D E M . Figure 5.37: Results of InSAR-Modified RADARSAT-1 data. 146 DTED-1 D E M experiment with (a) Input GTOPO30 D E M . (b) Output InSAR-GTOPO30 D E M . Figure 5.38: Results of InSAR-GTOPO30 D E M experiment with RADARSAT-1 data. (a) Input GTOPO30 Modified D E M . (b) Output InSAR-GTOPO30 (modified) D E M . Figure 5.39: Results of InSAR-Modified RADARSAT-1 data. 147 GTOPO30 D E M experiment with (a) Input DTED-1 D E M Errors, (b) Output DTED-1 D E M Errors. Figure 5.40: Errors of input and output D T E D DEMs (a) InSAR DTED-1 output. (b) InSAR GTOPO30 output. Figure 5.41: Difference of output DTED-1 and GTOPO30 InSAR DEMS with respect to the output InSAR-TRIM DEMs. Gray sale goes from -20 m to +20 m. 148 5.6 Summary of Experimental Results The D E M refinement algorithm was applied to ERS and RADARSAT-1 interferometric data sets using good (TRIM), moderate (DTED-1) and poor quality (GTOPO30) DEMs. Results for the D E M flattening algorithm and D E M updating algorithm were analyzed separately for each of the two datasets. In all cases, the D E M flattening algorithm was able to provide a residual interferogram with substantially less phase variability than the original data which facilitated phase unwrapping through filtering and subsampling of the original data. The accuracy of the estimated baseline parameters was strongly dependent on the input D E M quality versus the sensitivity of the interferometer to height changes. From the spectra of the residual interferogram, we saw that the ERS experiments with the GTOPO30 data and the RADARSAT experiments with the D T E D and GTOPO30 data could not yield valid baseline parameter estimates because of the large number of spectral peaks. The sensitivity of the D E M flattening algorithm to input D E M bias and trend errors was demonstrated in the ERS results for the D T E D and TRIM DEMs. The D T E D D E M used in the experiments differed from the TRIM D E M by an approximately linear trend which in the ERS case led to substantially different baseline estimates although the algorithm converged for both input DEMs. After removing the relative trend from the D T E D D E M , the baseline estimates were much closer to the TRIM parameter values as expected. Uncorrelated noise in the D E M or interferogram is not a significant problem for the flattening algorithm if the data is suitable for interferometry because of the large number of samples available for estimating each term in the residual phase eqn. (3.6). For the RADARSAT dataset which has relatively low SNR, the D E M flattening algorithm converged to a stable estimate of the baseline and provided a flattened residual interferogram for TRIM data set. The flattening algorithm is robust in the sense that it will always give aflattenedD E M that can be demodulated 149 and filtered in preparation for phase unwrapping. However, the baseline estimates derived from the flattening algorithm are likely to be in error and should be improved by fitting the InSAR D E M to the coarse low-quality input D E M . The D E M updating algorithm was successfully applied to most datasets except for the GTOPO30 data applied to the RADARSAT SAR data. For the RADARSAT-1 data, we were able to substantially improve the quality of the input DTED-1 (roughly halving the input standard deviation) in both the modified and unmodified D E M cases. For the GTOPO30 dataset, there were phase unwrapping discontinuities between different unwrapped phase groups which caused a major failure in the algorithm's operations. After the phase unwrapping errors had been corrected and the input D E M compensated for bias and trend errors, the results were on the order of a factor of 4 better than in the input D E M . For the ERS Tandem mission data, the results were similar though less impressive because of the small normal baseline which makes the ERS data very sensitive to noise. D E M quality improvement was marginal for the DTED-1 data and about a factor of five for the GTOPO30 data. 150 C h a p t e r 6 Conclusions We have developed an algorithm for including coarse, relatively low quality DEMs in the InSAR processing stream. The algorithm consists of 3 parts: 1. "Flattening" using a coarse, low-quality D E M . 2. Phase unwrapping of residual phase signal. 3. D E M updating of the coarse, low-quality D E M . In the first stage, the algorithm estimates the topographic phase associated with the coarse D E M without phase unwrapping by manipulating the baseline parameter estimates using non-linear optimization. The phase "flattening" afforded by this operation reduces the complexity of the phase unwrapping operation by minimizing the variation in the residual interferogram phase as a function of the coarse D E M data quality. In addition, this procedure facilitates simple filtering of the residual interferogram to ease the difficulties with phase unwrapping processing. In the last stage of the algorithm, the coarse input D E M is used as a template to re-estimate the baseline parameters using the unwrapped phase. In this way, the residual phase of the interferogram can be used to refine the quality of the input D E M by adding 151 detail at some small integer multiple of the SAR sample spacing. The accuracy of the topographic estimates of the algorithm are a function of the input DEM's quality. In particular the two baseline parameter estimation algorithms (one for D E M flattening, the other for D E M updating) are both sensitive to linear trends and biases in the input D E M . As the input D E M becomes coarser and of lower and lower quality, the likelihood of developing a spurious error trend in the output D E M becomes higher for a given data extent. It is therefore advantageous to process as large an input D E M area as possible. The performance of the D E M refinement algorithm was investigated for two data sets: one with high SNR (the ERS Tandem Mission data), the other with low SNR (the RADARSAT-1 data). The initial baseline parameter estimation without phase unwrapping facilitated aggressive lowpass filtering of the low SNR RADARSAT-1 data which resulted in DTED-2 quality topographic height accuracies as compared to the TRIM ground truth dataset. Substantial improvement in the InSAR D E M quality was demonstrated for input DTED-1 and GTOPO30 DEMs for modified input DEMs which had spurious trends removed. It is worthwhile noting that the ERS Tandem Mission data set is about 7% of a standard ERS scene (100 km by 100 km) while the RADARSAT data set was about 12% of a standard Fine beam scene (50 km by 50 km). There is substantial room for improvement in the algorithm results by processing larger data sets. Note that the baseline parameters can be mis-estimated to compensate for errors in other parameters such as phase offset and satellite position (altitude). We do not need accurate orbit data or accurate phase offset estimates to produce good quality topographic height estimates as shown by the RADARSAT-1 experiments. The baseline parameters can absorb the other errors to provide a good quality topographic estimate. 152 6.1 Thesis Contributions The contributions of this thesis are: • a method of baseline parameter estimation to facilitate using coarse, low quality DEMs for interferogram flattening, • a method for determining when the baseline parameters used in the flattening algorithm are accurate, • an analysis of the effect of D E M errors on the determination of the interferometer geometry, • an analysis of the effects of baseline parameter errors combined with other system errors on the output InSAR DEM, • a method for refining the estimate of the geometry of the interferometer after unwrapping the interferogram phase, • a method for calculating phase unwrapping group offsets to equalize independently unwrapped groups of interferogram phase, and • a method for predicting the lower bound on topographic accuracy of InSAR DEMs derived from coarse low-quality DEMs. 6.2 Future Work There are many possible directions to take with this work: • Testing with large datasets. The performance measured with the processed datasets was limited by computational resources. Better allocation of system resources and more use of disks for intermediate results would facilitate processing larger datasets. 153 Phase unwrapping using the D E M as a constraint. If the D E M is of fairly good quality (DTED-1), the D E M could act as a constraint on the unwrapped phase value. The results of this study used a crude method for defining relative phase offsets after phase unwrapping processing had been performed. Using the D E M to constrain the phase unwrapping processing itself might lead to better unwrapped phase estimates. Adaptive interferogram filtering using the input D E M as prior information. If the D E M is of good quality but not of high resolution, the interferometric SAR data can be used to sharpen the detail of the input D E M . In this case, simple low-pass filtering to minimize noise may actually smear some topographic details. A better approach might be to design adaptive bandpass filters with the input D E M determining the filter bandwidth. Investigating other optimization schemes such as absolute value minimization to minimize the effects of outliers in the input coarse D E M or due to phase unwrapping errors. Least-squares minimization was used to estimate the baseline parameters based on a model of independent, Gaussian errors. However, least-squares minimization is sensitive to outliers as it tends to fit to the mean value of all variables including outliers. An improvement on the technique implemented for this study would be using weighted least squares or minimization by absolute values in concert with a noise model more suited for coarse, low-quality D E M errors. Geocoding testing. The current software does not tie the SAR coordinate system to a well known cartographic reference system (i.e. geocoding). To perform geocoding an assessment of the placement accuracy of the SAR derived topographic map with respect to the existing input D E M would have to be performed. 154 • Addition of other information. Other information such as hydrologic maps showing the locations of rivers, lakes, and streams, coastal boundary maps, and thematic maps showing the scene cover could also be integrated into the D E M refinement algorithm to aid in constraining the optimization algorithm to areas of good data and also to aid in the geo-location accuracy of the output InSAR D E M product. 155 Bibliography [1] C. Prati, F. Rocca, A. Monti Guarnieri, and P. Pasquali, "Report on ERS1 SAR interferometric techniques and applications," tech. rep., Politecnico di Milano, Dip. di Elet., POLIMI, Piazza L. da Vinci 32, 20133 Milano Italy, June 1994. [2] D. Massonnet and H. Vadon, "Precision positioning of interferometric track," in Proc. of CEOS SAR Calibration Workshop, (Noordwijk NL), 1993. [3] D. Massonnet, "Producing ground deformation maps automatically: the DIAPASON concept," in IGARSS'97, (Singapore), 1997. [4] H.A.Zebker, C. Werner, P. Rosen, and S. Hensley, "Accuracy of topographic maps derived from ers-1 interferometric radar," IEEE Trans. Geoscience and Remote Sensing, vol. 32, pp. 823-836, July 1994. [5] D. Geudtner and M . Schwabisch, "An algorithm for precise reconstruction of InSAR imaging geometry: Application to "flat eath" phase removal, phaseto-height conversion, and geocoding of InSAR-derived DEMs," in EUSAR'96, (Konigswinter, Germany), pp. 249-252, 1996. [6] "ERS-1 system." ESA Publications Division, 1992. [7] R. Gens, Quality assessment of SAR interferometric data. PhD thesis, ITC: International Institute for Aerospace Survey and Earth Sciences, 1998. 156 [8] J. C. Curlander and R. N. McDonough, Synthetic Aperture Radar Systems and Signal Processing. Wiley Series in Remote Sensing, John Wiley and Sons. Inc., 1991. [9] G. Rufino, A. Moccia, and S. Esposito, "DEM generation by means of ERS tandem data," IEEE Trans. GRS, vol. 36, pp. 1905-1912, Nov. 1998. [10] P. Vachon, D. Geudtner, A. Gray, and R. Touzi, "ERS-1 Synthetic Aperture Repeat-pass Interferometry: Implications for RADARSAT," Canadian Journal of Remote Sensing, vol. 21, pp. 441-454, Dec. 1995. [11] A. Gabriel and R. M . Goldstein, "Crossed orbit interferometry: Theory and experimental results from SIR-B.," International Journal of Remote Sensing, vol. 9, no. 5, pp. 857-872, 1988. [12] E . Christensen, N. Skou, J. Dall, K. Woelders, J.H.Jorgensen, and S. Madsen, "EMISAR: An absolutely calibrated polarimetric L- and C-band SAR," IEEE Trans. GRS, vol. 36, no. 6, pp. 1852-1865, 1998. [13] H. Zebker and J. Villasenor, "Decorrelation in interferometric radar echoes," IEEE Transactions Geoscience and Remote Sensing, vol. 30, pp. 950-959, Sept. 1992. [14] H. Tarayre and D. Massonnet, "Effects of a refractive atmosphere on interferometric processing," in IGARSS'94, (Pasadena CA), pp. 717-719, 1994. [15] R. Goldstein, "Atmospheric limitations to repeat pass interferometry," Geophysical Research Letters, vol. 22, pp. 2517-2520, Sep. 1995. [16] A. Ferretti, C. Prati, and F. Rocca, "Multibaseline insar dem reconstruction: The wavelet approach," IEEE Trans, on Geoscience and Remote Sensing, vol. 37, pp. 705-715, March 1999. 157 [17] C. Reigber, "Impact of precise orbits on SAR interferometry," in FRINGE'96, (Zurich, Switzerland), pp. 223-232, ESA, 1997. [18] C. Clark and R. Michaud, "RADARSAT Mission Planning," in Proceedings of a Workshop on RADARSAT Data Quality, (St.-Hubert, Quebec, Canada), CEOS, February 1997. [19] U. Wegmiiller and C. Werner, "Gamma SAR processor and interferometry software," in Third ERS Scientific Symposium, (Florence Italy), pp. 1687-1692, 1997. [20] M . van der Kooij, B. Armour, J. Ehrismann, H. Schwichow, and S. Sato, "A workstation for spaceborne interferometric SAR data.," in 26th International Symposium on Remote Sensing if Environment, (Vancouver BC), pp. 330-33, Mar. 25-29 1996. [21] D. B. Gesch, "Accuracy assessment of a global elevation model using shuttle laser altimetry data," in IGARSS'98, Seattle 1998. [22] D. Massonnet and T. Rabaute, "Radar interferometry: Limits and potentials," IEEE Trans. GRS, vol. 31, no. 2, pp. 455-464, 1993. [23] P. Rosen, S. Hensley, H. Zebker, F . Webb, and E . Fielding, "Surface deformation and coherence measurements of Kilauea Volcano, Hawaii from SIR-C radar interferometry," Journal of Geophysical Research, vol. 101, pp. 23,109-23,125, October 1996. [24] L. C. Graham, "Satellite interferometer radar for topographic mapping," Proceedings of the IEEE, vol. 62, no. 6, pp. 763-768, 1974. [25] D. Massonnet, M . Rossi, C. Carmona, F . Adragna, G. Peltzer, K. Feigl, and T. Rabaute, "The displacement field of the landers earthquake mapped by radar interferometry," Nature, vol. 364, pp. 138-142, July 1993. 158 [26] R. Goldstein, H. Engelhardt, B. Kamb, and R. Frolich, "Satellite radar interferometry for monitoring ice sheet motion: application to an Antarctic ice stream," Science, vol. 262, pp. 1525-1530, December 1993. [27] K. E . Mattar, P. W. Vachon, A. L. G. Dirk Geudtner, I. G. Cumming, and M . Brugman, "Validation of ERS tandem mission SAR measurements of alpine glacier velocity," IEEE Trans, on GRS, vol. 36, pp. 974-984, May 1998. [28] J. Goodman, "Statistical properties of laser speckle patterns," in Laser Speckle and Related Phenomena (J. Dainty, ed.), pp. 9-75, Springer Verlag, 1984. [29] M . Seymour and I. G. Cumming, "Maximium likelihood estimation and Cramer-Rao bounds for SAR interferometry," Tech. Rep. CICSR-TR95-03, CICSR, University of British Columbia, CICSR, Rm. 289, 2366 Main Mall, Vancouver BC, Canada V6T 1Z4, 1995. A summary of this report was presented at IGARSS'94 in Pasadena, California. [30] F. K. Li and R. Goldstein, "Studies of multibaseline spaceborne interferometric synthetic aperture radars.," IEEE Transactions on Geoscience and Remote Sensing, vol. 28, no. 1, pp. 88-97, 1990. [31] R. Bamler and D. Just, "Phase statistics and decorrelation in SAR interferograms," in IGARSS'93, pp. 980-984, 1993. [32] D. Just and R. Bamler, "Phase statistics of interferograms with applications to synthetic aperture radar," Applied Optics, vol. 33, pp. 4361-4368, Jul. 1994. [33] F . Gatelli, A. Monti Guarnieri, F . Parizzi, P. Pasquali, C. Prati, and F . Rocca, "The wavenumber shift in SAR interferometry," IEEE Trans. Geoscience and Remote Sensing, vol. 32, pp. 855-865, July 1994. [34] M . Schwabisch and D. Geudtner, "Improvement of phase and coherence map quality using azimuth prefiltering: Examples from ERS-1 and X-SAR," in IGARSS'95, (Florence, Italy), pp. 205-207, 1995. [35] M . Catabeni, A. Monti-Guarneri, and F. Rocca, "Estimation and improvement of coherence in SAR interferometry," in IGARSS'94, (Pasadena CA), pp. 720722,1994. [36] U. Wegmiiller, C. Werner, D. Niiesch, and A. Sieber, "Interferometric signatures of temperate forest from ERS-1 SAR data," in IGARSS'94, (Pasadena CA), pp. 292-294, 1994. [37] J. Hagberg, L. Ulander, and J. Askne, "Repeat-pass SAR interferometry over forested terrain," IEEE Transactions on Geoscience and Remote Sensing, vol. 33, pp. 331-340, March 1995. [38] G. Davidson and R. Bamler, "Multiresolution phase unwrapping for SAR interferometry," IEEE Trans. GRS, vol. 37, pp. 163-174, Jan. 1999. [39] A. Ferretti, Generazione di mappe altimetriche da osservazioni SAR multiple. PhD thesis, Politecnico di Milano, 1997. [40] D. Massonnet, H. Vadon, and M . Rossi, "Reduction of the need for phase unwrapping in radar interferometry," IEEE Trans. GRS, vol. 34, no. 2, pp. 489497, 1996. [41] U. Spagnolini, "2-D phase unwrapping and instantaneous frequency estimation," IEEE Trans. Geo. Remote Sensing, vol. 33, pp. 579-589, May 1995. [42] G. A. Solaas, "ERS-1 interferometric baseline algorithm verification," Tech. Rep. ES-TN-DPE-OM-GS02 ver. 1.0, ESA, 1994. [43] D. Small, C. Werner, and D. Niiesch, "Baseline modelling for ERS-1 SAR interferometry," in IGARSS '93, (Tokyo), pp. 1204-1206, Aug 1993. [44] M . Schwabisch, "Large scale interferometric D E M generation using ERS tandem data: Example of the Czech Republic," in IGARSS'98, (Seattle, Wa), 1998. 160 [45] L. Scales, Introduction to Non-Linear Optimization. New York: SpringerVerlag, 1985. [46] M . Seymour and I. Cumming, "InSAR Terrain Height Estimation Using LowQuality Sparse DEMs," in ERS Symposium, (Florence, Italy), pp. 421-426, ESA, March 1997. [47] The Mathworks Inc., Matlab Reference Guide, July 1993. [48] D. C. Rife and R. R. Boorstyn, "Single-tone parameter estimation from discretetime observations," IEEE Transactions on Information Theory, vol. IT-20, pp. 591-598, September 1974. [49] B. G . Quinn and P. J . Kootsookos, "Threshold behaviour of the maximum likelihood estimator of frequency," IEEE Trans. Signal Processing, vol. 42, no. 11, pp. 3291-3294, 1994. [50] B. James, B. D. O. Anderson, and R. C. Williamson, "Characterization of threshold for single tone maximum likelihood frequency estimation," IEEE Trans. Signal Processing, vol. 43, no. 4, pp. 817-821, 1995. [51] H. Tarayre and D. Massonnet, "Atmospheric propogation heterogenities revealed by ERS-1 interferometry," Geophysical Research Letters, vol. 23, pp. 989-992, May 1996. [52] Surveys and Resources Mapping Branch, British Columbia Specifications and Guidelines for Geomatics. Ministry of Environment, Lands and Parks, Province of British Columbia, release 2.0 ed., January 1992. [53] D. S. Andrews, "Simplifying terrain models and measuring terrain model accuracy," Tech. Rep. TR-96-05, University of British Columbia, Dept. of Computer Science, May 1996. 161 [54] T. Flynn, "Two-dimensional phase unwrapping with minimum weighted discontuity," Journal of the Optical Society of America A, vol. 14, pp. 2692-2701, Oct. 1997. [55] D. C. Ghiglia and M . D. Pritt, Two-Dimensional Phase Unwrapping: theory, algorithms, and software. John Wiley and Sons, 1998. [56] D. of Defense, "Military Specification Digital Terrain Elevation Data (DTED)." [57] D. Ghiglia and L. Romero, "Robust two-dimensional weighted and unweighted phase wrapping that uses fast transforms and iterative methods," J. Opt. Soc. Am. A, vol. 11, pp. 107-117, January 1994. [58] W. Xu and I. Cumming, "A region growing algorithm for InSAR phase unwrapping," in IGARSS'96, (Lincoln, USA), pp. 2044-2946, May 1996. [59] I. Cumming, F. Wong, and R. Raney, "A SAR processing algorithm with no interpolation," in IGARSS'92, (Houston Texas), pp. 376-379, 1992. [60] R. Touzi, A. Lopes, J. Bruniquel, and P. Vachon, "Unbiased estimation of the coherence function for SAR interferometry," in IGARSS'96, (Lincoln, USA), pp. 662-664, May 1996. [61] E . Rodriguez and J. Martin, "Theory and design of interferometric synthetic aperture radars," IEE Proceedings-F, vol. 139, no. 2, pp. 147-159, 1992. [62] R. M . Goldstein, H. A. Zebker, and C. L. Werner, "Satellite radar interferometry: two-dimensional phase unwrapping.," Radio Science, vol. 23, no. 4, pp. 713-720, 1988. [63] C. Prati, M . Giani, and N. Leuratti, "SAR interferometry: A 2-D phase unwrapping technique based on phase and absolute values information," in IGARRS'90, 1990. 162 [64] B. Hunt, "Matrix formulation of the reconstruction of phase values from phase differences," J. Opt. Soc. Am., vol. 69, pp. 393-399, Mar 1979. [65] R. Bamler, N. Adam, G. Davidson, and D. Just, "Noise-induced slope distortion in 2-D phase unwrapping by linear estimators with application to SAR interferometry," IEEE Trans. GRS, vol. 36, pp. 913-921, May 1998. [66] M . Pritt, "Phase unwrapping by means of multigrid techniques for iterferometric SAR," IEEE Trans. Geoscience and Remote Sensing, vol. 34, no. 3, pp. 728-738, 1996. [67] M . Costantini, "A novel phase unwrapping method based on network programming," IEEE Trans. GRS, vol. 36, pp. 813-821, May 1988. [68] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling, Numerical Recipes The Art of Scientific Computing. New York: Cambridge University Press, 1986. [69] J. A. Rice, Mathematical Statistics and Data Analysis. Pacific Grove California: Wadsworth and Brooks/Cole Advanced Books and Software, 1988. [70] M . Srinath and P. Rajasekaran, An Introduction to Statistical Signal Processing With Applications. John Wiley and Sons. Inc., 1979. [71] W. Xu and I. Cumming, "A region-growing algorithm for InSAR phase unwrapping," IEEE Trans. GRS, vol. 37, pp. 124-134, Jan. 1999. 163 A p p e n d i x A Spectral Shift Theorem The monochromatic approach used to derive the relationship between interferogram phase and topography in Section 2.3 hides an important aspect of InSAR image data. SAR data spectra correspond to a finite bandwidth portion of the ground range data. When two sets of SAR image data are collected from spatially separated antennas, the ground range spectra corresponding to the SAR data will be shifted [33]. The degree of shift is controlled by the the normal baseline (B- ). 1 Interferograms are formed by complex-conjugate multiplication of two SAR images. Since conjugate-multiplication in the time/space domains corresponds to correlation in the frequency/wavenumber domains, forming an interferogram corresponds to correlating the two SAR image spectra. It is informative to examine a one-dimensional slant range model of the SAR image formation process as a function of sampled ground wavenumber [33]. The SAR system transfer function can be modeled in range as a demodulation followed by convolution with an impulse response function: J(r) = b(x(r))e- °***w(r), = b(x(r)) e~ °^ (A.l) jw juj 164 * *w(r), (A.2) = b(x(r)) e j k r r (A.3) * * w(r); where: b(x(r)) x(r) u 0 — ground backscattering coefficient, — ground range location, = 2TT f , SAR center frequency in radians/s, 0 2u> 0 4TT , k = ** = the convolution operator, w(r) = impulse response function. r c = — , slant range wavenumber A We assume that the ground backscattering function has no extent in the vertical direction. Physically, this corresponds to an arrangement of surface scatterers on the ground plane. The ground wavenumber (k ) sampled is (see Figure (A.l)) x sin 6. (A.4) In cross-track SAR interferometry, the second pass image is collected from a slightly different track which means that the off-nadir angle 9 is also slightly different. Differentially, one can compute the change in k as x Ak x = 8k 89 A6 = k cos 9 A9. r (A.5) One can project this change into the slant range frequency domain to calculate the equivalent range frequency shift as A/=-A-A0. (A.6) tan 0 Although both SAR images are collected at the same slant-range center frequency (f ), the ground range wavenumber sampled is different due to the change in off0 nadir angle. In the InSAR community, this concept is known as the spectral-shift 165 phenomenon [33]. Note that the change in look angle also "stretches" the ground range spectra in a slightly different manner for both SAR images. However, satellite SAR systems generally have a small system bandwidth (10's of MHz) compared with a microwave center frequency (1 to 10 GHz). Therefore one can consider the shift in ground wavenumber to be constant over the SAR system bandwidth. When forming an interferogram, the portion of the SAR image data in the two SAR images that correspond to the same ground wavenumber spectra contributes an impulse in the interferogram spectra at the local interferogram frequency. The non-overlapping portions of the two image spectra contribute noise to the interferogram phase. The noise contributed to the interferogram in this manner is known as baseline decorrelation because it is a function of the normal baseline size. With an estimate of the interferogram phase frequency one can effectively reduce the noise by prefiltering the SAR images before forming the interferogram. From geometry, one can derive the change in off-nadir angle as a function of the normal baseline magnitude: (A.7) r As the separation between the SAR images tracks increases, the spectral-shift in ground range frequency also increases until the point where the spectra are no longer correlated. This occurs when the spectral shift is equal to the system bandwidth (W) and the baseline at this point is known as the critical baseline. One can derive the critical baseline for SAR interferometry as j_ max B M^rtanf? ~ ~f Jo Artanf? = 7T *Pr • ( A - ) 8 The "spectral shift" approach arises because of the different mapping from ground range to slant range depending on the relative position of the SAR antennas. For areas of uniform slope, one can reduce the effects of baseline decorrelation by prefiltering the SAR images to eliminate the uncorrelated parts of the SAR image 166 Ground range Figure A . l : Diagram showing geometrical relation between ground wavenumbers and slant range wavenumbers. For a fixed set of ground scatterers and ground wavenumbers (kx), the slant range wavenumber (kr) is a function of angle of incidence. spectra as a function of the local frequency. Furthermore, volume scattering effects also contribute to a loss in coherence between the images [33] in a similar fashion because the vertical wavenumbers sampled by the SAR images used to form the interferogram are also shifted relative to one another by the different cross-track trajectories of the SAR systems. 167 A p p e n d i x B Outline of Satellite I n S A R A l g o r i t h m to Obtain Topography Although the basic principles of SAR interferometry can be concisely described, the algorithm is complicated by the effect of noise and terrain discontinuities in the interferometric observations. There are many different implementations of InSAR processors available [3, 44, 23, 20, 19] which we synthesize into a representative satellite InSAR processing algorithm in the following subsections. Detailed InSAR Processing Algorithm Description 1. Identify suitable datasets. One requirement of successful InSAR processing is that the two SAR images used to make the interferogram are correlated. In general this means that the images must satisfy two constraints: 168 • the normal baseline constraint as in Appendix A, • temporal decorrelation, i.e. the SAR data must not be decorrelated due to temporal changes in the backscatter, e.g. due to plants growing. 2. Process SAR signal (raw) data to SLC images. Careful SAR data processing is required to ensure that the interferometric phase is not compromised during the SAR image generation processing. The general requirement is that the SAR processor be "phase preserving" [59]. A study of the effect of processor errors on InSAR data is presented by Just and Bamler in [32]. 3. Oversample the SAR images. To avoid spectral wrap-around when forming the interferogram, it is necessary to oversample both SAR images. 4. Pre-filter images in azimuth. In two-pass interferometry, there is a high likelihood that the Doppler centroid of the two images will not be the same for both passes. This implies that the two SAR images have sampled different portions of the azimuth spectra. Uncorrelated spectral components generate noise in the interferogram. Azimuth filtering based on the Doppler centroids and antenna patterns of SAR system can be applied to increase the coherence magnitude of the data [34]. 5. Register the SAR images. For two pass satellite interferometry, the SAR images must be well-registered to ensure good coherence magnitude in the output interferogram. Because the data for the two SAR images are collected from nearly the same location in space, the main differences in relative pixel location for the two images is a range offset, a slight range stretch, and a slight rotation in the azimuth direction [1]. The SAR image registration can be effected by a linear stretch in the range direction followed by a range dependent azimuth shift. 169 The interferogram sensitivity to misregistration can be analyzed by investigating coherence magnitude as a function of misregistration in pixels (A) [1]: = (B.l) Generally, the images must be registered within 1/10 of an oversampled pixel to avoid significantly degrading the interferogram phase quality. There are a number of different methods of achieving the required accuracy. Except for the orbit data approach, the methods extract offset estimates between different chips extracted across the image and use these chips to create a registration model. The methods for registration are: Precision Orbit Data Approach [43] Precision orbit data are available for ERS 1/2 satellites from D L R or 1 TUDelft usually with time lags on the order of a few days. These data 2 are accurate enough to derive the relative position of the satellites using scene centers positions as a reference. Maximizing Correlation The image chip magnitudes are correlated to estimate the local chip offsets. This method is useful when the coherence magnitude is low. Maximizing Coherence Magnitude Mini-interferograms are formed at different offsets between the image chips. The offset with the maximum coherent sum is used to determine the chip offset. This method works well for interferograms without large frequency trends. Maximizing Fringe Frequency [11] Mini-interferograms are formed at different offsets between the image chips. The maximum peak spectra of each of the mini-interferograms is used to determine the chip offset. Performing the spectral estimation 'http://nng. esoc.esa.de/ers/ers.html http://deos.lr.tudelft.nl/orbits/ 2 170 is equivalent to compensating for the local phase trend in the data for coherence magnitude estimates. This method gives reasonable results in most cases but it is time-consuming. 6. Pre-filter the SAR images in range. Following the spectral-shift approach (Appendix A), the SAR images must be pre-filtered in range as a function of local interferogram range frequency to suppress uncorrelated parts of the image spectrum. In general, this is an iterative process which requires formation of the initial interferogram. One usual simplification is to filter the data to eliminate the positive or negative frequencies according to the average interferogram range frequency. 7. Form the interferogram. The interferometric image is formed by multiplying one image, pixel for pixel, with the complex conjugate of the other registered image: I n k k tJ'Ki."-^*), (B.2) 8. "Flatten" the interferogram. Flattening the interferograms corresponds to removing phase trends which come from the increase in slant range across the range swath of the image (see eqn. (2.20)). The phase trend is removed by multiplying the interferogram with a complex phase function. The term flattening comes from the use of a model based on locally flat terrain used to derive the phase multiply function. The precise form for the flattening function depends on whether a reasonable estimate of the baseline is available or not. If the baseline is not available, demodulating the interferogram using the largest frequency component will perform the same task. Flattening the interferogram is useful for two reasons. Firstly, the "flattened" phase approximately represents the difference between the actual topography 171 and the reference level used for the flattening function. The residual phase after the flattening function therefore looks very much like a topographic contour plot. Secondly, flattening moves the interferogram phase derivatives to lower values which helps to minimize phase unwrapping problems. See Appendix C and Section 2.6 for more details. 9. Calculate coherence magnitude. The maximum likelihood estimate of coherence magnitude is calculated by coherently summing the interferogram values and normalizing by the geometric mean of the average squared image amplitude. Care must be taken to demodulate the interferogram data to remove underlying phase trends which tend to bias the estimate of coherence magnitude low [29]. In addition, the coherence magnitude is biased in general. Correction for this bias may be calculated by spatially averaging the estimated coherence [60]. 10. Filter/subsample the interferogram. Filtering and/or subsampling are usually performed on the raw interferogram to: • eliminate phase trends associated with layover data by halfband filtering [1], • increase the SNR of the interferogram by filtering in azimuth, • reduce the amount of data for further processing stages. The interferogram phase is a noisy estimate of the slant range travel time phase difference. Phase smoothing or some estimation process is required to reduce the effects of the noise. The usual processing scheme is coherent averaging of the interferometric image amplitudes over some local region of the image: 172 where: arg(z) = ij) = N = the phase of the complex number z, the smoothed interferogram phase, the number of image samples averaged. This is a maximum likelihood estimate of interferogram phase assuming that the interferogram phase is constant over the area of summation [61, 29]. If local frequency modulation is present, band-pass filtering centered at the estimated local phase frequency is required. 11. Unwrap the interferometric phase. One of the most challenging parts of interferometric SAR processing for topographic applications is phase unwrapping. This process is extremely difficult for two-pass satellite InSAR and is discussed in detail in Appendix C. Phase unwrapping starts from the assumption that the phase is well-behaved and well-sampled and tries to drive the inter-pixel phase derivative to lie in the range (—7r,+7r] by using the following decision rule: = tffc-i + (ipk ~ V'ifc-i) + 2?r if = ip - ib -i < -7T, k k + (ipk - tf>k-l) if X>lh k i) -l k > -7T. (B.4) In actual InSAR processing, blindly applying this rule will not result in reasonable unwrapped phase for most two-pass interferograms. See Appendix C. 12. Calibrate the interferometer geometry. The baseline parameters B and O are needed to convert the unwrapped phase to height. See Section 2.6 for more details. 173 13. Convert unwrapped phase to height. Conversion of unwrapped phase to height is a simple numerical calculation using trigonometry (see eqn. (2.14)). 14. Resample SAR image D E M to mapping coordinates. The InSAR D E M natural coordinate system is the slant-range/along-track coordinate system of the SAR. For the D E M to be useful to end-users, the D E M samples must be placed into a well-defined reference coordinate system including Earth spheroid and geoid models. 174 A p p e n d i x C Phase Unwrapping Transformation of the smoothed interferogram phase to terrain height estimates is not straightforward because of the requirement for phase unwrapping or recovering the multiples of 2TT that are lost in the arg(z) calculation. In the absence of noise and phase aliasing, one could unwrap the interferogram phase by adding or subtracting 2ir as needed to keep the phase derivative between adjacent phase samples less than 7T rads/sample. However, there are two complications: 1. phase unwrapping for SAR interferometry is fundamentally a two dimensional process, 2. noise and abrupt terrain changes such as cliffs can introduce discontinuities of magnitude greater than n rads into the interferometric phase. One expects that a closed contour integral of phase differences along any valid path in a properly unwrapped interferogram phase image should result in a net value of zero residual phase. Zero residual phase from a closed path integral corresponds to returning to the same location and hence terrain height. This is a natural constraint on the two-dimensional phase unwrapping problem for SAR 175 Al = 9/10 7i - 7/10 TC = 2/10 TC 7/10 TC A4 = 7/10 TC - -2/10 TC = 9/10 TC4 -2/10 TC 1 2TC 3 9/10 TC 2 -2/10 A2 = -2/10 TC - 9/10 TC + 2 TC = 9/10 TC TC A3 =-2/10 .n- -2/10 7t = 0 Figure C . l : Example of residue calculation. interferometry that ensures that the same unwrapped phase (to within a constant) is calculated no matter what path is taken during the phase unwrapping process. However discontinuities in the image phase can yield anomalous paths for which the residual integrated phase difference will not be zero. The existence of anomalous paths in the interferogram phase is shown by "residues" [62]. Residues denote the ends of lines of continuous phase in the interferogram. They are identified by performing a "contour integration" of unwrapped phase derivatives around the intersection of four adjacent pixels. If there is an inconsistency in the phase present, the sum of unwrapped phase difference will not be zero but ±27r. Typically, positive residues are assigned for positive residuals while negative residues are assigned for negative residuals. A residue computation is illustrated diagrammatically in Figure (C.l) where a positive residue is detected between four adjacent pixel values of 7/10 TT, 9/10 TC, —2/10 TC, and — 2/10 TC rads respectively. One calculates the derivative between the samples enforcing the rule that the magnitude of the phase change between adjacent pixels should be'less than TC rads. This gives unwrapped phase derivatives of 2/10 TC,9/10 7r,0 and 9/10 TC rads. Summing up the unwrapped phase derivatives gives a residual phase of 27r around the four pixels or a positive residue. 176 Residues can occur because of structural discontinuities as shown in Figure (C.2) and Figure (C.3). In Figure (C.2), a simulated smooth unwrapped phase surface is shown with a discontinuity. In Figure (C.3), the wrapped phase values between ±7r are shown. A positive and negative residue pair are shown as text in the figure. Any phase unwrapping path that encircles only one of the residues will be inconsistent. 20 v • ' Figure C.2: Absolute phase surface with phase inconsistency. There are several different algorithms for dealing with phase inconsistencies : 1 cutline mapping The discontinuities in interferogram phase can be dealt with by cut-line mapping [62, 63]. Cut-lines limit the choice of paths that can be taken during the phase unwrapping process to create consistency in the unwrapped phase. One strategy used to define cutline locations is to place them between residues of opposite sign which are "close". This ensures that no isolated residues can be encircled during phase unwrapping. ' A fairly extensive review of the state of the art in InSAR phase unwrapping was published by Ghiglia and Pritt in 1998 [55]. 177 However in practice, there are many sets of equally valid cutline positions which yield consistent but possibly incorrect unwrapped phase images. Manual intervention is usually required to make a reasonable set of cut-lines for an entire image. Even with manual intervention, areas with high residue density cannot be unwrapped sensibly. These areas are masked, to be interpolated after the well behaved areas of phase have been unwrapped or defined by a differently oriented data set. Areas of unwrapped phase can become isolated by unfortunate combinations of cutlines and masked out areas [63]. Again, manual intervention is also required to offset the isolated areas of phase correctly. least squares phase unwrapping The problem of minimizing the sum of unwrapped phase derivatives and the reconstructed phase surface is referred to as least squares phase unwrapping [64] This problem may be rewritten as Poissons equation with von Neumann boundary conditions (i.e. zero values outside the 10 20 30 40 x index 50 60 70 60 Figure C.3: Wrapped phase and residues for the absolute phase image of Figure (C.2). Note the positive (+) and negative (-) residues which are signs of an inconsistency in the wrapped phase due to the discontinuity. 178 boundary) [57]. Solving Poisson's equation is a well studied numerical problem with many different techniques available. Ghiglia and Romero present a fast transform method in [57]. This is an interesting method because it avoids the troublesome question of which unwrapping path to choose as with the cutline method. However, it is well known that the least squares solution will be biased because the distribution of phase wrapping errors is not generally zero mean [65]. Practically, this means that the residual phase difference between the least-squares unwrapped phase and the original phase difference has a locally smooth bias which must be post-processed to reconstruct the true unwrapped phase. weighted least squares phase unwrapping Weighted least squares methods are an improvement on least-squares phase unwrapping methods that allow one to weight the importance of different phase derivatives in the processing scheme. Weighted least-squares phase unwrapping problems are more difficult to solve and they generally require iterative methods [57, 66]. The main problem in weighted least squares phase unwrapping is the choice of weights as the influence of the chosen weights on the unwrapped phase is difficult to define. Weighted least-squares phase unwrapping still suffers from the same bias problems as unweighted least squares (although the severity of the problem is reduced). region growing phase unwrapping The region growing technique unwraps the phase into regions from a set of seed pixels according to some criterion such as coherence magnitude[58] or local phase variance. When two regions start to overlap, they are merged if their overlap regions agree to within some confidence. The process continues until no more regions can be merged. The merging process is controlled by relaxing a criterion based on data quality; the aim being to unwrap and join high quality regions initially. This algorithm tends to be slow and somewhat sensitive to algorithm parameters [58]. 179 m i n i m u m c o s t f l o w p h a s e u n w r a p p i n g The minimum cost flow algorithm is a novel approach to phase unwrapping which maps the phase unwrapping problem to the problem of minimizing the flow through a connected network [67]. This algorithm will be used by the X-SAR team during the Shuttle Radar Topography Mission (SRTM). Minimizing the flow is equivalent to minimizing h norm with respect to unwrapped phase differences. The minimum cost flow algorithm requires the use of weight terms to guide the unwrapping process. Again the choice of applied weights is somewhat arbitrary. Another algorithm which minimizes the weighted l\ norm is Flynn's minimum discontinuity algorithm [54]. This algorithm works by minimizing the weighted phase differences between adjacent pixels. m u l t i - r e s o l u t i o n p h a s e u n w r a p p i n g Davidson and Bamler present a method for phase unwrapping which focuses on the phase derivative estimation problem [38]. Rather than making a pixel by pixel estimate of the phase derivative, they use a multiresolution approach to capture information from coarse to progressively finer resolution. They discontinue the multiresolution reconstruction of the phase derivative when the local phase derivatives' variability exceeds a certain threshold. The unwrapped phase is then reconstructed using the smoothed unwrapped phase estimates. Usually this method also requires some manual intervention to complete the phase unwrapping processing as there are generally some residual phase signals after removing the smooth estimate of the phase derivatives. Generally speaking, removing as much of the existing interferogram phase variability as possible increases the likelihood of successful phase unwrapping [40]. Because residue density varies with local interferogram frequency[41], reducing the interferogram frequency reduces the number of residues and hence the complexity of the phase unwrapping problem. The simplest approach as discussed in the InSAR 180 processing algorithm is to "flatten" the interferogram phase using a flat earth model. This type offlatteningdemodulates the interferogram globally and therefore reduces the average phase difference from having some non-zero mean to having a zero or close to zero mean. Another approach to reducing interferogram phase variance is demodulating one interferogram by an integer multiple of another interferogram's phase values [40]. This pre-processing methodology subtracts an integer multiple of one wrapped interferogram phase from the other interferogram which effectively reduces the interferogram phase variability without phase unwrapping. Although the idea is straightforward, practical application of this technique appears to require accurate orbit data in concert with accurate baseline parameter estimates. The authors mention the idea of using an existing D E M in the manner of differential interferometry before the integer multiplication process but do not discuss the algorithm for inclusion of coarse DEMs in detail. 181 A p p e n d i x D D E M Refinement A l g o r i t h m Details D.l Minimization Algorithm The basic premise of workable iterative minimization algorithms is that the function is relatively well behaved [45]. If so, then a small change in the function parameters given by the k step of magnitude a in the direction p will give a new function th k k value fk+i by Taylor series approximation with the value fk+i = fk{*k + Oi Pk) k ~ A(xfe) + cx V k T fk(pk), (D.l) where: x V/(-) = = the parameters of the model, gradient of the function at the current position, If the step reduces the value of F then it is clear that V jh Pk < 0. If one were to T perform an exact search to define the optimum value of a, then one can show [45] 182 that V / T f c + l P f c = 0, (D.2) or i. e. that the gradient vector at x^+i is orthogonal to the search direction. In Newton's method [45], the function is assumed to be locally quadratic so that one can predict the future gradient as V/ f c + i = V/fc + Dpjfe. (D.3) Since at a minima, V fk+i = 0, The exact step direction is given by Pfe = - D - V / (D.4) 1 f c ) where D is the Hessian matrix. Most non-linear minimization methods apply this equation in some form with approximations made in the definition of the Hessian (D ) or in the step magnitude (a^). In practise with non-linear function minimization, the quadratic approximation may not hold and the algorithm must be iterated to reach the minimum. A number of practical problems with numerical stability and other implementation issues have stimulated development of numerous algorithms [45]. One standard algorithm for dealing with minimization of non-linear quadratic functions when the first and second derivatives are available is the Levenburg-Marquardt (LM) algorithm (see section 14.4 of [68]). The basic premise behind the algorithm is that one can use either of two methods to determine where next to move the estimated values of the minima: • if the function is locally quadratic, than the minimum parameter values can be found as above from X in m = ^-current D Vf(x where: x = the parameters of the model, 183 c u r rent) i D = the matrix of second derivatives at the current parameters, V/(-) = gradient of the function at the current position, • if the function is not locally quadratic, the next best procedure is to a take a step in the "best" downwards direction using the gradient direction Xnerct — ^-current - constantsf(x. ) current where constant is chosen to not exhaust the downhill direction. The Levenburg-Marquardt method combines these two choices by solving a related equation at each step of the iteration: ro J2^k,iSxi = (3 , (D.5) k i=i where Sxi = the perturbation in the I th parameter, k = 1.. .m, m = the number of parameters in the model, I\/ = ^Dk,i, if k ^ /, Tfc,fc = fik = (1 + A), 2^^'fa-current) k(D.6) The value of A is varied to change between the two steps noted above. If A is small, the matrix T tends to mimic D. Conversely, if A is large, the matrix V is diagonally dominant and we get a gradient step scaled by the magnitude of the diagonal elements. The algorithm is: 1. Initialize A to some small value such as 0.001 184 2. Compute f{-),T,/3. 3. Solve eqn. (D.5) for the value of the optimal step <5x. 4. Evaluate /(x + <5x). 5. If the value of / ( x + 6"x) > /(x) then increase A by a factor of 10 and return to Step 3. 6. If the value of /(x+5x) < /(x) update x and compute /(•), T,/3. Set A = A/10 and return to Step 3. Note that if one is close enough to the minima and if the function is nearly quadratic, one can step to the minima in one iteration using Newton's method. D.2 Baseline Accuracy Requirements The most stringent constraints for baseline estimation come from topographic estimation so we shall consider the impact of baseline estimation errors on estimation of terrain height in terms of height error across the ground swath. If we consider the two dimensional problem of a range line of data, the terrain height h can be estimated as (D.7) Making the substitution of 6 = 0 - <j> where <f> is the angle estimated using the interferometric SAR data yields Substituting for the interferometric angle <f> yields: 1 r + B - r2 2 h= a? + r - 2 r a cos 2 0 - arccos 185 2 2 Vb 2 ) (D.8) (D.9) To estimate the magnitude of height errors due to baseline length and orientation errors, consider the gradient with respect to all the parameters: a, S(r), B and O: dh dh dh dh dB 9 0 9 a ' 9r2 asin( ~ 9) ( —B + cos( cj>) r) asin( r h^l-cos(<p) B 9) a - rcos(9) h 2 a sin(9) h h y / 1 - r2 cos{(f>y B . (D.10) 2 The perturbation at any point in the data can be evaluated using the error values and the gradient. For the moment we consider the errors in a and S(r) to be zero which leaves: asin(9) h ~ ~ ( —B + cos( <ft) r) AB h ^/l-cos( 4>) B 2 rasin(9) A0 h + ' [ ' ' One can choose to zero the height error at any point in the swath to eliminate one or the other of the perturbations: Q _ _ «sin (9 m i d ) + cos {-B hyjl (4> id) r id) m AB m - COS {4> id) 2 m r + m i a s m {9 d) d AQ mi B k We eliminate A 0 to get ( A0 — a sin (9 id) m V hy/l AB ^ a sin ( 9 - COS ((f)midf m i d ) AB cos (<f> i ) r j m h^Jl - COS (4> d) r id m asm mi v m \ d ^ B 2 —j2 : d • [9 id) . (D.ldJ m Substituting this relation into eqn. (D.ll) snd simplifying shows that there are two terms which contribute to the perturbation in height measurement: h A 1 ( h «sin(0)r 1 \r^l-cos(<f>) 2 r COs(<j>) Vl-cos(^)' + - m . ^ d c o s t ( f > mid j COS (<l> id) | 2 _ . m A_ 1 c o s u y ' ' K We can approximate the two differences using derivatives as h A h asin(9)r I d 1 AB \d<f> r y/1 - cos( <f>) 2 186 y (<f> mid AB 2 - (j)) A B . (D.14) Evaluating the derivatives and simplifying yields: A_,„ = J.(e-»)(-tf (sin(6-0) ) 2 D.3 3 / 2 r»)..i.(» + fi (r-cos(9-9)B)Bh Cramer-Rao Bounds for 2-D Frequency Estimates The Cramer-Rao (CR) bound gives a lower bound on the variance of any unbiased estimate [69, 70]. One can calculate the CR bounds from the Fisher information matrix J, with elements: ~d \nf(z\q) 2 Jij — E dqidqj (D.17) where: q = a vector of parameters to be estimated. (D.18) The minimum error variance or error covariance of the estimates is the corresponding element of the inverse of the Fisher information matrix: E[<?; Qj] > J»~/- (D.19) In the following, we review these calculations for interferometric SAR data. Complex valued focused SAR images of distributed target scenes are characterized as circular Gaussian random variables: variables with independent real and imaginary parts, each drawn from a Gaussian distribution with zero mean and the same variance [28, 8]. We assume the samples in each image are uncorrelated but individual returns in one image are correlated with the corresponding returns in the other image. We also assume that the same coherence magnitude, local frequencies, interferogram phase, and variance values apply to the pixels being considered in a 187 particular local area of the image. The probability density function (PDF) of the correlated returns can then be modeled as the product of the PDF of individually correlated values between the two images [28]: f(...I (x,y),I {x,y),Ci{x,y),C2{x,y)...) 1 = 2 N II fkih{x, y),Ci(»,y),h{x,y),y), I ^>M> 2> (D.20) CT where: constant interferogram phase, the coherence magnitude, 2 2 variance of the images, complex amplitude of the x, y pixel of the i th y/li(x,y)t x = image, -(Nx - l ) / 2 . . . + (Nx - l)/2, (D.21) - ( A y - l ) / 2 . . . + (JVy-l)/2, (D.22) M^l_ V i(«;»Ma,o.( J 2 C l («,»)-c (», y)+«»« a ,+«.,,, , , + » ) ' 2 (1-M ) 2 exp AO 167r cT cr|(l - fi ) 2 2 2 After some tedious algebra, it can be shown that the inverse of Fishers information matrix for this data model is 72 a-** ) n 2 NxNy(l+Nx Ny -Nx -Ny ) 2 2 2 2 n 6(i-/i ) 2 Nx Ny (Nx -1) n U u 2 n U 0 0 0 C 0 t 6(i-» ) /i2 ..7 NxNy(Ny AT— i\r.. / sr..2-l) 0 2 0 2 TT 2NxN j, J (D.23) 2 yf Therefore the Cramer-Rao lower bounds for the frequency terms in the residual interferogram are 72 (1 - fi ) > „2 li Nx Ny (1 + Nx Ny - Nx 2 w 2 2 188 2 2 - Ny )' 2 <J (D.24) D.4 ""•1 " ^ W s i V ^ W - l ) ' ( D ' 2 5 ) :•*] * ^ i v ' ^ w - D - ( D - 2 6 ) Conservation o f M e a n Values i n L S Problems The general least squares problem is min £ ( / ( * ! , . . . , z ) - z z ) . (D.28) 2 P The solution to the least squares problem is obtained by solving the set of equations: £(/,(*!,...,x ) - z,) P d ^ - - " x r ) = o i = 1,...,P. (D.29) If the <7 derivative is constant and non-zero over the / samples, clearly one can t/l make an equation: Y,(f (x ,...,x )-z ) l 1 P l = 0 (D.30) which implies that the mean of /(•) at the solution and the mean of the function to which we are fitting are equal. In the more general case, the derivatives are not likely to be constant over all / samples. However the same result can hold if some combination of the derivatives at the solution can yield a non-zero constant. In other words if }^}^{fi{xi,---,xp) i=i - zi) a / -r ° =0 (D.31) { at the solution such that P C{ ^ ^ ' —— = constant, 1 i=i d x i 189 (D.32) then the equality between the fitted function (27) and the output estimate holds. For baseline parameter minimizations in the D E M flattening or updating problems, the errors in the fitted function due to errors in B in 0 are approximately linear functions. The functions are not co-linear so they form a non-orthogonal basis for any linear or constant function. As such, they fulfill the requirements of eqn. (D.32) for maintaining the mean of the fitting function (z{) in the fitting process. D.5 Plane Fits to R a n d o m Noise One way to derive the best possible topographic accuracy achievable for an input D E M is to analyze the experiment of fitting a linear model of the D E M errors due to baseline parameter errors to a set of random, zero-mean, independent, and identically distributed random variables with variance a . 2 This models the ideal case of D E M fitting. The problem is to minimize: a b ^2(a + b y + cx + dx y - Ah) 2 = 11 A c d where: a models the constant error of the output D E M , 6 models the range slope error of the output D E M , c models the azimuth slope error of the output D E M , d models the bilinear trend error of the output D E M , Ah A models the errors of the input coarse D E M , = the matrix of fitting coefficients, 190 (D.33) h = the D E M errors in a vector. and V x (Ny - 1) = ,k = 0...N i_ = \l = l (D.34) - 1, y (D.35) Q...N -l. x The solution to the minimization problem is the usual least squares solution: a b = (A A) T A h. (D.36) T 1 c d Since we have taken care to choose ourfittingcoefficients to be orthogonal, ^ A A ^-l T takes a simple diagonal form: 12 N N (N-£-l) (A-A)- = 1 x y 12 N N,(N2-l) y AT| N*-NZ 144 N -N* y N +N x x N y J (D.37) The values of the cofficients are therefore: A = ]TVE M> A/I , 12 V - 2 d = Ny-1 , j ; ^ ! - ^ ) , 1 2 N N (N -l)^ y . x 144 o A , ,, Ny-1.. N - l iVj N* — N% N — N N + N N 3 y x x y k,i (D.38) 191 If the noise variables are independent and zero mean, then it is clear that the expected value of all the fitted variables is 0, ie E [a] = 0, E [b] = 0, E [c] = 0, E[d] = 0. A corollary of this result is that the output D E M will not (in the mean sense) have any trends. In addition, because the basis functions that are being fitted to are orthogonal, there is no correlation between the parameters, i.e. E [ab] = 0. The variance of the parameters may be calculated by taking the expectation of the above expressions. After some tedious algebra it can be shown that: E N Ny- L X 12 a 2 E[& I = 2 N N (N*-iy x y 12 a N N {N -\y 2 E 2 x y 144 cr N -N%N 2 N 2 JV3 - 7Y3 x + N Ny y x The variance of the output plane also establishes bounds on what the best fitting interferometric D E M could achieve. Because the fitted plane has zero mean, one must evaluate N~Ny~ [ E E ( ^ E + " + ( l ^) (D.39) ; Again after some algebra, one can show that 3a 2 L - E p r y ^ y + ca + dsy) ' 2 (D.40) N Ny N Ny X X The variance of the height difference across the data set is of interest because it gives an idea of the possible output trends in the InSAR D E M . One can calculate this as E[(h(x ,y 2 | a,b,c,d) - h(x yi \ a,b,c,d)) ] = 2 2 u [a + b y + c x + d x y - [a + b y + c x + d x y )}' 2 2 2 2 2 192 2 2 2 (D.41) Simplifying this expression by substituting for the differences yields [b (j/2 - E ^[bAy + cAx + d Axy] + c (x - zi) + d (x y - xi yi)f 2 2 2 2 , (D.42) where Ay = y - yi, Ax = x - xi, and Axy = x y 2 2 2 x\ y\. Expanding and using 2 the orthogonality of the basis functions yields [6 Ay + c Ax + dAxyf + E Ay 2 Ax 2 + E d? Axy . (D.43) 2 Substituting for the expectation operators yields: E [(h(x , y I a, b, c, d) - h(x yi 2 2 u 12 CT \ a, b, c, d)) } = 2 o 2 N N (N -l) x y +m A 12<r 2 * 2 Ax' N N ( A - 1) 2 y 144 a x 2 (D.44) •Axy N?, — N% N — A, N + N N, 2 3 v x x The maximum error possible occurs as a function of the relative position of the two points examined. There are six possible directions that could have the maximum values, (four directions from one corner to the next along an edge, and two differences diagonally). However one need only check diagonal direction and one edge to edge direction to get the form of the error. For the difference across the fitted D E M diagonally we have: Ax = Ay = Axy = Nx 1, Ny-l, 0, and (h(x , y I a, b, c, d) - h(x yi 2 2 u | a, 6, c, d))' _ diag ' 12 a (N - 1) 12cr (N - 1) 2 2 y N Ny(Ny + 1 ) X x + Ny N (N X X + l) " (D.45) 193 For the D E M difference across an edge we will evaluate the along the right edge of the D E M where we have: Ax 0, = Ay Ny-l, Axy A^-l = A _ r 1 2 _ (N -l)-(N -l) x y (N -1) = x (N -1) y { B A 6 ) and E [(h(x ,y 2 | a,b,c,d) - h(x ,y 2 l \ l a,b,c,d)) ] 2 right 144 a 12 a (N - 1) 2 2 y N N (N x y y + 1) + ' NiJV3 - Ni N - JV3 iV* + TV, N y This simplifies to E [(h(x , y I a, b, c, d) - h(x 2 2 For JVj, = N x y y right ' = 2 y x x 36a (iV - 1) (A^ - 1) 2 N N (N y | a, b, c, d)) ] 2 lt 12 a (N - 1) 4 ^ " y (D.48) y + 1) A^iV,, ( A ^ + JV + iV* + 1)" y y = iV, the edge difference will be more than the difference diagonally which suggest that the bounds for error performance of the estimate are really set along the edge of the D E M as _ 24a (2N-1) (N-l) 2 (h(x , y | a, b, c, d) - h(x 2 2 u y | a, b, c, d))' 1 J right (N + 1) N 2 2 (D.49) 194 A p p e n d i x E D E M Refinement A l g o r i t h m Definition The required inputs for the D E M refinement algorithm are • I e ~ ( ( ' ')) the raw interferogram (i.e. unflattened) at the range time r and J 5 r r ? ) the azimuth time 77, • /J(T, 77), the coherence magnitude at each interferogram sample, • s(rj), the estimated spacecraft position in meters at each azimuth time in ECR coordinates, • r(r), the slant range in meters at each range sample, • A, the radar wavelength in meters, • h(x,y,z), the input D E M positions in meters in ECR coordinates, 195 DEM Flattening DEM Updating Phase Unwrapping (Prepare DEM Data for USE / Register DEM I To Interferogram Generate Unwrapping Mask { j Estimate Baseline Parameters and Compute Residual Interferogram Filtered Residual Interferogram Phase / Filter Residual I Interferogram Unwrapped Residual Interferogram Phase Generate Group Labels V Re-estimate Baseline Parameters and 1 InSAR DEM ' Modify PhaseN i ^ Group Offsets y i y Unwrap Residual Interferogram Refined InSAR DEM Figure E . l : Top level block diagram for the D E M refinement algorithm. The algorithm proceeds in three stages: 1. D E M flattening, 2. phase unwrapping, and 3. D E M updating, as described in Chapter 3 and shown in Figure (E.l). The following list gives a concise description of the processing at each stage. • D E M Flattening Algorithm 1. Prepare the D E M data (h(x,y,z)) for use: (a) Convert the D E M sample positions to preliminary along-track time coordinates by finding the satellite position that minimizes the distance to D E M sample position. 196 (b) Grid the height data along lines of constant azimuth time and at azimuth increments that conserve the information in the input D E M data. (c) Search along each constant azimuth line to identify and discard layover points. Layover points are D E M samples where the derivative of the estimated slant range becomes negative. (d) Produce the input coarse D E M for processing, h(DEM), by gridding the D E M data at constant azimuth positions at the same azimuth time and slant range increments as the input interferogram data. 2. Register the D E M data to the interferogram by maximizing coherence magnitude over small regions of the input interferogram with distinctive features. The coherence is maximized as a function of the baseline parameters using the D E M flattening algorithm described in the next step. 3. Estimate the baseline parameters for the entire scene and compute residual interferogram. (a) Calculate/set the following values for the input coarse D E M in slant range coordinates: J(r), the target slant range difference for the fitting problem; B, the initial estimate of the baseline magnitude; 0#, the initial estimate of the baseline orientation; (b) while not converged and maximum number of iterations unreached do: i. Estimate B and 0 g by solving: min(fl, 9) £ > f ) k f c - (r2 (B, 0, a, h, r) - ii. Form the residual interferogram: 197 k r ))\ k iii. Estimate the frequency components: w , w ,w r a and constant ria phase offset ip from the residual interferogram, 0 iv. If the frequency components are all close to zero then the iteration has converged otherwise update S(r) as: 5{?) = (r2(B,Q,a,h,r) - r) + ^ w (T-T ) R (w , (n - r] )(T - T ) + r a + w (r)-r] ) 0 a 0 0 0 + Vo)- and return to Step i. 4. Filter the residual interferogram to prepare for phase unwrapping. Phase Unwrapping 1. Generate a phase unwrapping mask to eliminate low coherence, low backscatter, and high magnitude areas (i.e. layover areas) from the phase unwrapping processing. 2. Generate group labels (groups) for the unwrapped interferogram phase. groups is generated so that unwrapped phase values with the same group number can be identified. 3. Unwrap the residual interferogram phase. (a) Unwrap the bulk of the residual interferogram phase using weighted least-squares processing [57]. (b) Unwrap the residual from the weight least-squares phase unwrapping using region growing [71] or minimum discontinuity processing [54]. D E M Updating Algorithm 1. Generate an estimate of the total slant range difference using the unwrapped phase and the calulated 5(r) from the flattening algorithm as: ^ = %) + ! _ * . 198 (E.l) where: \I> = the unwrapped phase . 2. Redo the baseline parameter estimation by minimizing: min(B,0) M > I <*F),a,r)) . s J2( ( h D E M ) ~ 0 2 (E-2) where: h(B,Q h(DEM) = | S(r),a,r) = the input coarse DEM, the D E M calculated from the unwrapped phase. 3. Calculate terrain height for all processed data using the unwrapped flattened interferogram phase and the estimated baseline parameters from the minimization above. 4. Validate the phase unwrapping ambiguity for all contiguous groups of unwrapped phase. For each group of data (a) Compute the mean difference between the output height and the input D E M over the same pixels. (b) If the mean difference of a unwrapped phase group corresponds to greater than the half-ambibuity height for the interferometer, add/subtract the appropriate number of half wavelengths to/from the slant range difference estimate and return to step 2 of the D E M updating algorithm. 199 Appendix F ERS DEM Histograms 200 (a) Height differences between input TRIM and output InSAR-TRIM DEMs. -250 -200 -150 -100 -50 0 50 Height difference in meters 100 150 200 250 (b) Histogram of InSAR-TRIM and TRIM D E M differences. Figure F . l : InSAR-TRIM D E M error analysis for ERS Tandem data. 201 Height difference in meters (b) Histogram of input DTED-1 and input TRIM D E M differences. Figure F.2: Input DTED-1 D E M error analysis. 202 (a) Height differences between output DTED-1 and input TRIM DEMs. -250 -200 -150 -100 -50 0 50 Height difference in meters 100 150 200 250 (b) Histogram of output InSAR-DTED-1 D E M and input TRIM D E M differences. Figure F.3: Output InSAR-DTED-1 D E M error analysis for ERS Tandem Mission data. 203 (a) Height differences between modified DTED-1 D E M and input TRIM D E M . x10 10 j 9 4 1 - -250 -200 1 1 '• s -150 -100 1 1 1 1 \ • -50 0 50 Height difference in meters 100 1 1 ; 150 200 - 250 (b) Histogram of modified DTED-1 D E M and input TRIM D E M differences. Figure F.4: Modified DTED-1 D E M error analysis. 204 (a) Height differences between modified InSAR-DTED-1 D E M and input TRIM D E M . x 10 4 meai i : - 1 . 2 5 n l o : 26.9 m 9 0 * % of -250 -200 -150 -100 -50 0 50 Height difference in meters 100 I A h : : 43.2 m 150 200 250 (b) Histogram of modified InSAR-DTED-1 D E M and input TRIM D E M differences. Figure F.5: InSAR-Modified DTED-1 D E M error analysis. 205 (a) Height differences between input GTOPO30 and input TRIM DEMs. x 10* 6m -250 -200 -150 -100 -50 0 50 Height difference in meters 100 150 200 250 (b) Histogram of input GTOPO30 and input TRIM D E M differences. Figure F.6: Input GTOPO30 D E M error analysis. 206 (a) Height differences between output GTOPO30 and output InSAR-TRIM DEMs. -250 -200 -150 -100 -50 0 50 Height difference in meters 100 150 200 250 (b) Histogram of output InSAR-GTOPO30 D E M and input TRIM D E M differences. Figure F.7: Output InSAR-GTOPO30 D E M error analysis. 207 (a) Height differences between modified GTOPO30 and output InSAR-TRIM DEMs. -250 -200 -150 -100 -50 0 50 Height difference in meters 100 150 200 250 (b) Histogram of modified GTOPO30 D E M and input TRIM D E M differences. Figure F.8: Modified GTOPO30 D E M error analysis. 208 (a) Height differences between modified GTOPO30 and output InSAR-TPJM DEMs. -250 -200 -150 -100 -50 0 50 Height difference in meters 100 150 200 250 (b) Histogram of modified InSAR-GTOPO30 D E M and input TRIM D E M differences. Figure F.9: InSAR-Modified GTOPO30 D E M error analysis. 209 Appendix G R A D A R S A T D E M Histograms 210 (a) Height differences between input TRIM and output InSAR-TRIM DEMs. 18 x 10 16 14 h ean : -0.39 m 1 12 5 % of I A h ' : 18.6 m 10 6 m 1 1 L h -150 : 12.2 m -100 -50 Height difference in meters 50 100 150 (b) Histogram of InSAR-TRIM and TRIM D E M differences. Figure G . l : Output InSAR-TRIM D E M error analysis for RADARSAT-1 data. 211 (a) Height differences between input DTED-1 and input TRIM DEMs. x 10 4 I r 6 rrean : -16.44 m o : 35.6 m 4 9 3 %df I Ah | :61.3m m 2 0 -150 -100 -50 0 50 Height difference in meters 100 150 (b) Histogram of input DTED-1 and input TRIM D E M differences. Figure G.2: Input DTED-1 D E M error analysis. 212 (a) Height differences between output DTED-1 and input TRIM DEMs. x10 161 4 1 1 1 1 — 1 1 14h -150 -100 -50 0 50 Height difference in meters 100 150 (b) Histogram of output InSAR-DTED-1 D E M and input TRIM D E M differences. Figure G.3: Output InSAR-DTED-1 D E M error analysis for RADARSAT-1 data. 213 (a) Height differences between input modified DTED-1 D E M and input TRIM D E M . rr ean : -0.06 m o : 34.5 m 9 3*% of | "A"li 1':'56.9 rri 2 0 -150 -100 -50 0 50 Height difference in meters 100 150 (b) Histogram of input modified DTED-1 D E M and input TRIM D E M differences. Figure G.4: Modified DTED-1 input D E M error analysis. 214 -150 -100 -50 0 50 Height difference in meters 100 150 (b) Histogram of modified InSAR-DTED-1 D E M and input TRIM D E M differences. Figure G.5: Output InSAR-Modified DTED-1 D E M error analysis for RADARSAT1 data. 215 (a) Height differences between input GTOPO30 and input TRIM DEMs. x 10 4 rr ean : 7.83 m o : 98.3 m 9 3*% of -150 -100 -50 0 50 Height difference in meters I Ah I 43.7 : m 100 150 (b) Histogram of input GTOPO30 and input TRIM D E M differences. Figure G.6: Input GTOPO30 D E M error analysis. 216 (a) Height differences between output GTOPO30 and output InSAR-TRIM DEMs. -150 -100 -50 0 50 Height difference in meters 100 150 (b) Histogram of output InSAR-GTOPO30 D E M and input TRIM D E M differences. Figure G.7: Output InSAR-GTOPO30 D E M error analysis for RADARSAT-1 data. 217 (a) Height differences between input modified GTOPO30 and input TRIM DEMs. x 10 4 10 mean : 0.08 m a':: 91.1 m 9 0 * % of I Ah I 34.2 : m 5h -150 -100 -50 0 50 Height difference in meters 150 (b) Histogram of modified input GTOPO30 D E M and input TRIM D E M differences. Figure G.8: Modified input GTOPO30 D E M error analysis. 218 (a) Height differences between modified GTOPO30 and output InSAR-TRIM DEMs. 15 10 5 0 -150 -100 -50 0 50 Height difference in meters 100 150 (b) Histogram of modified InSAR-DTED-1 D E M and input TRIM D E M differences. Figure G.9: Output InSAR-Modified RADARSAT-1 data. 219 GTOPO30 D E M error analysis for
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Refining low-quality digital elevation models using...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Refining low-quality digital elevation models using synthetic aperture radar interferometry Seymour, Michael S. 1999
pdf
Page Metadata
Item Metadata
Title | Refining low-quality digital elevation models using synthetic aperture radar interferometry |
Creator |
Seymour, Michael S. |
Date Issued | 1999 |
Description | Two-pass synthetic aperture radar (SAR) interferometry (InSAR) is a technique for processing the phase difference between coincident SAR images to obtain the range difference from the two radars to a common point on the earth's surface. The accuracy of the range difference measurement is in the order of one millimeter, and this range information can be processed to obtain digital elevation models (DEMs) of the surface topography. The digital processing required to make the DEM is quite complicated, mainly due to two factors. Firstly, the phase information is obtained from complexvalued data and therefore lies between -π and +π whereas the complete phase information is needed. To obtain this, the phase must be "unwrapped" where the missing integer number of 2π are estimated for each data sample. Secondly, the geometry of the satellite passes relative to each other must be known to an accuracy of a few millimeters in order to obtain the surface height values to the required accuracy (about 10 m). Both of these steps require supplemental information and manual guidance to be performed correctly. Phase unwrapping is difficult because of noise and undersampling inherent in the measurements. The geometry estimates are difficult to make because the orbit is only known to an accuracy of a few meters, and the received phase data is a non-linear function of the satellite geometry. In the past, the geometry estimates have been made using known ground control points (GCPs), which requires a considerable manual effort, and has its own set of errors. The objective of this thesis is to use supplemental information in the form of a coarse DEM to make the InSAR processing more accurate and more automatic. We achieve this objective by developing a new algorithm which incorporates the coarse DEM directly into the processing stream, with the result that phase unwrapping and geometry estimation are performed accurately and reliably. In effect, the input DEM points serve as a large, dense set of GCPs. While the accuracy of each input DEM point is not very high, the large number of them provide adequate geometric accuracy, particularly as an automatic algorithm can register them directly to the radar data. There are two key steps in the new algorithm, which are interwoven in an iterative framework. First of all, the satellite geometry is estimated from the DEM and interferometric phase. This is done with a non-linear, iterative optimization algorithm without having to unwrap the phase. Avoiding phase unwrapping is important, as phase unwrapping errors can significantly bias the geometry estimates. Second, the input DEM along with the refined satellite geometry are used to create a model of the unwrapped interferogram phase that should be received from the two satellite passes. When this phase is wrapped, and compared with the measured phase, a differential interferogram is obtained which represents the difference between the coarse input DEM and the topography as measured by the satellite. This differential interferogram has a relatively low bandwidth, which means that it can be filtered and unwrapped reliably and accurately. Finally, the information in the unwrapped interferogram is used to refine the grid spacing and vertical accuracy of the coarse DEM. We have used mathematical analysis and simulation to develop the algorithm, to obtain statistical quality measures and to understand what system parameters affect the accuracy of the DEM results. We find that the main factors affecting accuracy are the interferometer's sensitivity of phase to height and the number of available DEM points, including the size and variability of the input DEMs' errors. We have successfully applied the DEM refinement algorithm to ERS Tandem Mission and RADARS AT-1 data. The generated InSAR DEMs had standard deviations of 12 to 20 meters compared to a control DEM with approximately 3 meters standard deviation. The output InSAR-enhanced DEMs had two to four times improvement in height accuracy compared with the input DEMs. In this way we have demonstrated that one can generate reliable estimates of topography for standard SAR scenes without having access to precision orbit data. Thus we have shown that the processing bottlenecks in dealing with repeatpass satellite InSAR data can be overcome, and useful topographic information can be obtained from the vast supply of existing InSAR data sets. |
Extent | 33487019 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-07-03 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065269 |
URI | http://hdl.handle.net/2429/10094 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1999-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1999-464210.pdf [ 31.94MB ]
- Metadata
- JSON: 831-1.0065269.json
- JSON-LD: 831-1.0065269-ld.json
- RDF/XML (Pretty): 831-1.0065269-rdf.xml
- RDF/JSON: 831-1.0065269-rdf.json
- Turtle: 831-1.0065269-turtle.txt
- N-Triples: 831-1.0065269-rdf-ntriples.txt
- Original Record: 831-1.0065269-source.json
- Full Text
- 831-1.0065269-fulltext.txt
- Citation
- 831-1.0065269.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065269/manifest