Estimating the 3-D Flow Rates of the Lowell Glacier with Using Spaceborne InSAR by Xiangzhou Joe Zhang B . A . S c , Tsinghua University, Beijing, China, 1990 M . A . S c , Chinese Academy of Sciences, 1993 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF Masters of Applied Science in T H E F A C U L T Y OF G R A D U A T E STUDIES (Department of Electrical and Computer Engineering) We accept this thesis as conforming to the required standard The University of British Columbia April 1998 © Xiangzhou Joe Zhang, 1998 In presenting t h i s thesis in p a r t i a l f u l f i l l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y available for reference and study. I further agree that permission for extensive copying of t h i s thesis for s c h o l a r l y purposes may be granted by the head of my department or by his or her representatives. It i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of F><< if,cv! E,U(.A etr't'a The University of B r i t i s h Columbia Vancouver. Canada Abstract Glacier motion data is very important for glaciology research, but traditional methods for collecting such data need to overcome staggering financial and logistical hurdles. During the period from July 1995 to May 1996, ERS-1/2 operated under Tandem Mission Mode with one satellite following the other one day apart. This gives a great opportunity to using Synthetic Aperture Radar (SAR) interferometry (InS AR) technology to measure the relative large displacement on the ground between the observation pair. Pioneer research work shows that this short time interval can maintain relatively high coherence for most glacier movement investigations. Some good results have been achieved to measure the velocity field of ice fields, ice streams and Alpine glaciers [1, 2, 3, 4]. In this study, the rationale and processing scheme of differential InSAR are discussed and examined first. Then 10 Tandem raw data sets (4 descending passes, 6 ascending passes) over the Lowell glacier at the boundary of B . C . and the Yukon, Canada, have been processed. 6 interferograms with enough coherence magnitude (3 ascending passes, 3 descending passes) have been generated. In concert with three different flow assumptions on the glacier movement direction, the line of sight (LOS) displacements along the glacier centreline measured from the ascending orbit and the descending orbit are converted into the 3-D velocity vectors. Our results show that the surface parallel assumption is more suitable for most part of the glacier. An overall accuracy around 4 cm/day rms has been achieved in measuring the surface ii motion of the Lowell Glacier. However, the 3-D projection model does not apply for the glaciers at all loca-tions on this planet. The glacier has to be away from the critical regions caused by the projection geometry. Also, in order to achieve high estimation accuracy, the flow direction of the candidate glacier has be approximately aligned with the radar LOS. But with specifical satellite missions, this technique can provide a feasible method for monitoring the global glaciers and ice sheets in the future. in C o n t e n t s A b s t r a c t i i Contents iv Lis t of Tables v i i i L i s t of Figures ix L i s t o f Acronyms xii i Acknowledgements x v 1 Introduct ion 1 1.1 Background 1 1.2 Objectives of Research 5 1.3 Outline of the Thesis 6 2 Rev iew of I n S A R Theory 8 iv 2.1 Introduction 8 2.2 Fundamentals of InSAR Theory 9 2.2.1 Basic principles of Synthetic Aperture Radar 9 2.2.2 InSAR geometry and height equations 14 2.2.3 The relationship of interferometric phase with baseline and to-pography 17 2.2.4 Differential interferometry 19 2.3 Correlation in Interferometric Images 23 2.3.1 Temporal correlation coefficient 26 2.3.2 Baseline correlation coefficient 27 2.3.3 Volume correlation coefficient 29. 2.3.4 Thermal correlation 30 2.3.5 Processing correlation 31 2.4 Multilook Technique 36 3 I n S A R Processing Algor i thms for Surface Displacement Measure-ment 40 3.1 Introduction 40 3.2 Registration - 43 v 3.2.1 Registration offset estimation 44 3.2.2 Interpolation 49 3.3 Flattening 51 3.4 Geocoding and Removing Topography Phase 53 3.5 Multilook Averaging 54 3.6 Phase Unwrapping 55 3.7 Converting Phase Information into 3-D Ground Displacement Vectors 57 3.8 Baseline Estimation 58 3.8.1 Relationship of topography and velocity accuracy with the base-line error 58 3.8.2 Baseline estimation 60 4 Glac ie r M o t i o n Es t ima t ion 63 4.1 Introduction 63 4.2 Study Area 65 4.3 ERS-1/2 Interferograms of the Lowell Glacier 69 4.4 Estimation of the LOS Displacement 72 4.5 Projection to Glacier Flow Direction 76 vi 4.5.1 Estimating the 3-D glacier displacement using single-track LOS measurements 81 4.5.2 Estimating the 3-D glacier displacement using dual-track LOS measurements 86 4.6 Discussion 92 4.7 Conclusions 95 B i b l i o g r a p h y 98 vii L i s t o f T a b l e s 2.1 Radar Parameters of ERS-1 SAR 10 4.1 Summary of ERS-1/2 passes over the Lowell Glacier 69 viii L i s t o f F i g u r e s 2.1 The geometry model of SAR 11 2.2 The geometry model of InSAR 15 2.3 Differential interferometry InSAR geometry 21 2.4 Baseline correlation coefficient as a function of Bn with typical ERS parameters 28 2.5 Correlation coefficient as a function of penetration depth L and normal baseline Bn 30 2.6 The thermal correlation coefficient, pthermai >as a function of SNR . . 31 2.7 Correlation coefficient, pr-registration, as a function of range direction registration error 33 2.8 Correlation coefficient, px-registration, as a function of azimuth direction registration error 34 ix 2.9 Correlation coefficient, pDe-mismatch, as a function of Doppler centroid mismatch 35 2.10 Phase error distribution for different number of looks with \p\ = 0.7 and 0 = 0 38 2.11 Phase standard deviation versus \p\ for different number of looks with 6 = 0 39 3.1 Block diagram of the InSAR processing steps to realize surface dis-placement measurement 41 3.2 The relationship of the cross-correlation peak with the slope of the phase ramp 46 3.3 Relationship of elevation error with parallel baseline component error. 60 3.4 Relationship of elevation error with normal baseline component error. 61 4.1 Latitude/Longitude of the related ERS-1/2 SAR frames 67 4.2 SAR magnitude image, Oct. 23,1995, showing an overview of the Lowell Glacier. Azimuth: top-to-bottom; Ground range: right-to-left. The scene covers an area of 23 km x 60 km 68 4.3 Canadian geographic map showing an overview of the Lowell Glacier. 68 4.4 Representative descending pass interferogram from October 22/23, 1995. 71 4.5 Representative ascending pass interferogram from January 12/13, 1996. 72 x 4.6 Elevation along the glacier centerline 74 4.7 Lowell Glacier LOS flow rate from descending orbit data 76 4.8 Lowell Glacier LOS flow rate from ascending orbit data 77 4.9 The geometry projecting LOS displacement to surface displacement. . 78 4.10 Plot of the critical angles in the calculation of the surface displacement for the descending and ascending passes of ERS-1/2 80 4.11 Velocity vectors in an idealized glacier 82 4.12 The vertical flow direction of the Lowell Glacier based on the surface-parallel and moraine-aligned assumptions 84 4.13 The horizontal flow direction of the Lowell Glacier based on the surface-paralle and moraine-aligned assumptions 85 4.14 The Lowell Glacier surface velocity along the glacier centreline us-ing single-track LOS measurements, under the surface-parallel and moraine-aligned assumptions 86 4.15 The vertical flow direction of the Lowell Glacier based on the surface-parallel and greatest-slope assumptions 87 4.16 The horizontal flow direction of the Lowell Glacier based on the surface-parallel and greatest-slope assumptions 88 xi 4.17 The Lowell Glacier surface velocity along the glacier centreline us-ing single-track LOS measurements, under the surface-parallel and greatest-slope assumptions 89 4.18 The horizontal flow direction of the Lowell Glacier derived from dual-track LOS measurements 90 4.19 The vertical flow direction of the Lowell Glacier derived from dual-track LOS measurements 90 4.20 Lowell Glacier 3-D flow rate derived from dual-track LOS measurements. 91 xii L i s t o f A c r o n y m s ATI Along-track Interferometry CCRS Canadian Centre for Remote Sensing dtSAR Desk-top SAR processor D E M Digital Elevation Model D P C A Displaced Phase Center Antenna D-PAF Deutsche Processing and Archiving Facility ERS -1 /2 The first/second European Remote Sensing Satellite ESA European Space Agency F F T Fast Fourier Transform F M Frequency Modulation GPS Global Positioning System IFFT Inverse F F T InSAR SAR Interferometry J P L Jet Propulsion Laboratory LOS Line of Sight N A S A National Aeronautics and Space Adminstration NSERC National Sciences and Engineering Research Council P D F Probability Density Function RRSG Radar Remote Sensing Group SAR Synthetic Aperture Radar SLC Single Look Complex product SNR Signal-to-Noise Ratio xiii S P E C A N Spectrum Analysis algorithm U K - P A F United Kingdom Processing and Archiving Facility U T M Universal Transverse Mercator X T I Cross-track Interferometry xiv A c k n o w l e d g e m e n t s I would like to begin by thanking my adviser Dr. I. G. Cumming. His inspiration, guidance, and support not only helped me with this research but also have prepared me well for what lies ahead. I also want to thank my colleagues for making Radar Remote Sensing Group (RRSG) a friendly and interesting place to spend my years at University of British Columbia. Special thanks to Michael Seymour and Wei X u discussions with them led to many of the insights needed to carry out this research. I also would like to thank MacDonald Dettwiler, the National Sciences and Engineering Research Council (NSERC) of Canada, the B C Advanced Systems In-stitute and the B C Science Council for funding this research. I am grateful to the European Space Agency (ESA) for providing ERS data. Finally, I want to thank my wife, L i , for her love and consistent support during my hard time. X I A N G Z H O U Z H A N G X V C h a p t e r 1 I n t r o d u c t i o n 1.1 B a c k g r o u n d Glaciers and ice-sheets form the largest component of perennial ice on the earth. Over 75% percent of world's fresh water is presently locked up in these frozen reservoirs [5]. The fluctuation of glaciers and ice-sheets from glacial periods to interglacial periods has been one of the most dramatic climate signals in earth's history. Recent evidence of the rapidity with which sea level has changed due to discharges of ice-sheets and glaciers has heightened our awareness of the dynamic nature of this icy component of our climate. Whi le the mass status of two of the largest ice-sheets in Greenland and Antarc-tica is s t i l l uncertain, the contribution of the mountain glaciers to the sea level rise is well known. The current rate of sea-level rise is estimated to be about 1.75mra/yr, 1 with large contribution from the melting of mountain glaciers (around 0.5mm/yr) [6]. In addition, mountain glaciers are important components of the hydrological cycle in local areas. So it is very important to determine whether the total' glacier mass is stable, shrinking or growing. High-resolution elevation data over the glaciers obtained through investigating techniques such as photoclinometry and field survey can be used to estimate glacier mass balance. However there are several weaknesses in these data. First, the num-ber of mountain glaciers monitored is too small to be representative, especially in remote areas. Secondly, field surveys are very costly and are usually limited by the weather. Finally, much of the data for glacier mass change actually are derived from length change other than changes in elevation which is more directly related to mass balance [7]. Interferometric Synthetic Aperture Radar (InSAR) can possibly provide infor-mation to remedy all of the above weaknesses. Since Graham [8] first demonstrated InSAR by using an airborne X T I SAR system in 1974, a lot of interesting research work has been done to examine the potential usefulness of InSAR systems. Major progress has been achieved in making Digital Elevation Models (DEMs) of the earth's surface by measuring the phase differences between two registered SAR images, taken with a shift in the viewing angle of the sensor [9, 10, 11, 12, 13]. A n accuracy of 10 m rms can be achieved for topography mapping over a quite large and relatively flat, stable area [4]. When the images are acquired separately from near-repeat orbits, the sources of interferometric phase are ground topography and surface displacement during this 2 time interval. By subtracting the phase component contributed by the topography from the repeat pass interferogram phase, we can estimate the surface displacement. Gabriel and Goldstein et al. [14] proposed a technique to map small elevation changes over large areas by using differential InSAR in 1989. Three SEASAT observations of an area in the Imperial Valley, California were used to make two interferograms and then a double-difference interferogram. Small motions (2 — 3 cm) of the fields due to swelling or shrinking associated with watering were detected. Massonnet and Rossi et al. [15] showed the usefulness of InSAR technology to monitor the earthquake-displacement by using ERS-1 satellite data obtained before and after the earthquake in Landers, California, June 28, 1992. This geodetic tool provides a denser spatial sampling (100 m per pixel) than traditional surveying methods and a better preci-sion (around 3 cm) than previous space imaging techniques. For deeper earthquakes associated with little or no surface rupture, this technique will also be a powerful tool for measuring surface displacement in the intermediate and far fields. Up to date, observable progress has also been made in estimating the velocity field of ice-sheets and glaciers. Goldstein et al. [16] examined the application of InSAR to the Rutford Ice Stream, Antarctic. Results show that the detection limit is about 1.5 mm for vertical motions and about 4 mm for horizontal motions in the radar looking direction. Kwok et al. [17] and Joughin et ai [18] demonstrated that repeat-pass interferometry with ERS-1 images can be used to measure small ice-sheet motion over the wide area of northern Greenland. Because the ice-sheet surface is relatively flat and stable, the topography can be derived from differential interferograms formed from sequential observations. With this measurement, a pure displacement field can 3 then be obtained by removal of the topographic contribution to the interferometric phase at each pixel. But for glaciers, high relief, high spatial variability, and rapid temporal varia-tion may complicate the acquisition and interpretation of the InSAR data. The as-sumption of constant surface displacement rate normally does not hold here. Direct abstraction of the topographic information from multiple inferograms is not possible. Thus an extra reference D E M is needed to provide the topographic phase. Further-more the fast movement of the glaciers sometimes can even cause the repeat-pass SLC SAR data pairs lose coherence with a few days acquisition interval, which dramati-cally limits the interpretation accuracy or even makes it impossible to measure the glacier motion. Hence shorter observation intervals are needed. During the period from July 1995 to May 1996, ERS-1/2 operated under Tandem Mission Mode with one satellite following the other one day apart. This short time interval can maintain relatively high coherence for most InSAR applications. So far, several glaciers in Greenland [1] and Antarctica [2] have been investi-gated using ERS- l /ERS-2 Tandem Mission data. Some promising results have been achieved in understanding the dynamic nature of the glacier system. Vachon et al. [3] have also demonstrated that space-borne InSAR can be a very feasible means to monitor the Saskatchewan and Athabasca glaciers. The potential accuracy of this technique could be 4 cm/day rms for displacement measurements [4]. 4 1.2 O b j e c t i v e s o f R e s e a r c h The major goal of this research is to further investigate the capability of spaceborne InSAR techniques for estimating the 3-D flow rate of glaciers. A typical mid-latitude alpine glacier, the Lowell Glacier, with a large size but moderate slope was chosen and analyzed. Ten raw data pairs (4 descending passes, 6 ascending passes) over the Lowell glacier have been processed and analyzed. Several ascending and descending pass interferograms with acceptable coherence level have been created. The main focus of this dissertation is interferometric estimation of glacier mo-tion. This includes an inquiry into the processing issues and the source of the com-plexity observed in motion-only interferograms, the results of which will lead to a better understanding of the glaciological information. In order to estimate the 3-D flow rates of the glaciers, the two measurements along radar LOS from descending and ascending passes are not enough. An additional assumption regarding the glacier flow direction is needed. In this study, the following three different assumptions have been used, and corresponding 3-D flow fields have been generated: • the glacier flow direction is aligned with the medial-moraine line; • the glacier flow direction is parallel to its surface; • the glacier flow direction is along the greatest slope. Finally, a discussion is made to show which assumption is more favourable under different circumstances. 5 1.3 O u t l i n e o f t h e T h e s i s The remainder of this thesis describes the theoretical background, approach, and results of the research work to realize the above goals. Chapter 2 reviews the basic principles and theory related to SAR and InSAR. Most of the notations in the rest of the dissertation are established in this chapter. The geometries of SAR and InSAR are used to derive the height and displacement equations. This is followed by InSAR correlation which determines the amount of phase error. A special treatment is given to InSAR phase statistics and its applica-tions to predict the interferometric-phase variation. Both single-look and multi-look interferometric-phase probability density functions (PDFs) are derived based on the echo variables distribution. Some results of simulations show the relationship of the standard deviation of phase with the number of looks and correlation coefficient. These are very explicit indicators illustrating the tradeoff between InSAR processing efficiency and accuracy. Chapter 3 describes the overall InSAR processing steps to generate interfero-grams due to surface displacement. Based on the rationale of InSAR systems and data formation, several popular algorithms and implementations for each of the processing steps are presented. A brief discussion of the accuracy and efficiency of each algo-rithm is also included. Especially, all the algorithms used in this study are examined in detail. Chapter 4 is the core of this dissertation. It focuses on the issues of estimating the 3-D surface motion of the alpine glacier. Several descending and ascending inter-6 ferograms are presented t o show the m o t i o n - i n d u c e d and t o p o g r a p h i c phase pa t te rns . A sma l l p a r t o f th is chapter is devoted t o e x p l a i n i n g the t h e sources o f these phase pa t te rns . A measurement a long the cent re l ine of the glacier is m a d e t o ref lect t o i ts surface d isp lacement . T h r e e reasonable f low assumpt ions are used t o p ro jec t the LOS disp lacements on to surface d isp lacement . T h i s is fo l lowed by an i n - d e p t h analysis o f the resul ts der ived f r o m di f ferent p r o j e c t i o n approach. F i n a l l y , chapter 5 presents a s u m m a r y o f the research f ind ings a n d the signif-icance o f the i nves t i ga t i on , and some issues re la ted t o i n t e r f e r o m e t r i c e s t i m a t i o n o f glacier m o t i o n t h a t requ i re f u r t h e r s tudy. 7 C h a p t e r 2 R e v i e w o f I n S A R T h e o r y 2.1 I n t r o d u c t i o n Interferometric Synthetic Aperture Radar (InSAR) is the process of extracting geo-physical information from the phase differences between two or more SAR images. Unlike conventional SAR systems which totally depend on the magnitude part of the complex images to present the illuminated scene, an InSAR system mainly depends on the phase information. It is now gaining increasing credibility as a technique for rapid, accurate topographic data collection. According to the relative orientation of the two antennas, InSAR can be classi-fied as along-track InSAR (ATI) in which the antenna separation direction is parallel to the flight track and cross-track InSAR (XTI) in which the antenna separation direc-tion is perpendicular to the flight track. The along-track antennas enable Displaced 8 Phase Center Antenna (DPCA) techniques [19], and ATI can be used to measure the velocity of the moving targets on the ground or the ocean currents with radial velocity errors ranging from 1 cm/s to 3 cm/s [20]; across-track interferometry can provide pixel height estimates accurate to between 3 m rms and 10 m rms for airborne SAR and spaceborne SAR, respectively. Another promising application of X T I is to measure the surface displacement at a large scale. In this chapter, emphasis is only placed on revealing the important aspects of X T I system principle. 2.2 F u n d a m e n t a l s o f I n S A R T h e o r y 2.2.1 Basic principles of Synthetic Aperture Radar Synthetic Aperture Radar (SAR) is a coherent microwave imaging system, which was first introduced by Carl Wiley of Goodyear Aerospace Co. in 1951 [21]. It can produce high resolution images of terrain and targets. The main mechanism of SAR is to improve the azimuth resolution through analyzing and processing the Doppler properties of the received signal. Due to the capability of all time, all weather sensing, SAR systems are specially useful for monitoring areas of persistent cloud cover, determining wave conditions under storm systems, disaster monitoring, and monitoring temporal change. In addition, SAR data contains phase information of the observed area, which is essential for SAR interferometry to derive DEMs. This cannot be achieved by other imaging systems such as optical systems, and infrared 9 systems. The imaging geometry for a SAR system is shown in Figure 2.1. From this figure, we can see that SAR is a side-looking sensor with an imaged swath width from 10-70 kilometers for airborne systems to 50-500 kilometers for spaceborne systems. The swath is the area scanned by radar beam along the moving direction called azimuth. The parameters of a typical satellite-borne SAR can be represented by those of the current ERS-1 satellite, which are approximately given by: Parameters Value Unit Altitude H 780 Km Orbit period 98 min Effective velocity Vr 7050 m/s Nominal slant range R 850 Km Slant range swath 40 Km Ground range swath 100 Km Range bandwidth 17 MHz Chirp Length 34 LIS Range F M rate Kr 0.5 MHz/LIS Radar wavelength A 0.0566 m P R F 1700 Hz Antenna length 10 m Doppler bandwidth 1410 Hz Azimuth F M rate 2052 Hz/s Table 2.1: Radar Parameters of ERS-1 SAR The antenna transmits chirp pulses in the slant range direction to obtain higher range resolution. The Chirp Signal, modulated by a high carrier frequency fo, has a large bandwidth which is inversely proportional to the resolution. 10 Sensor Trajectory Figure 2 . 1 : The geometry model of SAR The equation of the transmitted signal is: st(ri,T) = P(T)cos(27r/0T + 7 r / ^ ( T - T / / 2 ) 2 ) , r = [0,n] where • P(T) is the envelope of the range pulse (usually a rectangular function) • r is range time (s) 11 • KT is the range F M ra te of the ch i rp (HzIs) • TI is the t i m e d u r a t i o n of the ch i rp (s) • fo is the radar carr ier f requency (Hz) • 77 is a z i m u t h t i m e (s) T h e equa t ion of the ideal received signal af ter coherent d e m o d u l a t i o n t o base-b a n d can be expressed as a comp lex signal w h i c h conta ins b o t h m a g n i t u d e and phase i n f o r m a t i o n o f the ref lector : sd(ri,r) = A(ri - n c ) P ( r - Td)exp(-j2iv f0rd + jirKr(r - rt/2 - rd)2), (2.2) where r is i n the in te rva l [rd — T ; / 2 , rd + 77/2]; 77 is [r/c — r/i/2, rjc + r/i/2]; j2 = —1, and rd = 2R(r))/c (2.3) R(V)2 = Rl + V 2 ^ - ^ ) 2 (2.4) w. here rd is the t i m e delay of the echo f r o m ref lector (s) R(rj) is the range f r o m an tenna t o ref lector at t i m e 77(771) c is the speed of l i gh t (m/s) 12 • A is the radar wavelength (rn) • Ro is the range to the reflector at time of closest approach (m) • Vr is the effective platform or radar velocity (m/s) • rjc is the azimuth time that the beam center crosses the target(s) • r/o is the azimuth time that the antenna is closest to the target (s) • rji is the synthetic aperture duration (s) • A(rj) is the two-way azimuth antenna gain pattern • cj) is the azimuth angle measured from the beam center {radians) • k is the beam width factor From equation (2 .2) , the phase of baseband received signal contains two parts, jirKR(T — T / / 2 — TA2 and — J ^ T T / O ^ • Because 77 changes much slower than r , we can think is approximately a constant in JTxKr{r — TI/2 — TJ) 2 , which represents a linear FM signal in range direction. Considering the phase part of —j2nfoTli and expressing the range equation approximately as a second order polynomial in 77, we can see the azimuth signal is also a linear FM signal with a negative FM rate. Thus, the SAR signal data can be viewed as a two-dimensional array of linear FM signals, indexed by slant range and azimuth time. Usually the major processing occurs on the ground after the received signal has been down-linked to Earth. So far, several different techniques have been developed to process the SAR raw data into an image, such as Range/Doppler, Chirp Scaling, 13 and S P E C A N [22, 23, 24]. In all techniques, matched filtering (or pulse compression) operations are applied in both azimuth and slant range directions. Pulse compression in the range direction is the concept used in traditional radar systems. Pulse compres-sion in the azimuth direction is achieved through matched filtering the azimuth phase information, using the phase history which a pixel on the ground has experienced as it is swept by the radar beam. The output of SAR processor is a two-dimensional complex image corresponding to the radar reflectivity of the observed area on the ground. 2.2.2 InSAR geometry and height equations The theory of InSAR has already been presented in detail by Rodriguez et al. [11]. Here we will only summarize the main results and establish notations. For repeat-pass InSAR systems, on each pass the radar acts as both a transmitter and receiver, therefore the total path difference for each radar observation to a given pixel on the ground is twice the range from the antenna to the ground. So, the relationship of the phase value with the slant range will differ by a factor of two, <& = 2R/X. While only a single set of raw data is used in conventional SAR, InSAR systems make use of two sets of raw data. We can consider the InSAR system geometry in Figure 2.2. Arbitrarily, we call one of the two images the master image (M) and the other the slave image (S). We can see that two pixels located at the same range in the master image do not have the same slant range as in the slave image if the elevation of the corresponding ground points are different. Because the phase value </5 is strictly 14 related to the range which the microwave travels, we can estimate the elevation of the ground pixel by using the phase and the geometry information. Figure 2.2: The geometry model of InSAR To analyze the system geometry, we assume the ground surface has topography described by Z(y), the baseline length is B, the altitude of the master track is H, the incidence angle of master track is 9, the distance between the master track to the observed pixel is R, the elevation angle of the baseline is a, the difference of the distance from master track and slave track to the observed pixel is 5. Ignoring phase 15 noise, the phase value of the master image is: 4>M = 2R 2TT T (2.6) where A is the wavelength of the radar system. Similarly, the phase value of slave image can be expressed as: 4 > s = 2(R + 6) 2TT T ' (2.7) where S is the difference in slant range between the master and slave images. The phase difference of the two S L C images can be obtained from equation (2.6) and equation (2.7): In interferometry, the phase difference is the phase of the product between one image and the complex conjugate of the other image. From the above equation, we can see that the phase difference for the observed point P is proportional to the range difference. Apply ing the law of the cosines to the triangle formed by the master track M , slave track S, and pixel P results in: In equation (2.9), all the parameters except 9 can be obtained from I n S A R system. Hence we can express 9 as a function of the other parameters. After 9 is determined from the above equation, we can easily get the elevation of the observed pixel P: (2.8) a) (2.9) (2.10) 16 2.2.3 The relationship of interferometric phase with baseline and topography To i l l u s t r a t e the effect o f baseline and topog raphy on i n t e r f e r o m e t r i c phase, several a p p r o x i m a t i o n s have to be used. I n th is sect ion, t h e basel ine is separated i n t o t w o t e r m s , n o r m a l basel ine, Bn, and para l le l basel ine, Bp, by us ing the d i r e c t i o n o f the radar cen t ra l look angle as a reference. T h e t o p o g r a p h y i n f o r m a t i o n is s t i l l def ined by t h e h o r i z o n t a l d is tance between the sate l l i te a n d t h e ta rge t , y , a n d the target he igh t , z. Here the hor i zon cou ld be a curve para l le l t o the surface o f the e a r t h i f a more accurate m o d u l e is used t o describe t o geomet ry . F r o m Figure 2.2, the baseline components are re la ted to the t i l t angle and t h e t o t a l basel ine l e n g t h b y Bn = Bcos(6c-a) (2.11) and Bp = Bsin(9c-a) (2.12) Equation 2.9 can be expressed as Bsm(0-a) = ^ + ^ ~ ^ ( 2 -13 ) Because S is always far smal ler t h a n t h e d is tance f r o m t h e sate l l i te t o the ta rge t , t he above equa t ion can be a p p r o x i m a t e d as Bsm(O-a) « 5 - — (2.14) 2 r t T h e b e a m e levat ion angle can be expressed as the s u m m a t i o n : 6 = 0C + 9d, where 9d is t h e e levat ion angle d e v i a t i o n about a n o m i n a l e levat ion angle 6C, w h i c h 17 yields Bsm((9c-a) + 9d) = Bv cos 9d + Bn sin 9d (2.15) Combining Equation 2.14 and Equation 2.15 and solving for the approximation range difference, we have B2 S « — + Bnsm6d + Bpcos6d (2.16) ZK For most spaceborne SARs as E R S l / 2 , the deviation angle, 9d, only varies around several degrees across a swath. Thus the range difference is almost a linear function of 9d, with its slope determined by Bn. Referring to Figure 2.2, nc is defined as the distance normal to the radar central slant-range direction from the target point. This distance is related to the deviation angle and group range by the following two equations: sin8d = ^ (2.17) anc nc = (y - yc)cos9c + zs'm9c (2.18) Inserting these two equation into Equation 2.16 and combining Equation 2.8 yields 2TTB2 ATTB $ « -j^- + ^^((y -yc)cos9c + zsm9c) + Bpcos9d (2.19) This equation shows that the effect of topography and ground range on the interferometric phase is dominated by the normal baseline length, Bn, with only small, nonlinear terms affected by the parallel baseline length, Bv. From this equation, it 18 is also obv ious t h a t the i n te r fe romet r i c phase is g iven by t h e f o l l o w i n g f o r m u l a w h e n the ta rget is a f la t surface (i.e. z = 0) . 47T Bn $|z=0 « " y cos 0C + Bv cos Od + const (2.20) I t is apparen t t h a t a a p p r o x i m a t e l y l inear phase r a m p exists i n a i n t e r f e r o m e t r i c phase image even on a flat surface. I n order to o b t a i n t h e phase i n f o r m a t i o n caused by the topography , th is phase r a m p m u s t be sub t rac ted t o e x t r a c t the phase due to t o p o g r a p h y : 47TB n s in 9C , n . $2 « ^- -z + const (2.21) AR F r o m the above a p p r o x i m a t i o n , i t is clear t h a t t h e sens i t i v i t y o f t h e in te r fe ro-m e t r i c phase t o topog raphy is d i r e c t l y p r o p o r t i o n a l t o Bn, and inversely p r o p o r t i o n a l to R. I n order to achieve more accurate e s t i m a t i o n of t h e topography , longer n o r m a l basel ine are needed. I f we wan t to m a i n t a i n the same sens i t i v i t y , an in te r fe romete r o p e r a t i n g at larger range requires a longer basel ine. I t is w o r t h t o keep i n m i n d t h a t long baselines or sma l l R can m a k e phase u n w r a p p i n g d i f f i cu l t , a n d exaggerate the effect o f topography . 2.2.4 Differential interferometry Di f fe ren t ia l i n t e r f e r o m e t r y is a very in te res t ing top ic due t o i ts a b i l i t y t o measure t h e surface d isp lacement w i t h qu i t e a h i g h accuracy. U n d e r t h e assumpt ion t h a t t h e d ie lec t r ic character is t ics o f the g r o u n d r e m a i n constant a n d the o rb i t s sat isfy t h e cond i t ions necessary for coherence, the phase i n f o r m a t i o n i n an i n t e r f e r o g r a m 19 contains two parts: (1) topographic information, and (2) any changes in position of the ground reflector between the acquisition times of the two SAR images. The task of differential interferometry is to separate the phase parts due to topography from those due to ground displacement, then to get the information about surface movement of the observed area. An accuracy of 2 — 3 cm has been achieved by using SEASAT L band SAR [25] or ERS-1/2 [26] in measuring surface movement. There are two methods to realize differential interferometry: (1) Removing the phase part due to topography by using elevation information from a D E M , where the D E M is transferred into an purely topographic interferogram and differentiated with the interferogram obtained from SAR images; (2) Subtracting topographic phase information through double-differencing, where two interferograms made from three (or more) SLC SAR images collected at different times are differentiated once again to produce a third interferogram. In the first case, the D E M of the observed area is first transferred into a slant range interferogram using the same geometric parameters as the InSAR system where the other interferogram was obtained. Then these two interferograms can be differentiated to get the phase information due to the pure displacement on the ground. The key points are have an accurate D E M and to obtain fine co-registration between the two data sets. But the DEMs of the examined areas are not always available or with required accuracy. Then the second method can be used to derive the topographic information if the surface velocity field is constant over the span of the observation sequence. The problem for this method is that the current spaceborne SAR orbit repeating time is 20 at least one day. So this method is only suitable for the use in areas with a relatively high surface stability. For the case of mountainous glaciers as the Lowell, which is moving quite fast, we have to rely on the available DEMs to extract the pure surface motion phase. The general theory of detecting surface change from multiple orbits was described by Gabriel, Goldstein and Zebker [14]. A simple geometry associated with three orbits observation is shown as Figure 2.3. to target \ to target Figure 2.3: Differential interferometry InSAR geometry 21 Assuming that the radar rays of the three antennas are parallel, and the surface of observed area moves with a constant velocity during the multiple passes over the flow measurement, then the phase difference between first and second observation can be approximately presented as [13]: $ 1 2 = ^ p s i n ( 0 - a ) + ^ t ; i i 2 (2.22) A A Similarly, the phase difference between the second and the third is: $ 2 3 = y g s i n ( 0 - / ? ) + y i ; < 2 3 ( 2.23) where 9 is the incidence angle of radar LOS; v is the surface velocity along LOS; ti2 is the time difference between the first and the second observations; £23 is the time difference between the second and the third observations. Now consider the situation of two interferograms acquired over the same re-gion with the same time difference (i.e. ti2 = £23)- The difference interferogram gives a new interferogram, $ 1 3 , where the contribution of the displacement has been removed. Here we use a prime to indicate that there is no phase component due to displacement in this resultant interferogram. $ ' 1 3 = $ 1 2 - $ 2 3 ( 2-24) 47T = —(psm(9-a)-qsm(9-f3)) (2.25) A 47T = — r s i n ( 0 - 7 ) (2.26) A From the above equation, it is clear that a topography-only interferogram can be derived from two mixed topography/motion interferograms. -The sensitivity 22 of the interferometer to topographic variation is now dependent on r. Of course, differencing the two interferograms requires the phase be unwrapped to resolve the 2ir phase ambiguity. Then $ 1 3 is to be subtracted from one of the mixed interferogram to obtain a pure motion interferogram. 2.3 Correlation in Interferometric Images In InSAR technology, one of the basic processing steps is to create a phase difference image called an interferogram by multiplying one single look complex (SLC) SAR image with the conjugate of another well-registered SLC SAR image. Research work done by Rodriguez and Martin [11], Bamler and Just [27], Zebker et al. [28] show that the accuracy of elevation and differential-motion estimates made from this inter-ferogram depends on the amount of correlation between the two SLC SAR images. It is very important to understand the relationship between system parameters and the overall InSAR image correlation or coherence magnitude. In an InSAR system, there are many factors causing decorrelation in InSAR images. Zebker and Villasenor [29] list five major multiplicative sources, each char-acterized by a corresponding correlation coefficient: 1. Temporal decorrelation, ptemporai, is caused by the changes in the relative posi-tions of the scatterers between passes for repeat-pass InSAR. 2. Baseline decorrelation, pbaseiine-, results when the target is viewed from slightly 23 different positions. 3. Volume scattering decorrelation, pvoiume, occurs when the penetration depth of radar beam through the target medium varies from different incidence angles. 4. Thermal decorrelation, pthermal, is due to the presence of uncorrelated-noise sources. 5. Processing decorrelation, ppr0cessing, is introduced by the processing errors such as misregistration and interpolation. The overall correlation coefficient can be presented as the product of these individual correlations [4]: P = Ptemporal X Pbaseline X Pvolume X Pthermal X Pprocessing (2.27) Rodriguez et al. [11] and Joughin [30] gave the derivation of the total correlation coefficient in detail based on the following reasonable assumptions: • Each resolution cell consists of large number of discrete scatterers. • Resolution cell size is much larger than the radar wavelength. • The reflectivity function has a flat spectrum (spatially white). • The SAR imaging function is separable in the range and azimuth directions. • The average backscattering cross section varies only as a function of depth. 24 Based on these assumptions, the received signal for each pixel is given by V = \V\e~^ = J2^e-j4'' (2.28) i where V{ and fa are the amplitude and the phase of individual scatterer respectively. Under the Central Limit Theorem, V is a complex Gaussian random variable with Rayleigh-distributed amplitude, | V | , and uniformly distributed phase This is the well-known speckle phenomenon of SAR images. For a distributed target, a pixel in a complex image can be represented as: Vx = expi-jAnro/X^V^-^, (2.29) where r0 is the range for the radar sensor to the ground pixel; $1 is the phase noise; exp( — jAnro/X) is due to the additive phase caused by the round trip of the electro-magnetic wave. A complex interferogram is formed as: ViV* = exp(-jATrr0/X+jAw(r0 + 8)/X)\V1\\V2\e^l+^ (2.30) where 5 is the distance difference between the two passes along radar LOS. The phase of the above interferogram is: $ = (47rA/A + ( $ 2 - $ 1 ) ) r n o d ( 2 T ) . (2.31) Normally phase noise $ i and $ 2 are uniformly distributed. But if the two SAR images are correlated, the phase noise difference, ($ 2 — $ i ) , will not be uniformly distributed any more. In fact the distribution can be quite sharply peaked at the 25 point of ($2 — $1) = 0 if the two images are highly correlated. This means that the interferogram phase mainly reflects the range difference of the same target between the two views. The higher the correlation coefficient of the two images is, the sharper the peak of the distribution of the variable ($ 2 — $1) w m be. The correlation coefficient of the two SAR images will definitely affect the accuracy of interferogram phase, thus the measurement accuracy of topography and motion. Normally, in order to achieve a subcentimeter accuracy in measuring the displacement along the radar LOS, the overall correlation coefficient should be greater than 0.4 [31]. The factors affecting the correlation of the two SAR images are discussed in the following subsections. 2.3.1 Temporal correlation coefficient Temporal correlation coefficient shows the time-variant volume-backscattering prop-erty of the ground targets [29]. qe(0) / ae(z)dz Ptemporal - ^ J ^ ) where a is the volume-backscattering coefficient, ae is the effective volume-backscattering coefficient, which is also a function of time implicitly, z is the penetration depth of the radar microwave. The volume-backscattering coefficient is affected by the positions and backscat-tering properties of the scatterers within the imaged medium between passes. En-vironment and characteristics of the imaged area dominate temporal decorrelation. Zebker and Villasenor [29] have analyzed the InSAR images generated from SEASAT 26 data and found that decorrelation coefficient of vegetated area is far more sensitive to the lapse between passes than desert area. They also found that 2 to 3 cm ran-dom motion is enough to cause complete decorrelation in C-band images. Recent researches by Valero [32] and Joughin [33] show that high correlation still can be achieved with several tens of centimeters motion between two passes. The main rea-son is that glaciers and ice sheets scatterers roughly remain in the relative positions during the movement. Due to all of these, the temporal decorrelation coefficient is not easy to be quantified accurately. But an approximate conclusion can be drawn here that corre-lation coefficient will decrease if the time interval between two observations increases to a certain extent. 2.3.2 Baseline correlation coefficient It is not hard to see that the existence of a non-zero baseline will add an additional phase to the interferogram. Baseline decorrelation occurs because the different view-ing angle of the radar sensor causes a change of the effective scattering center of each pixel. Further study conducted by Joughin [33] shows that this additional phase part is proportional to the normal baseline, Bn. A general presentation of pbaseline can be: P b a s e H n e ~ fWrA(8)W*2(8)dB [ Z - 6 6 ) where Wr is the range part of the SAR imaging function; exp(—j^j~^jf^) 1S the additive phase part due to normal baseline; r = r i ^ r 2 is the average of the slant ranges from the two sensor positions to the ground pixel, 3 = rj — ro; ro is the 27 nominal slant range in the resulting complex image. From the above equation, it is clear that the baseline correlation strongly de-pends on the length of normal baseline, Bn, and range resolution (width of Wr). Figure 2.4 shows the relationship of correlation coefficient with the length of nor-mal baseline with typical ERS parameters [34]. The correlation coefficient almost decreases linearly when the normal baseline length increases. Thus in order to obtain high enough overall correlation coefficient, ERS-1/2 normal baseline should be less than 300 m. 28 2.3.3 Volume correlation coefficient Volume scattering decorrelation occurs because the penetration depth of radar rays into a 3D medium varies with the incidence angle of the radar beam. The volume correlation coefficient, pvoiUmei is mainly dominated by the effective volume-scattering coefficient, ae, and the normal baseline, Bn, of the InSAR geometry. The general expression reflecting the relationship is as follows [34]: _ fqe'z)exp(-jg£sz)dz fae(z)dz The volume-scattering coefficient usually depends on the properties of the scat-tering medium in a very complicated way. But a reasonably simple model based on a statistically-homogeneous random medium can be [35]: O-AZ — < 0 z > 0 where L is the maximum penetration depth of the radar rays. Using the above model of effective volume-scattering coefficient, we can derive the simplified expression of volume correlation coefficient as: Pvolume — . _ • 2nLBn (2.35) 1 3 A r t a n f l Figure 2.5 represents the volume correlation coefficient as a function of penetration depth and normal baseline. Typical ERS-1 parameters (r = 860 km, 6 = 23°, and A = 5.6 cm) are used to generate these plots. Usually the penetration depth varies from several tens of centimeters for wet surfaces to several meters for dry surfaces. 29 It is easy to observe that volume scattering decorrelation for normal baseline longer than 300 meters with the penetration length of 5 m. This also serves as a guideline to what the suitable normal baseline we can use. O 5 10 15 20 25 30 35 40 45 50 L(m) Figure 2.5: Correlation coefficient as a function of penetration depth L and normal baseline Bn. 2.3.4 Thermal correlation Thermal correlation, by the name, is mainly due to the thermal noise of the radar receiver. The thermal-noise correlation coefficient can be expressed in the form [29]: Pthermal Z ' J (2.36) i SNR where signal strength is determined by the integrated scattering coefficient and the transmitter power. The relationship of pthermal and SNR is shown as Figure 2.6. 30 0 5 10 15 20 25 30 35 40 SNR Figure 2.6: The thermal correlation coefficient, pthermal ,as a function of SNR 2.3.5 Processing correlation There are many factors during the processing which will affect the correlation of the two SAR images. We can generally classify them as registration decorrelation, pregistrati, and Doppler centroid mismatch decorrelation, pmiSmatch-Regis t ra t ion correlation When the two SAR images are not perfectly registered, decorrelation results. The cor-responding correlation coefficient can be separated into coefficients for misregistration 31 in range and in azimuth. The coefficient is given by [35]: iWrMW^B - 8r)exp{-j •_47LBn 3) dp J ^ , 1 ( 7 ) ^ ( 7 -Sx) d1 Ar tan 6 fWrtl(B)W*2(8)exp(-j B)d3 / W , l l ( 7 ) W * 2 ( 7 ) d 7 (2.37) Ar tan 8 The curve reflecting the relationship of the correlation coefficient and the reg-istration error is a sine function, which reflects the shape of the SAR image func-tion [21] . Using the typical parameters of ERS, results are plotted as Figure 2.7 and Figure 2.8 [35]. Hagberg et al. [34] gave the similar results for the processor used by United Kingdom Processing and Archiving Facility (UK-PAF) . In U K - P A F , "raised-cosine" spectral weighting function is used to depress the relatively large sidelobes of sine function. In the meantime the main lobe of the sine is widened. As result, the correlation coefficient will be higher for the same small registration errors (i.e. 8X and 5r are less than 0.5 pixels). Of course, this correlation gain is at the expense of the resolution. From these figures, we can see that in order to achieve a reasonable level of correlation, registration errors should be kept less than a few tenths of a pixel. Doppler centroid mismatch It is a common practice to estimate the Doppler centroid during SAR processing for proper focusing the SAR data. The Doppler centroid depends on the squint angle of the SAR, which is the angle by which the direction of the beam center differs from the direction normal to the flight path. Different squint angles for each orbit result in different Doppler centroids, thus different azimuth center frequencies for the received 32 8 r (pixels) Figure 2.7: Correlation coefficient, Pr-registration, as a function of range direction registration error data. During the processing of InSAR, this spectral misalignment is related to the imaging functions by [35]: WXi2(x) = WxA(x)exp(-j(QDCil - flDC}2)x) (2.38) where ODC,I and 0/jc,2 are the Doppler centroids for the first and the second images. The correlation coefficient is given by [35]: /\WXil(j)\2exp(-j(nDCtl - 0 D c , 2 )7) d l Pmismatch (2.39) Figure 2.9 shows this function for the square weighting function. It is not hard to see 33 1 0.9 0.8 c 0.7 g trat 0.6 w 'cn CD 0.5 a 0.4 0.3 0.2 0.1 0 ~ ~ ~ ™ — ' ' 1 1 1 " 1 " i 1 1 " T 1 i 1 1 T 1 i i r " I i " r i I I n I I - I -\ 1 - \ 4 -1 - -L 1 1 J-1 I i- -i - u I - i i I 1 " " "1 1 \ T \ 1 1 T 1 i r _ I - r i ~i i 1 J - 1 - J -1 \ 1 l \ 1 4-1 I L -I - i _ I j i 1 " I " 1 " 1 ~ 1 1 \ - T 1 1 V r I r • I - r i 1 J . 1 . _1 . 1 . 1 1 i \ I L _ . I . u j i 1 ~1 " 1 " 1 -1 1 - T 1 i r \ I - r i - i I 1 J . 1 . _1 . 1 1 . 1 1 i i I \ _ L I j i — " j — " j -1 f - -1 1 1 1 1 i - t -i i I _ _ ! I I i 3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 c 8 x (pixels) Figure 2.8: Correlation coefficient, px-registration, as a function of azimuth direction registration error that the shape of this curve is similar to that of baseline correlation in Figure 2.4. The reason is that baseline decorrelation is also a form of spectral misalignment. For baseline decorrelation the misalignment is caused the by the exponential factor, which causes a relative shift of the spectra in the range direction. For either baseline decor-relation or spectral misalignment, the expression for the correlation coefficient can be expressed as the area under the product of the misaligned spectra normalized by the area under product of the un-shifted spectra. As the misalignment increases, the area under the product decreases, causing a decrease in correlation. For a square-weighting function, this decrease is linear. If we the raised cosine weighting function, the de-34 crease will be slow at first. Then we will suffer less decorrelation when the spectral mismatch is not very big. Of course, the negative effect is the loss in resolution. 1 0.9 0.8 0.7 0 0.6 (0 E 0.5 CO E 0.4 a 0.3 0.2 0.1 0 \ l I N . I \ I I i i 4 1 1 1 1 1 i i i 4-1 1 I I T I I 1 1 1 X ! 1 1 1 1 1 1 -1 l \ 1 1 \ 1 1 r L 1 1 T 1 1 -L 1 1 1 r i i i 1 1 1 1 " 1 1 1 1 1 1 1 1 r ^ 1 1 v t \ ^ 1 1 \ 200 400 600 800 1000 1200 ^ D c , 1 " ^Dc ,2 c-mismatch, a s a function of Doppler centroid Figure 2.9: Correlation coefficient, prj mismatch Bamler and Just [27] showed that decorrelation caused by spectral misalign-ment or baseline decorrelation can be reduced or eliminated at the expense of resolu-tion. This is accomplished with additional filtering to remove the spectral components that do not overlap. This means that only the coherent parts of the spectra are re-tained, increasing correlation. The filter on the other hand, reduces the bandwidth of the signal, and hence the resolution also. Usually we can always choose a SAR image pair with normal baseline less than 200 meters so that the decorrelation due 35 to the normal baseline can be ignored. But the Doppler centroid mismatch can be as serious as 500 Hz [27], thus this type of decorrelation must be considered. The current matlnSAR processor in the RRSG of U B C has used the azimuth filter to reduce the spectra decorrelation. As mentioned in the previous section, multilook processing technique is usually ap-plied to reduce the statistical noise speckle at the expense of losing resolution. In InSAR processing, it is very important to smooth the interferogram before phase un-wrapping. Given a level of coherence magnitude, multilook processing will reduce the phase noise significantly. The main idea of multilook processing is to average several complex pixels in the interferogram to create a new complex pixel. Lee et al. [36] gave the expression for normalized estimation of corelation as: where the number of looks, n, is the number of complex interferogram pixels averaged. The general expression of the probability density function (PDF) for a n look interferogram as a function of phase error was also given by Lee et al. [36]: 2.4 Multilook Technique p = z\n=xSi(k)s;(k)\ (2.40) p(V0 r(n + i/2)(i-H2r/? ( i - H 2 ) 2 v ^ F r ( n ) ( l - / ? 2 ) n + l / 2 + 27T F(n,l;l/2;82) 7T < < 7r; (2.41) with 8 = \p\ cos(tp — 9) (2.42) 36 where F(n, 1; 1/2; (32) is a Gaussian hypergeometric function, 8 is the mean phase for the averaging region. An alternative expression for the phase-error distribution in higher order hypergeometric functions was independently derived by Lopes et al. [37]. The P D F of equation (2.41) depends only on the number of looks and the correlation coefficient. The peak value of the distribution is located at ip = 8. As 8 varies, the center pointer (i.e. the mean value ) of the distribution will shift with the shape of P D F remaining the same. Actually, a fringe on the interferogram is defined by 8 varying from — ix to + 7 r . For small n, the Gaussian hypergeometric function can be expressed by trigonometric and algebraic functions. Lee et al. [36] gave the close-form phase-difference PDFs for n = 1,2,3,4. Figure 2.10 shows distributions for n = 1,2,3, and 4 with \p\ = 0.7 and 8 = 0. We can see that multilook processing improves the phase accuracy. It can be shown that when \p\ = 0, the P D F is uniformly distributed, which implies the interferogram phase does not contain any information of the observed area (no phase fringes appear on the interferogram); also when \p\ = 1, the P D F becomes a Dirac delta function, which means the phase error is zero. To see the relationship between the phase error and the magnitude of correla-tion coefficient, a plot of the phase-error standard deviation versus \p\ with different number of looks is shown as Figure 2.11 [36]. It is evident that when the number of looks increases, the phase-error standard deviation will decrease monotonically. 37 38 39 Chapter 3 I n S A R Processing Algor i thms for Surface Displacement Measurement 3.1 Introduction The complexity and number of processing steps of an InSAR system depend on the applications and required accuracies. Because the purpose of this study is to inves-tigate the feasibility of measuring the 3-D motion vectors of the Lowell glacier using InSAR technology, we will focus on the related processing procedures and algorithms in this chapter. Figure 3.1 shows the processing scheme to realize InSAR surface displacement estimation and multilook averaging is an optional processing step. 40 Cadidate Image Reference Image Registration Offset Estimating Resampling The Candidate Image Forming Interferogram Removing Phase Ramps Geocoding/ Removing Topo Phase Multilook Averaging Calibrat ion ( op t iona l ) Phase Unwraping 3-D Motion Vectors Figure 3.1: Block diagram of the InSAR processing steps to realize surface displa ment measurement. 41 The algorithms to generate interferograms to be examined in the next chapter are discussed here following the same order of Figure 3.1. As shown in Figure 3.1, the very first step of InSAR processing is to extract two SLC SAR images. This normally involves an extensive amount of computation and data manipulation for the simple reason that the formation of each pixel involves the coherent combination of data gathered along the flight path. SAR processing algorithms such as Range/Doppler, Chirp-Scaling [21] are commonly used to gener-ate SLC SAR images. An ideal SAR processor should have both high computation efficiency and processing accuracy. But for InSAR application, phase preservation in the processing is of utmost importance. In this research, Range/Doppler is used to obtain the SLC SAR images. Because the emphasis of this study is on the interfero-metric processing of SAR information, the detail of Range/Doppler algorithm is not presented here. After the two correlated SLC SAR images are available, the InSAR processing starts with the registration of these two images. In this step, one image is considered as the reference while the other is considered as the candidate. Then a certain al-gorithm is needed to estimate the registration offsets, which are the shifts needed to register the images of an interferometric pair. From Figure 2.7 and 2.8, in order to achieve a acceptable level of correlation, the registration errors should be kept under a few tenths of a pixel. These offsets are used to interpolate the candidate image to align with the reference. This is realized efficiently by a finite impulse response (FIR) interpolator. Section 3.2.1 is devoted to discussion of the common algorithms dealing with the registration offset estimation. Section 3.2.2 shows the implementation of 42 the FIR interpolator. It is always necessary to remove the azimuth and range phase ramps (i.e. flattening) from the interferograms before further analysis. Azimuth phase ramp is caused by the gradual baseline length change along the satellite track, while the range phase ramp is due to the flat earth surface. This is investigated in Section 3.3. Because the phase variance in a single-look interferogram is usually too large, multilook average is applied at the expense of resolution. In this study, D E M information is used to remove the phase due to topography. First the D E M information in a certain map format is converted back to a simulated slant-range interferogram, then it is registered to the candidate interferogram. A geocoding algorithm is used to deal with this conversion. Special attention should be paid to the accuracy of baseline length and orientation estimation. Not only does it affect the removal of the phase ramps, but also the interpretation of the interfero-grams. By selecting a certain geographic area with accurate D E M information as a reference, InSAR baseline can be estimated. Section 3.6 describes the phase unwrapping algorithms used to remove the modulo-27r ambiguity. Several popular algorithms are discussed in this section. 3.2 Registration Usually the two SLC SAR images are not aligned perfectly. To satisfy the requirement of interferometry, two steps are needed to register the two SLC SAR images within 43 a fraction of a,pixel. Estimation of registration offset and interpolation are discussed in detail in this section. 3.2.1 Registration offset estimation Due to the random factors of the satellite mechanical system, the track and the look angle between ERS -1 and ERS-2 can experience minor changes. Offsets up to a few hundreds of pixels between a pair of images always exist. As E R S -1 /2 are side-looking SAR systems, each image is obtained from a slightly different look angle, thus the size of the ground-range cells differ by a small amount at each point. As a result, the offset in range direction varies with range. Azimuth offset can also vary with range (i.e. range-dependent azimuth offset) if the orbits of ERS -1 and ERS-2 are not parallel. A l l these complicate the offset estimation procedure. In order to achieve maximum computational efficiency, we use a multi-scale approach to realize the whole registration process. First of all, large offsets between an image pair are manually estimated by marking a few identifiable features found in each image. Normally after this procedure, the pair of images can align with each other to within a few pixels. Of course, when the investigated area is featureless terrain such as a flat ice-sheet, the method does not apply. After the pair of SLC SAR images are roughly aligned, the more accurate esti-mation of the offset is normally realized by locating the peak of the cross-correlation between the image pair. In order to determine the offset dependence on range and azimuth, this procedure should be repeated at a few points evenly distributed through-44 out the image. Then a linear or quadratic fit is used to find the offset function between this image pair. The effect of phase ramps on cross-correlation peak According to the discussion in Section 2.2.3, there is an almost linear phase ramp in the range direction contributed by the flat Earth surface. It is also common that the length of Bp will change by a few meters over the frame of 100 Km. This variation of Bp can also create a significant phase ramp in the azimuth direction because a 2.8 cm change in the radar line of sight yields a 27r phase difference. Experiments show that it is good enough to just approximately consider the azimuth phase ramp as a linear phase ramp [2]. In this section, we will show that the presence of these phase ramps make it difficult to detect the cross-correlation peak between an image pair. To illustrate this effect, we can analyze the cross-correlation between two one-dimensional complex signals: fi(t) and /2(f) = / i (£)e J U ' t . Here to denotes the slope of the phase ramp. Then the cross-correlation of these two signals can be expressed as where L is the computation length of the signals. Figure 3.2 shows how the cross-correlation peak is blurred with an increasing phase ramp slope. So the phase ramp in the InSAR image will definitely make it more difficult to accurately estimate the shift between the two SAR images, and result in a registration error. R(t) (3.1) L/2 45 J I I I I I L 0 1 2 3 4 5 6 7 8 9 10 Slope of phase ramp, oo (n) Figure 3.2: The relationship of the cross-correlation peak with the slope of the phase ramp. Offset est imation When the phase ramps are known and significant, these phase parts should be sub-tracted from the SAR images. If accurate baseline length and orientation are available, the phase ramps can be estimated using Equation 2.20 and Equation 3.8. Then one can use the cross-correlation function to accomplish the registration. Normally the baseline is either unknown or inaccurate. Without compensation for phase ramps, many ERS-1/2 images cannot be correctly registrated using the cross-correlation function [38]. ESA publishes a list of baselines for ERS-1/2 interfer-ometric pairs. This baseline can be used to estimated the range and azimuth phase 46 ramp. In this study, a two-stage offset estimation is used to register the images for the interferograms. At the first stage, both the reference and candidate data is subsampled by the factor of 2 in the range direction and 8 in the azimuth direction respectively. The reason for these unbalanced subsampling ratios is that the azimuth resolution is 4 times higher than the range resolution for ERS-1/2 SARs. A reference chip with the size of 128 x 128 and a candidate chip with the size of 65 x 65 are selected. The initial central points of these two chips are overlapped. Then the candidate chip is shifted inside the boundary of the reference chip. 64 x 64 cross-correlations are calculated to estimate the offset. The accuracy of this stage is ± 2 pixels. Once the rough offsets have been determined to the nearest pixel, the second stage is used to estimate the offset to the accuracy of a quarter pixel. First, the input reference and candidate data are oversampled by the factor of 2 in both azimuth and range direction. Then the candidate data is compensated with the estimated offset. After all these, the registration is performed by estimating the shift between smaller chips (20 x 80) evenly distributed through the images. Four different rules are used to decide whether the correct estimation has been made: 1. Maximizing correlation of image magnitude. First of all, interpolate the two images. Then several chips of images are extracted at regular intervals across the scene. The chips are shifted with respect to each other by fractions of a pixel and the correlation of the image magnitude is calculated every time. We take the estimated shift as the needed shift when the local correlation is 47 maximum. The more the extracted chips, the higher the accuracy, but the longer the processing time. 2. Maximizing coherence magnitude. This approach will consider both the mag-nitude and phase information of the images. Coherence magnitude can be measured using equation 2.40 with n = 1. The rest of the processing is the same as method one. 3. Minimizing residue counts. Residues refer to points in an interferogram where there are phase discontinuities. Under ideal conditions, no noise or layover, the number of residues should be zero. But it is not the case for real InSAR systems. Different kinds of sources for phase noise and rough terrain will always cause residues. It is not hard to see that the number of residues will be fewest when the two images are exactly registered with other conditions remaining constant. So, very similar to method one, we count the number of the residues for each small shift then can find the best shift. 4. Maximizing fringe visibility. As the phase values of interferogram are measured with a modulo of 2TT, SO the interferogram before phase unwrapping is composed by a series of fringes. These fringes should be most visible when the noise due to misregistration is canceled (two images are correctly registered). So following the same step as method one and using the peak value of the Fourier transform of the phase fringes as the quality measure, we can get the fine estimation of the shift. Although the criteria of finding offsets in each program are different, the final 48. candidate image shifts are all found via a grid search where the optimal conditions are met. Of course, denser searching grid will increase the registration accuracy. But Figure 2.7and Figure 2.8 show that no significant improvement is gained after the accuracy reaches 0.1 pixel. So in order to get best compromise between efficiency and accuracy, it is normally enough to register the two image to around 0.1 pixel. 3.2.2 Interpolation After the fine offset estimation and the determination of the range/azimuth offset functions, the candidate image is shifted to align with the reference image. This is realized by resampling the candidate image. The offsets can also be separated into integer part and fractional part. The integer shift can be done by direct shifting the rows and columns of the image without any additional decorrelation. The fractional part must be done by using interpolation. According to sampling theory, as long as a bandlimited analog signal is sampled above the Nyquist rate, the original signal can be reconstructed perfectly from its samples. Thus it is possible to interpolate the signal without error. Assuming x{nT) is the samples of the signal x{t) sampled with a ratio of fs = ^ which equals or exceeds the Nyquist rate. The reconstructed signal xr(t) relates to the samples through the following infinite summation. 00 t — nT xr(t) = x(nT) sine (—-—) (3.3) 49 A resampled version of x(nT) can be derived from xr(i) with shift of ST. xs(mT) xr((m — S)T) (3.4) x(nT) sine [( oo (m — 8 — n)T f (3.5) = x(m) (^) hs(m) (3.6) where hs(m) sine (m — 5) (3.7) From the above equation, it is clear that in order to obtain the resampled version of x(m), an infinite convolution is necessary. This is obviously infeasible in practice. So hs(m) is always truncated to a certain length. Here comes the typical engineering trade-off: A too long truncated length will increase the computation burden and the additional computation yields little gain; and a too short truncated length will cause significant decorrelation. Jougnhin et al. [35] shows that for a certain length of kernel length, the interpolation decorrelation increases with a larger kernel length. Normally when the interpolation kernel length is larger than 12, the decorrelation caused by the truncation can be ignored for fractional shift up to 0.5 pixel in both azimuth and range directions [35]. Because a shorter kernel length can be tolerated for smaller fractional shifts, it is important to note that a truncated sine function may not be the optimal interpolation kernel for a fixed length. As a result, the kernel length can be a function of the fractional shift value to maintain a constant level of correlation. In this study, the interpolation filter length for both azimuth and range di-rections is 15. Because of the investigated area is only a small part of a swath, 50 only one range interpolation filter is used to handle all the range lines. Different azimuth interpolation filters are generated for every azimuth line to best compensate the range-dependent azimuth shift. After the candidate image has been registered to the reference image, an inter-ferogram is formed by multiplying each pixel of the reference image with the complex conjugate of the candidate image. 3.3 Flattening As discussed in the last section, the phase ramp caused by the flat Earth has an adverse effect on the registration procedure. This part of phase has also to be re-moved before we can proceed to obtain the topography only or motion only phase image. Phase ramps are removed by multiplying the interferogram with a complex exponential that has the phase ramp as its exponent. The azimuth phase ramp slope can be estimated from the parallel component of the baseline at two different points along track. This slope is given by Equation 3.8: nx « ^Bv(x2) - B(Xl) A (x2 - xx) and the phase ramp for range direction can be calculated from the Equation 2.16. 4ir B2 $r(0djlat,x) ~ -r(^75 + Bn(x) SmOdjlat + Bp(x) COS 6d< flat) (3.9) Here we express the range phase ramp also as a function of azimuth position x to compensate the effects caused by the change of baseline length and its orientation along the track. 51 The above two equations show that an accurate baseline length and its orien-tation must be known in order to thoroughly remove the phase ramps. The linear approximation for the azimuth phase ramp is sufficient in most cases as long as a parallel baseline is accurate enough. The range phase ramp is derived assuming a flat earth. A more realistic approximation should use an ellipsoidal earth model. The ellipsoid is parameterized by its major axis, Remajor, and minor axis, Reminor. The local geoid height, Re, at the geodetic latitude of 9\at is determined by the following formula. Re Re j-i ±L^-manor±L^minor /o i n \ Re = = (3.10) y (Re m a j o r sin 8iat)2 + (Reminor cos 9lat)2 Over the area covered by a SAR image, it is reasonable to further approximate the ellipsoidal earth as a sphere with the above local radius. The flat look angle, 0jiat, refers to the radar look angle to a certain target where its topographic height is zero. Now the look angle for a point at zero elevation on a spherical earth is given by the following equation. Bfiat = arccosf ^ + R ) ], (3.11) where Re denotes the radius of the local sphere, H the satellite height measured from the flat earth's surface, and R the distance from the radar sensor to the target with zero elevation. 52 3.4 Geocoding and Removing Topography Phase After the phase ramps are removed from the interferogram, the remanent phase is still the summation of the phase due to topography and that caused by surface displace-ment. The topographic phase part has to be subtracted before the surface change can be interpreted. Normally the D E M data is in Universal Transverse Mercator (UTM) projection or other standard map projections such as latitude/longitude. An algo-rithm is needed to generate two simulated complex SAR images by illuminating the map from the ERS-1/2 imaging geometry, and at the same time, these two simulated images are precisely registered to the actual interferogram. In this study precise or-bits data in concert with Canadian Topographic Map (1 : 50,000) are used to present the InSAR geometry. Shi et al. [39] of RRSG in U B C have developed an algorithm to convert U T M grid to slant range mapping. The details of the algorithm are beyond the scope of this research. The very basic processing steps of Shi's algorithm are briefed below: The first step of the algorithm is to roughly calculate the maximum and min-imum northings and eastings of the U T M grid for the given slant range frame. Then an U T M grid is denned at giving space in easting and northing. For each point on the vertices of U T M grid, the corresponding geodetic latitude, longitude and earth-centered, earth-fixed (ECEF) coordinates (XQ, yo, Zo) on some reference ellipsoid can be obtained by using equations provide by Snyder [40]. There must be an unique slant range target with the E C E F coordinates of (x,y,z) whose ground projection coordinates are (XQ, yo, ZQ), while both (x,y,z) and (xo,yo,Zo) are 53 located in the same zero-Doppler plane. This zero-Doppler plane is not necessarily one of the samples in the zero-Doppler plane. Normally it is located between two sampled zero-Doppler planes. Thus interpolation should be performed on the sampled zero-Doppler planes. After the interpolation in the azimuth direction, the distance between (x, y, z) and the radar sensor is calculated. Then a search is done to locate (x,y,z) between two sampled range pixels. Interpolation in range direction is performed to obtain the simulated SAR image. Then the topography phase can be derived from this simulated complex image pair by multiplying one image with the complex conjugate of the other. This phase is then removed from the interferogram using the same method. 3.5 Multilook Averaging The phase variation of a single look complex interferogram normally is too large for accurate estimation of ground displacement. To reduce the phase variance, the in-terferogram can be multilook averaged at the expense of resolution. This is realized by replacing groups of adjacent pixels in the interferogram with their averages. Av-eraging the complex interferogram yields better results than directly averaging the phase [30]. In this study 2 pixels in range by 8 pixels in azimuth averaging is used. This yields approximately square multilook pixel due the the resolution in the azimuth 54 direction is 4 times greater than that in the range direction. 3.6 Phase Unwrapping After the phase ramps and topographic phase have been removed, a pure displacement interferogram is created. The phase value of a pixel is a function of the surface change experienced between two consecutive SAR passes. Large variations in the surface displacement will cause the phase value wrap over multiples of 2n. Phase unwrapping, the procedure of resolving these 2n phase ambiguities, is necessary to extracting the displacement information from the interferogram. Due to residues caused by low SNR or layover, phase unwrapping is usually the most difficult step of InSAR processing. The main idea of phase unwrapping is to choose some paths across the fringes, then apply integration along these paths to unwrap the phase. So the success of phase unwrapping mainly depends on choosing these paths. Valero [32] gave a survey of phase unwrapping techniques with applications to InSAR. He analyzed tens of phase unwrapping methods and found that cut-line methods and binary mask methods [41] are most promising. During the processing, residues are identified at first, then branch cuts are defined in conjunction with the identified residues and serve to interdict integration paths. Ghiglia and Romero [42] reported that 2-D phase unwrapping was related to the discrete version of an elliptical partial differential equation. But it was reported 55 that this method had problems dealing with irregular boundaries. Prati et al. [10] introduced the concept of ghost lines which delineate a border between residues, and presented a hybrid method from cut-line method and elliptical partial differential equation method. This method is based on a region partitioning of interferogram. Each of the partitioned regions is unwrapped according to the coherence magnitude of each region. X u and Cumming [43] developed a region growing algorithm which is claimed to be able to handle noisy SAR interferograms more than other other unwrapping algorithms. The region growing algorithm minimizes unwrapping errors by starting at pixels with high coherence magnitude, and proceeding along paths where unwrapping confidence is high. The algorithm is also be able to correct unwrapping errors to some extent, and stop their propagation. Because this algorithm is used in this study to unwrap the pure displacement phase, more detail is given here. • Unwrapping is carried out in several regions, which is started from a seed where data has relatively high coherence magnitude and allowed to grow outwards during the unwrapping procedure. • Each pixel is unwrapped based on prediction made from its neighbors. This allows phase change larger than ir exist between two adjacent pixels. • Information from as many directions as possible is used to unwrap each pixel. This mitigates the effect of errors in the individual prediction directions. • A reliability check is applied when each pixel is unwrapped. 56 • The reliability tolerance is gradually relaxed to allow as many pixels as possible to be unwrapped while still maintain the reliability. • Regions are merged with a reliability check when they overlap, and new regions are initialized whenever qualified data and computer resources are available. Also the structure of this algorithm is suitable for distributed computation. To the areas with very high concentration of residues, it is usually impossible to realize unwrapping using any algorithm. Thus manual processing is needed to handle these areas. 3.7 Converting Phase Information into 3-D Ground Displacement Vectors After using a known spot as zero-speed reference, the interferogram phase is con-verted into slant-range displacement. In order to obtain 3-D surface displacement vector, both measurements from ascending pass and descending pass together with an assumption on the surface target are needed. The conversion and projection from slant-range displacement to the 3-D surface movement vector are given in Chapter 4. 57 3.8 Baseline Estimation For all currently operating SARs, baseline estimated from satellite ephemeris data are not sufficiently accurate. In order to achieve accurate velocity and topography estimates, the baseline must be known to centimeter accuracy or even better. In this section, the accuracy of topography and velocity caused by the baseline error is presented first, then the baseline estimation algorithms are discussed. 3.8.1 Relationship of topography and velocity accuracy with the baseline error The sensitivity of height error to baseline estimation can be derived from Equation 2.9 and 2.10: 8ZB = ^ tan(0 - a) sin 0<J£ (3.12) B This error varies greatly with the tilt angle, a. When a « 8, baseline is perpendicular to the radar look direction (i.e. Bp —>• 0), relatively very small errors occur. While much larger errors occur when (a — 8) « n (i.e. Bn —¥ 0). The ESA regularly publishes the parallel and normal components of the base-line of ERS-1/2. It is more convenient to express the elevation error as a function of these two components. Considering 8d = 8 — 8C and B2 = B2 + B2, the elevation error contributed by the parallel component baseline can be derived from the above equation: xy flcosfljsinfl . . Bp sm 0d - Bn cos 8d 58 Elevation error is related to errors in Bn by i?sin 9d s'mO 8B, (3.14) Bp sin 9d - Bn cos 9d n Figure 3.3 shows 8ZBP as a function of look angle for a parallel baseline com-ponent error, SBP, of 1 rn. In this simulation, Bv is fixed with a value of 100 m, while Bn is either 50, 100 or 150 rn. A similar plot is included in Figure 3.4 showing 8ZB„ as a function of look angle with 8Bn = 1 m. It is clear that the error caused 8Bp and 8Bn varies with the incidence angle. Error in Bp can cause a significant elevation error, although Bp has little effect on the way the interferometer responds to changes in elevation. For Bn the plots of error pass through zero at the center of the image (i.e. 9 = 23°). According to Equation 2.21, the sensitivity of the interferometer to height is proportional to Bn, elevation errors are larger for smaller baselines. Once again, this shows that for the purpose of topography mapping, longer normal baseline is preferred. The error of motion estimation caused by the baseline error is determined by the residual phase ramps. This phase is due to the ground range and topography are removed incorrectly. From Equation 2.16, the phase ramp error can be expressed as: Because the phase error is the form of a ramp, the differential motion error between points increases with the range distance between them. So a higher accuracy can be achieved if the investigated area has small extension in range direction. From the magnitude of the errors in either baseline component, it is clear that 8$ ramp 47T T sm9d8Bn + cos9d8Bp + -{8Bn + 8BP) . (3.15) 59 E l e v a t i o n e r r o r w i t h B = 1 0 0 m , B = 5 0 , 1 0 0 , 1 5 0 m , a s a f u n c t i o n o f 6 n p -2600 -2800 P -4200 23 6 (degrees) Figure 3.3: Relationship of elevation error with parallel baseline component error. the baseline estimates must be extremely accurate. The required accuracy could exceed the accuracy of baseline estimates derived only from the satellite orbit data. As a result, a baseline calibration procedure is necessary by using the known elevation of a certain area. 3.8.2 Baseline estimation There are several well-known algorithms to calibrate the baseline. The most common used algorithm is baseline estimation using tie points (points of known elevation) [41]. Due to the convergence of ERS-1/2 orbits, baseline length and orientation vary along the satellite track. This algorithm assumes a linear baseline length change model for 60 E l e v a t i o n e r r o r w i t h B = 1 0 0 m , B = 5 0 , 1 0 0 , 1 5 0 m , a s a f u n c t i o n o f 0 p n 22 23 24 0 (degrees) Figure 3.4: Relationship of elevation error with normal baseline component error. both Bn and Bp. After establishing the relationship between unwrapped phase with the unknown baseline parameters, at least 4 tie points are needed to solve for the 4 unknown variables: center scene normal baseline component length, B„, the change rate of Bcn along track, center scene parallel baseline component length, B^, and its change rate along track. Zebker et. al. [38] shows that an accuracy of elevation with a standard deviation of 5 meters can be achieved by using around 8 tie points. Of course, the accuracy increases with the number of tie points. In this study, precise (PRC) orbits data are used to convert the topography information into interferogram phase. Reigber et. al. [44] shows that P R C data is accurate up to 10 cm which can yield an accuracy of several millimeters for surface displacement. The details of baseline estimation algorithms are beyond the scope of 61 this research and will not be discussed here. 62 Chapter 4 Glacier M o t i o n Est imat ion 4.1 I n t r o d u c t i o n Glaciers are important components of the hydrological cycle in mountainous areas. They are key indicators of climate change, and also are important natural resources for areas such as western and central Canada. For example, the runoff from glaciers dominates the flow of many western rivers during summer, and provide an important buffering effect on the entire river system by providing more runoff during warm dry summers than wet cool summers. An understanding of the process which could lead to such changes is hindered by the inability to make accurate measurements of the state of the glaciers. For ice mechanics and glacier mass balance studies, the glacier velocity field is always desired. The existing velocity measurement methods are accurate to a few millimeters per observation interval, but are capable of sampling only limited locations at favorable times of the year and the costs of field logistics 63 and personnel are not always affordable. In this chapter, we detail the InSAR technique to realize the estimation of glacier surface displacement at a large scale. C-band ERS1/ERS2 Tandem Mission data with one-day repeat coverage were used to carry out the experiment. In this study, we focus on the Lowell Glacier near the border of British Columbia and the Yukon in Canada (./V60.30, W^ISS.S0) as our research site, because it is a typical mid-latitude glacier. Furthermore, the approach can be easily extended to monitor a large group of alpine glaciers which are impossible to measure by conventional ground-based techniques. Two different approaches are used to convert the radar LOS displacements into the glacier surface 3-D velocity field. The first attempt is using only a single LOS measurement together with two assumptions pertaining to the glacier flow direction. The second approach is combining ascending pass and descending pass LOS mea-surements with only one flow assumption to resolve the 3-D velocity vector. This is followed by an evaluation on which approach and assumption are more suitable to achieve the velocity estimation. 64 4.2 S t u d y A r e a The ERS-1/2 operation is divided into a series of missions with differing coverage and repeat cycles. Between 16 August 1995 and mid-May 1996, ERS-1 and ERS-2 were being operated simultaneously. This is the first time ever that 2 identical SARs were operating in tandem. The orbits of the two satellites were carefully phased to provide a 1-day revisit interval; this enables collection of (interferometric) SAR image pairs which can be used for the generation of a global Digital Terrain Model (DTMs) set. In addition, the tandem mission data can be used in many other novel applications such as differential interferometry for geo-hazard risk assessment arising from earthquakes, volcanic eruptions, landslides and glacial surges. In most cases, glacier surface changes do not cause significant temporal decorrelation during a 1-day interval. This makes it practical to use the interferogram to interpret the surface motion pattern of the glacier. ERS-1/2 coverage is broken up into individual 100 km x 100 km frames. Images are usually processed to cover portions or the entire area of a frame. The location of the ascending pass and descending pass ERS-1/2 frames which cover the Lowell Glacier is shown in Figure 4.1. The black box in the figure represents the extent of the Lowell Glacier. We chose frame 2385 from the descending pass data set, and frame 1215 from the ascending pass data set because they cover the whole glacier and have enough surrounding areas to aid the registration process. Figure 4.2 is a descending-pass SAR magnitude image taken on October 23, 1995 which gives an overview of the Lowell Glacier and its surroundings from the radar line-of-sight 65 perspective, while Figure 4.3 shows the overview of the Lowell Glacier on a 1 : 250,000 Canadian geographic map published in 1978. Each contour line indicates a 200 meter change in elevation. The SAR image has been stretched in the range direction to give it the same aspect ratio as the geographic map. A small tilting angle exists between the orientation of the descending pass SAR image and that of the geographic map. This is because the satellite track has an 16° angle away from the longitude. It is clear that the Lowell Glacier runs towards the north-east with an angle of around 20° in the left side of the image, then bends with almost a right angle, where it momentarily joins with the Dusty glacier. The orientation of the remaining lower part of the Lowell Glacier is nearly eastward, with a heading of about 100°. This means both ascending and descending pass LOS measurements have relatively good viewing directions for the lower part of the glacier. On the right of the scene, the glacier flows into a terminal moraine and lake, and then the Alsek River. The upper portion of the glacier is unfortunately almost perpendicular to the radar viewing direction for both descending pass and ascending passes, which makes it impossible to measure the surface displacement of this part of the glacier. So in this study, we only investigate the feasibility of our approach along the lower portion of the glacier below the sharp bend. 66 Figure 4.1: Latitude/Longitude of the related ERS-1/2 SAR frames. 67 Figure 4.2: SAR magnitude image, Oct. 23,1995, showing an overview of the Lowell Glacier. Azimuth: top-to-bottom; Ground range: right-to-left. The scene covers an area of 23 km x 60 km. Figure 4.3: Canadian geographic map showing an overview of the Lowell Glacier. 68 4.3 ERS-1 /2 Interferograms of the Lowell Glacier Ten tandem mode data pairs over the Lowell Glacier have been collected. First of all, the raw data was processed to generate SLC images using the M D A dtSAR (Desk Top Synthetic Aperture Radar) processor. Then each tandem data pair was co-registered following the processing steps discussed in Chapter 3. The measured coherence magnitude (p) are listed in Table 1. The 22/23-Mar-1996 raw data pair contains some unknown problems which prevents successful registration. In examining the beam geometry and the images, we found that layover is not a big problem on the glacier surface in both the ascending and descending data. In the case of the ascending pass, some parts of the glacier have a relatively low coherence magnitude. This will lower the measurement accuracy by a small amount. Pass Date RO Bn (m) P Pass Type 17/18 Sep 1995 300 -82 0.18 Descending 22/23 Oct 1995* 300 -87 0.73 Descending 31 Dec 95/01 Jan 96* 300 135 0.54 Descending 04/05 Feb 1996* 300 -185 0.51 Descending 29/30 Sep 1995 464 218 0.21 Ascending 08/09 Dec 1995* 464 -80 0.47 Ascending 12/13 Jan 1996* 464 -113 0.68 Ascending 16/17 Feb 1996* 464 -169 0.62 Ascending 22/23 Mar 1996 464 -108 N / A Ascending 26/27 Apr 1996 464 -14 0.27 Ascending Table 4.1: Summary of ERS-1/2 passes over the Lowell Glacier The processing steps to achieve the interferogram from the ERS-1/2 tandem single look complex data are also detailed in Chapter 3. Basically the SLC images 69 are first co-registered to an accuracy of one-tenth of a pixel. To optimize coherence, the images are then filtered in the range and azimuth directions. The range filtering selects only overlapping portions of the object spectrum, which varies with beam elevation angle. The azimuth filtering suppresses non overlapping portions of the azimuth Doppler spectra resulting from processing with different Doppler centroids. These filtering operations improve the scene coherence magnitude at the expense of the spatial resolution. The filtered images are then oversampled by a factor of 2 in both the range and azimuth directions, and the interferogram formed. The coherence is estimated using an averaging window size of 3 (range) by 15 (azimuth) samples and corrected for finite signal-to-noise ratio [4]. Since the azimuth and range filterings are to used avoid decorrelation due to spectral misalignment and the measured coherence has been corrected for the finite signal-to-noise ratio, the resulting coherence estimate can be interpreted as measuring the temporal scene coherence only. From the discussion in Chapter 2, the interferogram coherence magnitude should be high enough to achieve an acceptable estimation accuracy. In this study, we use those interferometric pairs with adequate coherence magnitude (about 0.45 or greater) to continue our experiment. Thus only three ascending pairs and three descending pairs, marked by * in each in the table, are used as our candidates. A representative descending-pass interferogram from the October 22/23, 1995 data pair of the Lowell Glacier is shown in Figure 4.4. In order to achieve accurate registration, enough mountain areas full of features are included in the SLC images. The lower part of the Lowell Glacier runs from the left to the right. A representative 70 ascending-pass interferogram from the January 12/13, 1996 data pair of the Lowell glacier is shown in Figure 4.5. The orientation of the glacier is a quite bit different from the descending orbit due to the different look angles of the radar sensor. In both figures, the top panels show the interferogram intensity and scene coherence magni-tude; the bottom panels show the raw interferogram phase and flat-earth corrected phase respectively. A l l the phase images are naturally wrapped in interval of 27T ra-dians, shown as a fringe. The coherence magnitudes of these two interferograms are quite high (greater than 0.7) on most part of the glacier except at the areas near the toe of the glacier, especially in the water area of the Alsek river. It is also noticeable that the scene coherence is a little bit lower at the upper portions of the glacier where ice accumulation might occur. (c) raw phase (d) flattened phase Figure 4.4: Representative descending pass interferogram from October 22/23, 1995. 71 (a) interferogram magnitude (b) coherence magnitude (c) raw phase (d) flattened phase Figure 4.5: Representative ascending pass interferogram from January 12/13, 1996. 4.4 Estimation of the LOS Displacement The interferograms achieved in the previous section contain phase information due to topography and the glacier surface motion. In order to isolate these mixed phase information, accurate topography information over the test area is needed. There are three different approaches that can be used to achieve the phase separation. • if the motion mode of the glacier is constant and the scene is quite coherent over 3 or more consecutive observations, we can estimate both the topography and the displacement of the observed area by combining the InSAR measurements with different baselines [18]. 72 • if the surface topography and the geometry of the satellite orbits are known, it is possible to convert the surface topography information into interferogram phase, then we can subtract it from the interferogram with mixed phase to get pure motion-induced phase. • obtain an InSAR pair with near-zero baseline, in which case the topogrphy does not induce significant phase as the parallax is small. Normally, most glaciers move at a speed from several tens centimeters to several meters a day. The ERS-1/1 and ERS-2/2 interferometric pairs will dramatically decorrelate with the second shortest repeat interval (35 days) during the period when the tandem mission is in operation. In our case, we use the second approach, as neither zero-baseline data nor data with coherence over 3 passes was available. Then, we need topography information over the test area. Glaciers are poly-crystalline ice masses combined with minor amounts of rock and other impurities. They flow downhill as a visco-plastic-elastic material. The status along the centreline of a glacier is the most important component to revealing the glacier mass balance. So here we will only concentrate on estimating the 3-D velocity along the centreline of the Lowell Glacier. Because no D E M data is available on the Lowell Glacier, we use the 1:50,000 Canada topographic map (115B7-8/115C5) published in 1987. The contour line interval of this map is 40 meters. The elevation along the glacier centreline is read from this map and illustrated in Figure 4.6. It is inevitable that the elevation information contains errors. But if the relative elevation accuracy read from the map is 10 m, only a relative displacement 73 error of 0.27cm is made with a 100 m normal baseline [31]. So we think the map elevation accuracy is adequate for this study. 1 2 0 0 1 1 0 0 1 0 0 0 8 0 0 > CD 6 0 0 5 0 0 — ^ ^+^-— i i i i 1 1 1 ! ! — i J -1 u 1 j — -i _ i i \ v 1 L i j — -I i i i 1 i T 1 ^ * N . r | i T 1 1 r i r ^ 1 1 i i j_ - V i 1 1 1 1 1 i i i i ^: 0 5 1 0 1 5 2 0 2 5 3 0 G l a c i e r c e n t r e l i n e ( k m ) Figure 4.6: Elevation along the glacier centerline. At the toe of the glacier, we find a rock area and consider the displacement to be zero at this point. The published ERS baselines are accurate to within a few meters [45]. This is not accurate enough to remove the effects of the topography which leaves a residual phase ramp. For example, a 1 meter error in Bn introduces a velocity error of approximately 10.7 cm/day. In this study, precise orbit data is used to estimate the normal baseline to an accuracy of 10 cm. The flat-earth corrected phase is first unwrapped using the region growing phase unwrapping algorithm [46]. Then the elevation along the glacier centreline is converted into topographic phase and registered to the interferogram by using the precise orbit data. The phase difference can be calculated by subtracting the 74 topographic phase from the interferogram. Finally, the phase change is converted to radar line of sight displacement. Figures 4.7 and 4.8 show the LOS displacement SR along the centreline of the glacier for the descending and ascending passes respectively. From these two figures, we note the existence of significant differences (about 3 cm to 5 cm/day) in LOS displacement between the ascending and descending results due to the differing viewing angles. Even within descending or ascending orbit set, there is also a discrep-ancy of around 2 cm among various measurements. These nondeterministic phase components could be caused by several different factors: 1. the radar receiver noise contributes phase error 2. atmospheric disturbances can create phase error 3. the registration step in our processing is not perfect 4. the inaccuracy of the precision orbit data 5. error in reading the elevations from the topo map 6. the glacier motion may not be completely stationary over the 4 months of the study. Mattar et al. [47] and Reigber et al. [44] found that the inaccuracy of P R C orbit data will cause a linear phase trend of up to 1 radian, or 0.5 cm LOS, in range and azimuth directions over a 10 km scene. Mattar et al. [47] also showed that the differential interferometric phase can also have a component due to atmospheric 75 inhomogeneities, typically on the order of 1 — 2 radians, or 0.5 — 1.0 cm LOS. How-ever, despite these sources of error, the overall agreement between the various LOS displacement measurements in Figures 4.7 and 4.8 is encouraging. Lowell Glacier LOS Flow Rate - descending orbits 25 20 b >% cC •o E H- 15 cz CD E CD O cc Q . 10 CO CO o - 1 5 22/23 Oct 1995 31 Dec 95/01 Jan 96 04/05 Feb 1996 10 15 20 Glacier centreline (km) Figure 4.7: Lowell Glacier LOS flow rate from descending orbit data. 4.5 Projection to Glacier Flow Direction The LOS displacement must now be projected to an assumed glacier surface flow direction. The basic geometry which relates the LOS displacement and the glacier surface displacement is illustrated in Figure 4.9. whose axes are given by the satellite track vector x, the cross-track or ground range vector y and the local vertical z. The radar LOS direction is shown, along with the measured LOS displacement 7?, which 76 30 Lowell Glacier LOS Flow Rate - ascending orbits 2 5 5^ CO ~0 2 0 o a. CO CO O 12/13 J a n 1996 16 /17 Feb 1996 0 8 / 0 9 Dec 1995 \ \ \ 10 15 20 Glacier centreline (km) 2 5 3 0 Figure 4.8: Lowell Glacier LOS flow rate from ascending orbit data. are assumed to lie in the y-z plane1 with an angle 9 from the vertical. The surface displacement, D, can be derived from the LOS displacement magnitude, R, using the following formula [4]: D R sin LI cos 9 -f sin 7 cos LI sin 6 \ (4.1) where • LI is the angle between the vertical and local surface slope normal, varying between 0° and 90°. 9 is the radar incidence angle, changing from 20° to 26° for ERS. 9 is positive 1Because E R S - 1 and 2 were operated in "yaw-steering" mode, and the data were processed to a near-zero Doppler centroid. 77 for left-looking and negative for right-looking radars. In the case of the Lowell Glacier, 8 is —22.5° for descending orbits and —21.5° for ascending orbits. • 7 is the angle between azimuth and the projection of the local surface normal onto the horizontal, which varies between —180° and 180°. Vertical Radar LOS Z \ Satellite Motion Figure 4.9: The geometry projecting LOS displacement to surface displacement. Once a flow direction assumption is made, the horizontal flow direction 7 is given by: 78 7 = v + 2.2° + 6 - 90° (4.2) where v is the angle from grid east to the horizontal component of the displacement, measured clockwise, the 2.2° covergence angle converts the U T M grid north to true north for the location of the Lowell Glacier, and S is the platform track angle, —16° and 196° for the ascending passes and descending passes respectively. In equation 4.1, it is not hard to find that the denominator can be zero with certain critical geometry angles of fi and v. This happens when the surface dis-placement is perpendicular to the radar LOS, in which case it is impossible for the radar to measure the surface displacement. Figure 4.10 shows the critical angles for the descending and ascending ERS-1/2 passes. Good projection accuracy is achieved when the angles /J, and v for the measured sites are well away from these critical angles. The solid block in Figure 4.10 shows the range of displacement angles for our study area, indicating that our projection angles are favourable. The derivatives of the surface displacement D respect to both fi and 7 are: 3D R cos fi sin 9 cos 7 <?7 (sin LI cos 9 + sin 7 cos LI sin 9)2 dD i ?(cos / i cos# — sin LI sin 7 sin 9) dfi (sin fi cos 9 + sin 7 cos fi sin 9) 2 (4.3) (4.4) where the denominators are the square of the denominator of equation 4.1. As fi and 7 are close to the critical curves, the surface displacement are very sensitive to the accuracies of fi and 7 (i.e. | ^ and | ^ become larger). 79 Critical angles of projecting from LOS to surface displacement -150 -100 -50 0 50 100 Displacement direction: CW from Grid East, v (degrees) 150 Figure 4.10: Plot of the critical angles in the calculation of the surface displacement for the descending and ascending passes of ERS-1/2. There are four unknown variables, R, fi, 9 and 7 in equation 4.1. In order to determine the 3-D surface velocity vectors, three possible approaches can be used. First, if only one InSAR LOS measurement is available, only one degree of freedom, 9 can be resolved in the estimated flow. Thus an assumption on the 3-D glacier flow direction, or the angles of \x and 7 , must be made. Then the magnitude of the surface displacement can be derived from equation 4.1. Second, if two InSAR measurements along different LOS are available, two degrees of freedom can be resolved in the estimated flow. A n assumption on 7 or fj, is needed to determine the 3-D surface displacement vectors. Because it is not easy to know the angle between local surface normal and the vertical, ti, of the glacier, 80 this study focuses on examining different assumptions on the glacier horizontal flow direction, 7 . The third approach relies on having three non-coplanar lines-of-site. The extra InSAR line-of-site measurement could be provided by an airborne SAR or another spaceborne SAR with different viewing angle from ERS-1/2. Then the 3-D glacier surface displacement can be finally determined without any assumption on the glacier flow direction. Of course, all of the above approaches should be based on the premise that the glacier flow rate is constant over the whole interferometric measurement set. In this following part of this section, the first two approaches are used to resolve the 3-D glacier flow rate. 4.5.1 Estimating the 3-D glacier displacement using single-track LOS measurements The glacier local surface does not necessarily lie in the plane of the glacier surface. Normally, the glacier flow direction points slightly downward with respect to the glacier surface at accumulation area, and slightly upwards in ablation area. Fig-ure 4.11 shows the velocity vectors in an idealized glacier in steady-state motion. However, the radar is measuring what happens on the surface during the 1-day obser-vation interval, and if snow falls or blowing snow or significant melting occurs during this interval, it will be noted by a reduction in coherence, possibly to the point of ruining the measurements. 81 The scene coherence magnitude will suffer in regions where significant accu-mulation or ablation takes place. Significant accumulation and ablation between the data sets that form the interferogram would result in a significant change in path length and the surface characteristics, thus a loss of coherence. In this study, all the candidates of interferometric pairs from both descending and ascending passes have quite high coherence magnitude (greater than 0.5) along the whole glacier. It is rea-sonable to assume that both accumulation and ablation could not significantly affect the measurement accuracy, and the glacier flow direction is parallel to its surface. In order to estimate the 3-D glacier displacement using a single LOS mea-surement, two assumptions have to be made. One is the glacier horizontal flow direction, v, the other is the angle between local surface slope normal and the verti-cal, ji. In this section, in addition to the surface-parallel assumption, moraine-aligned and greatest-slope assumptions are used respectively to resolve the values of v and Ll. 82 Under the assumptions of surface-parallel and moraine-aligned for the glacier flow direct ion From the intensity SAR image of the Lowell Glacier, we found that there are several distinguishable moraine lines on the surface of the glacier. These lines show the horizantal direction of the surface displacement, assuming the present flow directions have remained stable for a long time. The longest moraine line is close to the glacier centreline. We consider this line as the horizontal flow direction of the glacier. The medial-moraine line ends at the place around 20 km from the the starting point. In order to compare the results derived from different assumptions on the glacier horizontal flow direction, we further limit our attention on the velocity vectors only along this 20 km centreline in the rest of this study. The 1:50,000 Canadian Topographic Map of the Lowell Glacier was used to measure the values of /J, and v. To reduce the root-mean-square (RMS) errors, five measurements along the medial-moraine line have been made and then averaged. The measured values of /i and v are shown in Figure 4.12 and Figure 4.13. Using these values together with equation 4.1 and 4.2, the surface displacement is estimated in turn using the descending and ascending LOS measurements. Because the combined values of \x and v are well located outside the critical curves of both descending and ascending passes, it is expected both projections can give reasonable results if the flow direction assumptions are accurate. Each set of descending and ascending LOS measurements is averaged first to continue the projection calculation. The derived surface displacements over the 1-day interval which separates the tandem 83 2.5 V e r t i c a l F l o w D i r e c t i o n o I I I I I I I I I I I 0 2 4 6 8 10 12 14 16 18 20 G l a c i e r c e n t r e l i n e ( k m ) Figure 4.12: The vertical flow direction of the Lowell Glacier based on the surface-parallel and moraine-aligned assumptions. data pair along the medial-moraine line are plotted in Figure 4.14. The results from both descending and ascending LOS measurements agree with each other quite well except for a large discrepancy at the starting point, implying that the flow direction assumption does not hold at this region. The glacier velocity reaches the maximum around 65 cm j day at the starting point region, then decreases almost monotonically to 35 cm/day at the 20 km point of the medial-moraine line. Under the assumptions of surface-parallel and along the greatest-slope for the glacier flow direction For the surging type glaciers like Lowell, it is also reasonable to assume the glacier flow direction is surface-parallel and along the greatest downhill slope. Based on 84 H o r i z o n t a l F l o w D i r e c t i o n 25. 1 1 1 1 1 1 0 2 4 6 8 10 12 14 16 18 20 G l a c i e r c e n t r e l i n e ( k m ) Figure 4.13: The horizontal flow direction of the Lowell Glacier based on the surface-paralle and moraine-aligned assumptions. these two assumptions, the vertical angle, LI , and horizontal angle, v , of the glacier flow direction can also be measured from the Canadian Topographic Map over the Lowell Glacier. The value of LI is plotted in Figure 4.15, and the value of v is shown in Figure 4.16. Using the same projection equation as the previous section, the corresponding surface displacement over the 1-day interval can also be derived. Figure 4.17 shows the results from descending pass and ascending pass data separately. Apparently, a significant discrepancy exists between the results from descend-ing pass and ascending data around the sharp bending area of the Lowell Glacier. It is partially because the composed values of JJ, and v at this region make the the actual surface displacement almost perpendicular to the radar LOS for the descending track. This is clearly shown in Figure 4.10. The right side of the gray box is very close to the 85 Surface Flow Rate Derived from One Type Orbit Data (A) 80 CO . c CD E CD O C L C/3 T3 CD O CO tz co 1 1 1 1 1 ~~ " ~ \ \ • — -. . , \ . -- : — - - • . . . -— ^ \ \ . \ \ 1 i I I N 1 6 8 10 12 14 Glacier centreline (km) 16 18 20 Figure 4.14: The Lowell Glacier surface velocity along the glacier centreline using single-track LOS measurements, under the surface-parallel and moraine-aligned as-sumptions. critical curve for descending pass. In this case, we should mainly rely on the result derived from ascending pass data. The other explanation is the flow assumptions are probably not correct for this region. 4.5.2 Estimating the 3-D glacier displacement using dual-track LOS measurements With both descending and ascending LOS measurements, only one flow direction assumption is needed to resolve the 3-D displacement if the glacier flow rate is constant over the period between the two LOS measurements. The following three different 86 Vertical Flow Direction 2.5 I ; — i 1 1 r 0 2 4 6 8 10 12 14 16 18 20 Glacier centreline (km) Figure 4.15: The vertical flow direction of the Lowell Glacier based on the surface-parallel and greatest-slope assumptions. assumptions can be taken into consideration: • The glacier flow direction is aligned with the medial-moraine line; • The glacier flow direction is parallel to its surface; • The glacier flow direction is down the greatest slope. Under the first or the third assumption, the glacier horizontal flow direction, v can be measured from the Canadian topographic map. We use the same testing points as in the single type LOS projection attempt. Thus we are able to take advantage of using the same values of v. After obtaining the value of u , we can use equation 4.1 twice. Then the surface displacement magnitude and the value of fi can be calculated. 87 H o r i z o n t a l F l o w D i r e c t i o n 6 8 10 12 14 16 18 20 G l a c i e r c e n t r e l i n e ( k m ) Figure 4.16: The horizontal flow direction of the Lowell Glacier based on the surface-parallel and greatest-slope assumptions. If the glacier flow direction is parallel to its surface zs(x,y), we separate the surface velocity vector into three components Vv = Vxx + Vyy + Vzz. (4.5) Then the vertical velocity, Vz , can be related to the horizontal velocity, Vx and Vy by [48]: d d Vz = Vx—zs(x,y) + Vy—zs(x,y) (4.6) By using the above equation and the projection equation, The relationship of Vx, Vy with the LOS displacements from descending pass, Rd, and ascending pass, Ra, can be established [48]. We can then measure the values of v , fi and the surface 88 Surface Flow Rate Derived from One Type Orbit Data (B) 0 2 4 6 8 10 12 14 Glacier centreline (km) Figure 4.17: The Lowell Glacier surface velocity along the glacier centreline using single-track LOS measurements, under the surface-parallel and greatest-slope assump-tions. displacement magnitude. The measured values of v, / i and the surface displacement from the above three assumptions are illustrated in Figure 4.18, Figure 4.19 and Figure 4.20 respec-tively. 89 50 CD CD 40 D) CD "D Horizontal Flow Directions Using Different Assumptions 30 I 20 o CD 5 , w o 0 -*—• c o N •s -10 x -20 / / / \ \ - / • • / / / - - . \ \ \ / - / ' 1 /' ' /' \ • A \ \ / / / 1 ^ \ X r Jj \ V ' \ ' / 1 ' 1 / 2 4 6 8 10 12 14 16 18 20 Glacier centreline (km) Figure 4.18: The horizontal flow direction of the Lowell Glacier derived from dual-track L O S measurements. 4 3.5 3 CD CD 2.5 cn CD S 2 CD g.1.5 Vertical Flow Directions Using Different Assumptions cn 1 0.5 0 moraine-aligned • — • — • surface-parallel — — — greatest-slope A / \ ' \ \ \ \ \ . . \. . . . .1 / 1 '• < /. .... ' V \ \ \ X/: \ \ / \ / / N / : V v A \ \ \ / • ^ v / y / / / \ \ ~ "r' / \ - -~ 1 i / 2 4 6 8 10 12 14 16 18 20 Glacier centreline (km) Figure 4.19: The vertical flow direction of the Lowell Glacier derived from dual-track L O S measurements. 90 Lowell Glacier Surface Flow Rate Derived from Dual Orbit Data 100 6 8 10 12 14 Glacier centreline (km) Figure 4.20: Lowell Glacier 3-D flow rate derived from dual-track LOS measurements. 91 4.6 Discussion Two different approaches have been used to convert the L O S displacement to the glacier surface displacement. In both attempts, we depend on the Canadian to-pographic map to obtain the information related to the glacier flow direction after certain assumptions related to the flow direction are made. In the attempt to resolve the glacier 3-D flow rate using only single L O S mea-surement, we examined two assumption combinations. Under the flow assumption combination of surface-parallel and moraine-aligned, the results achieved from de-scending and ascending L O S displacements generally agree wi th each other except a significant discrepancy exists around the starting point. Because the values of fi and v are well away from both descending and ascending passes cri t ical values, this discrepancy suggests that at least one of the values of JJL and v does not reflect the real flow direction of the glacier in this region. The results derived from the assumptions of surface-parallel and greatest-slope show a quite significant difference at the bending portion of the glacier between ascending and descending pass projections. As explained in the previous section, the glacier flow orientation, v and \x , is not favourable for the descending pass projection at this region. Even though we only reply on the result from ascending pass data to reveal the glacier flow rate, it st i l l exagerates the flow rate by around 7 cm at this region. In the second attempt, all three assumptions surprisingly give a very simi-lar velocity magnitude at the starting point. But a further investigation into the 92 other product of the estimation, fi, and the coherence magnitude of the interfero-gram, a conclusion can be drawn on which assumption is the best. First of all, the moraine-aligned assumption gives an inconsistent value of / i . Thus, this assumption is not proper in this region. On the other hand, from the descending and ascending pass coherence magnitude images in Figure 4.4 and 4.5, we find that the coherence magnitude is slightly lower at this location. This is probably caused by a moderate accumulation process when extra ice mass comes down from the Dusty glacier at the joint. The value of p = 2.1° derived from the surface-parallel assumption is a little bit larger than the average normal direction angle (0.96) we made along the glacier moraine-line. This coincidently shows that the glacier flow direction is slightly downwards at this area. So the greatest-slope is the most suitable assumption at this area. Another lower coherence area is found covering the range from 12 km to 17 km. This might be a very complicated flow mode combining ablation, accumulation and deformation here. Using only the measured velocity magnitude and the vertical flow angle, it is impossible to explain the glacier mechanism. A fairly good overall agreement has been achieved for the velocity vectors from all three assumptions along the lower portion the glacier. The down-slope assump-tion gives an exaggerated velocity magnitude at the region just below the starting point. This is because the glacier is pushed by its momentum along the centreline of the glacier other than in the direction of the greatest downslope gradient. So each assumption has an advantage and disadvantage. In general, the moraine-aligned as-sumption and the surface-parallel assumption appear to be good for the investigated 93 20 km region. 94 4.7 Conclusions 3-dimensional velocity vectors are important information to glaciological research, but collection of such data with traditional methods is constrained both logistically and technically. Recent research shows that satellite SAR interferometry can pro-vide a promising new method to determine the velocity field of glaciers due to its global-wide coverage and all-weather capability [3, 18, 48]. Although the accuracy of this technique is affected by many different factors, a measured accuracy of sub-centimeter was demonstrated by Vachon et al. in the application of the Athabasca and Saskatchewan Glaciers [3]. In this dissertation, a general introduction of SAR and InSAR theory is pre-sented. Some simulation results are given to illustrate InSAR phase statistics. This is followed by the discussion of the overall InSAR processing procedure. Some of the major InSAR processing steps and corresponding algorithms are examined in detail. A 3-dimensional geometry module is used to derive glacier surface displacement (D) as a function of LOS displacement (R) and the orientation of the glacier flow direc-tion (/i, 7 ) . Critical curves, which show surface displacement is perpendicular to the radar LOS, for both ascending and descending passes of ERS-1/2 are plotted. In order to maintain an acceptable level of projection accuracy, the combination of \i and 7 should not be in the vicinity of these critical curves. In our investigation of the Lowell Glacier, we found these angles are well away from these critical curves except for the case under the assumptions of surface-parallel and greatest-slope. This makes it feasible to measure the surface motion vectors of the glacier using either track data. 95 After a brief discussion on the flow mechanism of glacier, we think the following three assumptions to the flow direction of the glacier are reasonable. 1. The glacier flow direction is aligned with the medial-moraine line; 2. The glacier flow direction is parallel to its surface; 3. The glacier flow direction is down the greatest slope. Based on these assumptions, we were able to examine all the possible projec-tion approaches using the ascending/descending pass data. In the attempt of using a single-track direction data, two combinations of assumptions have been used to resolve LI and 7: the first combination is surface-parallel and moraine-aligned assumptions; and the second combination is surface-parallel and greatest slope assumptions. Obvi-ous discrepances exist between the results derived from these two assumption combi-nations. This shows different flow assumptions give different interpretation of radar LOS displacement. Surprisingly, based on the same combination of assumptions, the results from ascending pass and descending pass data are quite consistent except for the area near the starting point. This is probably because interferograms acquired along a single track are only sensitive to the displacement along one direction. The other explanation for this is the flow assumptions are not proper at this area. Under the assumptions of surface-parallel and greatest-slope, a poor alignment between the radar LOS and the glacier flow direction was formed for descending pass measure-ments. Thus we should mainly use the ascending pass measurement to interpret the magnitude of the glacier motion. 96 Combining the ascending and descending track data and using only one of the above three assumptions, we also have measured the surface displacement of the Lowell Glacier. The results show that each assumption has its advantage and disad-vantage. In general, moraine-algined assumption and the surface parallel assumption appear to be good for the investigated 20 km region, especially for the region just below the starting point. Our results show ERS-1/2 tandem phase SAR interferometric image pairs col-lected with a 1-day repeat-pass interval can be used to measure the surface dis-placment field of a mid-latitude alpine glacier using differential InSAR techniques. Although the LOS displacements change over the time span from October 1995 to February 1996, there is no evidence of a seasonal change in glacier flow rate. In the future work, more representative glaciers with accurate velocity metrics should be chosen. The comparison between the differential InSAR results and the known velocity metrics can help us to understand how accurate this technique can be. It also can provide valuable information to reveal the flow mechanism for a certain type of glaciers. The ERS-1/2 SARs cannot image polar areas of the earch. R A D A R S A T will image all of Antarctica during the R A D A R S A T Antarctic Mapping Mission, but the 24-day time interval is too long for velocity mapping. So in order to be able to monitor global ice-sheets and glaciers, special missions are needed to cover the whole Antarctica in the future. 97 B i b l i o g r a p h y [1] J. Jacob and S. N . Madsen, "Multi-pass Interferometry for Studies of Glacier Dynamics," in ESA Workshop on Applications of ERS SAR Interferome-try, FRINGE'96, (Zurich, Switzerland), September 30 - October 2, 1996. http://www.geo.unizh.cn/rsl/fringe96/. [2] K . H . Thiel and X . Wu, "The Use of Tandem Data in the Antarc-tic Area," in ESA Workshop on Applications of ERS SAR Interferome-try, FRINGE'96, (Zurich, Switzerland), September 30 - October 2, 1996. http://www.geo.unizh.cn/rsl/fringe96/. [3] P. W. Vachon, D. Geudtner, K . Mattar, A. L. Gray, M . Brugman, and I. G. Cumming, "Differential SAR Interferometry Measurements of Athabasca and Saskatchewan Glacier Flow Rate," Can. J. Remote Sensing, vol. 22, pp. 287-296, Sept. 1996. [4] P. W. Vachon, D. Geudtner, A. L. Gray, K . Mattar, M . Brugman, I. G. Cumming, and J . L. Valero, "Airborne and Spaceborne SAR Interferometry: Application to the Athabasca Glacier Area," Proc. IGARSS'96, pp. 2255-2257, May 1996. [5] G. Ostrem and M . Brugman, "Glacier Mass-Balance Measurements," tech. rep., NHRI, Saskatoon, Saskatchewan, Canada, 1991. [6] M . F. Meier, "Contribution of Small Glaciers to Global Sea Level," Science, vol. 226, pp. 1418-1421, 1984. [7] W. S. B. Paterson, The Physics of Glaciers. Elsevier Science Ltd., 1994. [8] L. C. Graham, "Satellite Interferometer Radar for Topographic Mapping," Proc. IEEE, vol. 62, pp. 763-768, 1974. 98 [9] C. Prati and F. Rocca, "3-D Synthetic Aperture Radar Surveys," SEP, vol. 57, pp. 463-477, 1988. [10] C. Prati and F. Rocca, "Limits to the Resolution of Elevation Maps from Stereo SAR Images," Int. J. on Remote Sensing, vol. 11, pp. 2215-2235, 1990. [11] E. Rodriguez and J. Martin, "Theory and Design of Interferometric Synthetic Aperture Radars," IEE Proceedmgs-F, vol. 139, pp. 147-159, 1992. [12] M . D. Thompson and J. B. Mercer, "Digital Terrain Models," EOM, pp. 22-26, Mar. 1996. [13] H. A . Zebker and R. M . Goldstein, "Topographic Mapping from Interferometric Synthetic Aperture Radar Observations," J. of Geophysical Research, vol. 91, pp. 4993-4999, 1986. [14] A. K . Gabriel, R. M . Goldstein, and H. A . Zebker, "Mapping Small Elevation Changes Over Large Areas: Differential Radar Interferometry," J. Geophysical Research, vol. 94, pp. 9183-9191, 1989. [15] 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, Jul. 1993. [16] R. M . Goldstein, H. Engelhardt, B. Kamb, and R. M . Frolich, "Satellite Radar Interferometry for Monitoring Ice Sheet Motion: Application to an Antarctic Ice Stream," Science, vol. 262, pp. 1525-1530, Dec. 1993. [17] R. Kwok and M . A. Fahnestock, "Ice Sheet Motion and Topography from Radar Interferometry," IEEE Trans, on Geoscience and Remote Sensing, vol. 34, pp. 189-200, 1996. [18] I. Joughin, D. Winebrenner, M . Fahnestock, R. Kwok, and W. Krabill , "Mea-surement of Ice Sheet Topography Using Satellite Radar Interferometry," J. of Glaciology, vol. 42, pp. 10-22, 1996. [19] A . L. Gray, K . E. Mattar, and M . W . A. van der Kooij, "Cross-track and Along-Track Airborne Interferometric SAR at CCRS," in 17th Can. Symp. Rem. Sens., (Saskatoon, Sask.), pp. 232-237, Jun. 1995. 99 [20] T. Ngo and I. Cumming, "SAR along-Track interferometry for enhanced vessel wake detectoion," tech. rep., Radar Remote Sensing Group, Dept. of E. E. , U B C , Vancouver, British Columbia, Canada, 1996. [21] 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. [22] C. Wu, K . Y . Liu, and M . Y . Jin, "Modeling and a Correlation Algorithm for Spaceborne SAR Signals," IEEE Trans. Aerospace and Electronics Systems, vol. 18, pp. 563-574, Sept. 1982. [23] H. Runge and R. Bamler, "A Novel High-precision SAR Focusing Algorithm Based on Chirp Scaling," Proc. IGARSS'92, pp. 372-375, 1992. [24] M . Sack, M . R. Ito, and I. G. Cumming, "Application of Efficient Linear F M Matched Filter Algorithms to SAR Processing," IEE Proceedings-F, vol. 132, pp. 45-57, 1985. [25] A . K . Gabriel and R. M . Goldstein, "Crossed Orbit Interferometry: Theory and Experimental Results from SIR-B," International Journal of Remote Sensing, vol. 9, pp. 857-872, 1988. [26] D. Massonnet and T. Rabaute, "Radar Interferometry: Limits and Potential," IEEE Trans, on Geoscience and Remote Sensing, vol. 31, pp. 455-464, Mar. 1993. [27] R. Bamler and D. Just, " Phase Statistics and Decorrelation in SAR Interfero-grams," Proc. IGARSS'93, pp. 980-984, Aug. 1993. [28] H . A . Zebker, T. G. Farr, R. P. Salazar, and T. H. Dixon, "Mapping the World's Topography Using Radar Interferometry: The T O P S A T Mission," Proc. of IEEE, vol. 82, pp. 1774-1786, Dec. 1994. [29] H. A . Z. J . Villasenor, "Deccorelation in Interferometric Radar Echoes," IEEE Trans. Geosci. and Remote Sensing, vol. 30, pp. 950-959, Sept. 1992. [30] I. Joughin, D. Winebrenner, and D. B. Percival, "Probability Density Functions for Multilook Polarimetric Signatures," IEEE Trans, on Geosci. and Remote Sensing, vol. 32, pp. 562-574, May 1994. 100 [31] I. G. Cumming and J. X . Zhang, "Estimating the Flow Rate of the Lowell Glacier Using ERS Tandem Mode Data," in Geomatics in the Era of RADARSAT, GER'97, (Ottawa), May 1997. On C D - R O M . [32] J. L. Valero, "A Survey of Phase Unwrapping Techniques, with Application to InSAR," tech. rep., University of British Columbia, Vancouver, B C , Canada, 1995. [33] I. Joughin, D. Winebrenner, and D. B. Percival, "Observations of Ice-sheet Mo-tion in Greenland Using Satellite Radar Interferometry," Geophysical Research Letters, pp. 562-574, May 1995. [34] J. 0 . Hagberg, L. Ulander, and J. Askne, "Repeat-Pass SAR Interferometry Over Forested Terrain," IEEE Trans, on Geosci. and Remote Sensing, vol. 33, pp. 331-340, Mar. 1995. [35] I. Joughin, "InSAR Applications in Ice-Sheet Motion Estimation," Glaciological Applications of Satellite Radar Interferometry Meeting, pp. T1-L19, March 28-29 1996. [36] J . Lee, K . W. Hoppel, S. A . Mango, and A. R. Miller, "Intensity and Phase Statistics of Multilook Polarimetric and Interferometric SAR Images," IEEE Trans, on Geosci. and Remote Sensing, vol. 32, pp. 1017-1027, Sept. 1994. [37] A . Lopes, "Phase Difference Statistics Related to Sensor and Forest Parameters," Proc. IGARSS'92, pp. 779-781, Aug. 1992. [38] H. A . Zebker, C. L. Werner, P. A . Rosen, and S. Hensley, "Accuracy of Topo-graphic Maps Derived from ERS-1 Interferometric Radar," IEEE Trans. Geosci. and Remote Sensing, vol. 32, pp. 823-836, Jul. 1994. [39] M . S. Seymour, I. G. Cumming, and L. Gray, "The Generation of Geocoded Map-Like Products from CCRS Airborne InSAR Data," Personal Communication, RRSG, UBC, pp. 1-9, 1994. [40] J . P. Synder, "Map Projections Used by the U.S. Geological Survey," U.S. G.P.O., pp. 53-73, 1983. [41] R. M . Goldstein, H. A. Zebker, and C. L. Werner, "Satellite Radar Interferome-try: Two-dimensional Phase Unwrapping," Radio Science, vol. 23, pp. 713-720, 1988. 101 [42] D. G. Ghiglia and L. A . Romero, "Robust Two-dimensional Weighted and Un-weighted Phase Unwrapping That Uses Fast Transforms and Iterative Methods," J. Opt. Soc. Am. A, vol. 11, pp. 107-117, 1994. [43] W. X u and I. G. Cumming, "Region Growing Algorithm for InSAR Phase Un-wrapping," Proc. IGARSS'96, pp. 2044-2046, May 1996. [44] C. Reigber, Y . Xia , H. Kaufmann, F. H. Massmann, L. Timmen, and a. M . F. J . Bodechtel, "Impact of Precise Orbits on SAR Interferometry," in ESA Work-shop on Applications of ERS SAR Interferometry, FRINGE'96, (Zurich, Switzer-land), September 30 - October 2, 1996. http://www.geo.unizh.ch/rsl/fringe96/. [45] G. A . Solaas and S. N . Coulson, "ERS-1 SAR Interferometric Baseline Algorithm Verification," ES-TN-DPE-OM-GS02, ESA/ESRIN, vol. v. 1.1, Mar. 1994. [46] W. X u and I. G. Cumming, "A Region-Growing Phase Unwrapping Algorithm," IEEE Trans. Geoscience and Remote Sensing, Accepted for publication, January 1998. [47] K . E. Mattar, P. W. Vachon, D. Geudtner, A . L. Gray, I. G. Cumming, and M . Brugman, "Validation of ERS Tandem Mission SAR Measurements of Alpine Glacier Velocity," IEEE Trans, on Geoscience and Remote Sensing, vol. 36, pp. X 1 - X 2 , May 1998. [48] I. Joughin, R. Kwok, and M . A. Fahnestock, "Interferometric Estimation of Three-Dimensional Ice-Flow Using Ascending and Descending Passes," IEEE Trans, on Geosci. and Remote Sensing, vol. 36, pp. 25-37, Jan. 1998. 102
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Estimating the 3-D flow rates of the Lowell Glacier...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Estimating the 3-D flow rates of the Lowell Glacier using spaceborn InSAR Zhang, Xiangzhou Joe 1998
pdf
Page Metadata
Item Metadata
Title | Estimating the 3-D flow rates of the Lowell Glacier using spaceborn InSAR |
Creator |
Zhang, Xiangzhou Joe |
Date Issued | 1998 |
Description | Glacier motion data is very important for glaciology research, but traditional methods for collecting such data need to overcome staggering financial and logistical hurdles. During the period from July 1995 to May 1996, ERS-1/2 operated under Tandem Mission Mode with one satellite following the other one day apart. This gives a great opportunity to using Synthetic Aperture Radar (SAR) interferometry (InSAR) technology to measure the relative large displacement on the ground between the observation pair. Pioneer research work shows that this short time interval can maintain relatively high coherence for most glacier movement investigations. Some good results have been achieved to measure the velocity field of ice fields, ice streams and Alpine glaciers [1, 2, 3, 4]. In this study, the rationale and processing scheme of differential InSAR are discussed and examined first. Then 10 Tandem raw data sets (4 descending passes, 6 ascending passes) over the Lowell glacier at the boundary of B.C. and the Yukon, Canada, have been processed. 6 interferograms with enough coherence magnitude (3 ascending passes, 3 descending passes) have been generated. In concert with three different flow assumptions on the glacier movement direction, the line of sight (LOS) displacements along the glacier centreline measured from the ascending orbit and the descending orbit are converted into the 3-D velocity vectors. Our results show that the surface parallel assumption is more suitable for most part of the glacier. An overall accuracy around 4 cm/day rms has been achieved in measuring the surface motion of the Lowell Glacier. However, the 3-D projection model does not apply for the glaciers at all locations on this planet. The glacier has to be away from the critical regions caused by the projection geometry. Also, in order to achieve high estimation accuracy, the flow direction of the candidate glacier has be approximately aligned with the radar LOS. But with specifical satellite missions, this technique can provide a feasible method for monitoring the global glaciers and ice sheets in the future. |
Extent | 6013534 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-06-15 |
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. |
IsShownAt | 10.14288/1.0065028 |
URI | http://hdl.handle.net/2429/9202 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1998-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1999-0308.pdf [ 5.73MB ]
- Metadata
- JSON: 831-1.0065028.json
- JSON-LD: 831-1.0065028-ld.json
- RDF/XML (Pretty): 831-1.0065028-rdf.xml
- RDF/JSON: 831-1.0065028-rdf.json
- Turtle: 831-1.0065028-turtle.txt
- N-Triples: 831-1.0065028-rdf-ntriples.txt
- Original Record: 831-1.0065028-source.json
- Full Text
- 831-1.0065028-fulltext.txt
- Citation
- 831-1.0065028.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-0065028/manifest