Local Fringe Frequency Estimation in Synthetic Aperture Radar Interferograms Using a Multiband Pre-Filtering Approach by Diego E. Perea-Vega B.Sc. in Electronic Eng. Pontificia Universidad Javeriana, 1991 B.Sc. in Physics Universidad Nacional de Colombia, 1996 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR THE D E G R E E OF MASTER OF APPLIED SCIENCE in THE F A C U L T Y OF G R A D U A T E STUDIES D E P A R T M E N T OF E L E C T R I C A L A N D COMPUTER ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH C O L U M B I A May 18th 2000 © Diego E. Perea-Vega In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of ^ ' ^ / W C o ^ £ - \ y \ e c The University of British Columbia Vancouver, Canada Date _ "*i DE-6 (2/88) Abstract Synthetic aperture radar (SAR) interferometry is a technique for obtaining accurate elevation maps of the Earth from radar images. Local fringe frequency estimates are needed in several stages of the interferometry process. They are used to correct the effect of the topographic slope in the estimation of the interferogram coherence. They are needed to define the frequency center for adaptive band-pass filtering or to model the local phase in slope-correction filtering. Finally, accurate fringe frequency estimates facilitate the phase unwrapping process. In this work, I propose a new algorithm for local fringe frequency estimation in which the SAR interferogram signal is pre-filtered before the local frequency estimation is performed. This allows the use of a simpler and more efficient frequency estimator that operates at the pixel level. The proposed scheme shows advantages over other schemes because it achieves a better space-frequency resolution and therefore tracks the topographic changes of the scene more accurately. The filters used in this work are modulated Gaussian functions with variable spatial aperture and bandwidth. In this way, the analysis window is adapted to the local characteristics of the signal at all samples. The variable-aperture filters are similar to the variable space-frequency domain filters used in wavelet analysis. I present results for synthetic and real SAR interferograms, as well as the performance of the proposed algorithm. An application of the frequency estimation method is developed for the noise filtering stage. A non-linear phase model is built to locally flatten the phase allowing the averaging of a higher number of samples without significant distortion. Simulations show that this alternative method achieves a better performance than two other reported methods when the interferogram coherence is moderately high. However, the alternative method does not solve the problem of phase discontinuities in the presence of topographically induced residues. Table of Contents A B S T R A C T II T A B L E OF CONTENTS III LIST OF T A B L E S VI LIST OF FIGURES VII LIST OF A C R O N Y M S IX A C K N O W L E D G E M E N T S X C H A P T E R 1 1 INTRODUCTION 1 1.1 THESIS CONTRIBUTION 2 1.2 ORGANIZATION OF THE THESIS 3 C H A P T E R 2 5 S A R INTERFEROMETRY BACKGROUND 5 2.1 S A R INTERFEROMETRY PRINCIPLES 5 2.2 INTERFEROGRAM PROCESSING CHAIN 11 2.3 T H E U S E OF FREQUENCY ESTIMATES IN COHERENCE ESTIMATION 12 2.4 PHASE NOISE FILTERING APPROACHES 15 2.4.1 Constant Phase Estimation Approach 15 2.4.2 Frequency Estimation Approach 18 2.4.3 Post-filtering Approach 19 2.5 MAIN APPROACHES TO PHASE UNWRAPPING 20 C H A P T E R 3 22 L O C A L FRINGE F R E Q U E N C Y ESTIMATION - EXISTING METHODS 22 3.1 SHORT TIME FOURIER TRANSFORM 23 3.1.1 The S T F T as a Signal Decomposition Transform 24 3.1.2 Use of the Periodogram in Loca l Frequency Estimation 28 3.1.3 Reported Works Using this Approach 30 3.2 Two DIMENSIONAL M U S I C ALGORITHM 31 3.3 MULTI-GRID METHOD 36 3.4 SUMMARY AND CONCLUSIONS 38 i i i C H A P T E R 4 41 L O C A L FRINGE F R E Q U E N C Y ESTIMATION - PROPOSED M E T H O D 41 4.1 A M - F M SIGNAL REPRESENTATION M O D E L 42 4.2 PROPOSED FRINGE FREQUENCY ESTIMATION METHOD 44 4.3 FILTERS PARAMETERS SELECTION 53 4.3.1 Wavelet-l ike Filters Scheme 53 4.3.2 Advantages of the Wavelet-l ike Filters Scheme. 55 4.3.3 Filters Parameters Selection Procedure 60 4.3.4 Appl icat ion of the Wavelet-l ike Filter Scheme 60 4.4 APPLICATION OF THE PROPOSED METHOD TO 2D INTERFEROGRAM SIGNALS 61 4.4.1 Appl icat ion to Flattened Interferograms 62 4.4.2 Practical Procedure 63 C H A P T E R 5 _69 INSTANTANEOUS F R E Q U E N C Y ESTIMATION 69 5.1 INSTANTANEOUS FREQUENCY ESTIMATION 69 5.2 ENERGY SEPARATION ALGORITHM 72 5.3 D E S A - 1 ALGORITHM PERFORMANCE 76 C H A P T E R 6 78 RESULTS ON SYNTHETIC INTERFEROGRAMS 78 6.1 FREQUENCY ESTIMATION RESULTS 78 6.1.1 IF Estimation Error Performance 82 6.1.2 Comparison Against Window-based Methods. 86 6.2 PHASE NOISE FILTERING RESULTS 92 6.2.1 Comparison between Phase Noise Filtering Methods 93 6.3 CONCLUSIONS 95 C H A P T E R 7 J 0 3 RESULTS ON R E A L INTERFEROGRAMS 103 7.1 RESULTS ON E R S TANDEM MISSION D A T A 103 7.1.1 Limitat ion of Noise Phase Filtering Method in the Presence of Residues. 112 7.2 RESULTS ON RADARSAT DATA 119 C H A P T E R 8 131 CONCLUSIONS 131 8.1 SUMMARY 131 8.2 CONCLUSIONS 132 8.3 PRACTICAL APPLICATION OF THE PROPOSED METHOD 133 8.4 FUTURE WORK 133 iv A P P E N D I X A 136 I N T E R F E R O G R A M R E S I D U E S 136 A . l DEFINITION 136 A . 2 D E T E C T I O N A L G O R I T H M 137 A P P E N D I X B 140 R E L A T I O N O F T H E P R O P O S E D M E T H O D A N D T H E S T F T 140 A P P E N D I X C 142 E V A L U A T I O N O F T H E D E S A - 1 I F E S T I M A T I O N A L G O R I T H M 142 C . l F R E Q U E N C Y DISTORTION O F T H E D E S A - 1 E S T I M A T O R 142 C . 2 DISTORTION O F D E S A - 1 WITH F A S T I F C H A N G E S 144 C .3 N O I S E R E D U C T I O N B Y PRE-FILTERING 145 C . 4 D E S A - 1 P E R F O R M A N C E A F T E R M U L T I - B A N D F I L T E R I N G 152 C . 5 S U M M A R Y A N D C O N C L U S I O N S 155 A P P E N D I X D 159 P H A S E N O I S E F A L T E R I N G M E T H O D S 159 D.1 Standard Mul t i -Look Averaging 159 D . 2 Slope-Correction Filter 160 D . 3 Non-Linear Phase Model ing Filter (Proposed Method) 165 R E F E R E N C E S 171 List of Tables Table 2.1 Example of interferogram acquisition parameters 8 Table 4.1 Filters Parameters of Figure 4-6 56 Table 4.2 Filters Parameters of Figure 4-12 64 Table 7.1 ERS Interferogram Processing parameters 105 Table 7.2 ERS Interferogram Acquisition parameters 106 Table 7.3 Radarsat Acquisition parameters 120 Table 7.4 Radarsat Processing parameters 121 vi List of Figures Figure 2-1: Interferogram acquisition geometry 8 Figure 2-2: Magnitude of ERS Tandem mission interferogram 9 Figure 2-3: Phase of ERS Tandem mission interferogram 10 Figure 2-4: Interferogram Processing Chain 13 Figure 2-5: Interferogram average magnitude spectrum in the range direction 14 Figure 2-6: ID slice of the phase in the transect line A of Figure 2-3 14 Figure 3-1: One-dimensional Synthetic interferogram signal 26 Figure 3-2: Space -Frequency representation of signal in Figure 3-l(b) 27 Figure 3-3: ML Frequency estimator based on the local periodogram 33 Figure 3-4: Adaptive filtering using 32 pixels FFT 34 Figure 3-5: Window splitting of different methods 40 Figure 4-1: One-dimensional Synthetic interferogram signal 47 Figure 4-2(a): Real part of Gabor filters in the space domain 48 Figure 4-2(b): Magnitude Spectrum of Gabor filters 48 Figure 4-3: Block diagram of the proposed method 50 Figure 4-4: Outputs of the multi-band filters of Figure 4-2 for the input signal in 4.1(b) 51 Figure 4-5: Constant bandwidth filters results 52 Figure 4-6(a): Real part of Gabor wavelet-like filters in the space domain 57 Figure 4-6(b): Magnitude Spectrum of wavelet-like Gabor filters 57 Figure 4-7: Constant bandwidth space-frequency partition 58 Figure 4-8: Wavelet-like filters space-frequency partition 59 Figure 4-9: Outputs of the multi-band filters of Figure 4-6 for the input signal in Figure 4-1 (b) 65 Figure 4-10: Wavelet-like filters results 66 Figure 4-11(a): Center frequencies of the filters for the two-dimensional case 67 Figure 4-11(b): Center frequencies of the filters for the two-dimensional case 67 Figure 4-12: Symmetric filters distribution for interferograms flattened with a pre-existing DEM 68 Figure 5-1: Example of IF estimation by the DES A-1 algorithm 77 Figure 6-1: Synthetic interferogram example 80 Figure 6-2: Outputs of three of the Gabor filters 81 Figure 6-3: Output of the frequency estimation process: 83 Figure 6-4: ID slices of estimated local frequency 84 Figure 6-5: ID slices of estimated IF when no noise is present in the signal 85 Figure 6-6: Estimation error of the proposed method for different values of a 87 Figure 6-7: Synthetic interferogram example with rough topography (a=20) 88 Figure 6-8: ID slice of estimated IF when no noise is added to the signal and a = 20 90 Figure 6-9: Performance of proposed method vs. Hypothetical Estimator 91 Figure 6-10: Phase distortion by different filters using 11 pixel window 97 Figure 6-11: Unwrapped Filtered Phase - 1.0 coherence 98 Figure 6-12: Range Frequency estimates 99 Figure 6-13: Phase noise filtering using 11 pixel windows 100 Figure 6-14: Unwrapped Filtered Phase - 0.5 coherence, 11 pixel window 101 Figure 6-15: MSE of the phase difference error for 5 and 11 pixel windows 102 Figure 7-1: Magnitude of ERS Tandem mission interferogram 107 Figure 7-2: Phase of ERS Tandem mission interferogram 108 Figure 7-3: Magnitude Frequency Spectrum of ERS Tandem mission interferogram 109 Figure 7-4: Gabor Filters Frequency centers in 2D 110 Figure 7-5: Gabor Filters in ID 111 vii Figure 7-6: IF in the range direction - Estimated on the interferogram of Figure 7.2 113 Figure 7-7: IF in the azimuth direction estimated on the interferogram of Figure 7.2 114 Figure 7-8: Filtered interferogram using the IF estimates of Figure 7.5 and 7.6 115 Figure 7-9: Filtered interferogram using the standard multilook averaging filter 116 Figure 7-10: Residues of input interferogram (Figure 7-2) 117 Figure 7-11: Residues of filtered interferogram (Figure 7-8) 118 Figure 7-12: Magnitude of Radarsat interferogram : 122 Figure 7-13: Phase of Radarsat interferogram 123 Figure 7-14: Magnitude Frequency Spectrum of Radarsat interferogram 124 Figure 7-15: IF in the range direction - Estimated on the interferogram of Figure 7.13 125 Figure 7-16: IF in the azimuth direction - Estimated on the interferogram of Figure 7.13 126 Figure 7-17: Filtered interferogram using the IF estimates of Figure 7.15 and 7.16 127 Figure 7-18: Filtered interferogram using the standard multilook averaging filter 128 Figure 7-19: Residual phase map of input interferogram (Figure 7-13) 129 Figure 7-20: Residual phase map of filtered interferogram (Figure 7-17) 130 Figure A-1: Example illustrating naturally occurring residues 139 Figure C-l: DESA-1 IF estimation for IF within [n/3,2n/3] 147 Figure C-2: DESA-1 IF estimation with modulation 148 Figure C-3: DESA-1 RMS estimation error when IF varies fast 149 Figure C-4: DESA-1 Performance for Q m o d = n/l00 150 Figure C-5: DESA-1 Performance for £2 m o d = nl\6 151 Figure C-6: Effect of the Gabor filtering for Q M O D =TZ7100 153 Figure C-l: Effect of the Gabor filtering for & M O D = TT74 154 Figure C-8: Influence of the Gabor filtering on the IF Estimation error in the absence of noise 156 Figure C-9: Effect of Gabor filtering for Cfcl , Q m o d =n/4 157 FigureC-10: Effect of Gabor filtering for a=3, Q m o d =nlA 158 Figure D-l: Input and filtered Synthetic interferogram a=50, 0.8 coherence 161 Figure D-2: Linear Phase Model 162 Figure D-3: ID slices on Figure D-2 163 Figure D-4: Noisy case 164 Figure D-5: Integration paths 167 Figure D-6: Linear vs. non-linear phase modeling - Theoretical IFs 169 Figure D-7: Linear vs. non-linear phase modeling - Estimated IFs 170 viii List of Acronyms A M - F M amplitude modulation - frequency modulation B W bandwidth D E M digital elevation model DESA-1 discrete energy separation algorithm DFT discrete Fourier transform ERS (1 or 2) Earth remote sensing satellite FFT fast Fourier transform GPS global position system LF instantaneous frequency InSAR Interferometric synthetic aperture radar JPL Jet Propulsion Laboratory M L maximum likelihood M S E mean square error MUSIC multiple signal classification algorithm pdf probability density function Q (filter) quality factor radar radio detection and ranging Radarsat Canadian radar remote sensing satellite RMS root mean square SAR synthetic aperture radar SLC single look complex-value (image) S R T M shuttle radar topographic mission STFT short time Fourier transform ID one-dimensional 2D two-dimensional i x Acknowledgements I would like to thank my supervisor Dr. Ian Cumming for his support and for a wonderful time in Europe during the conference on SAR Interferometry. Special thanks to Dr. Michael Seymour for his support throughout all the stages of the project, and to Jeremy Benson from MacDonald Dettwiler for his valuable comments on the thesis report. This work was financially supported by MacDonald Dettwiler, NSERC, BC ASI and the Science Council of BC. Disconcerted, knowing that the children were waiting for an immediate explanation, Jose Arcadio Buendia ventured to murmur: "It's the largest diamond in the world". "No," the gypsies countered. "It's ice." Jose Arcadio Buendia, without understanding, stretched out his hand toward the cake, but the giant moved it away. "Five reals more to touch it," he said. Jose Arcadio Buendia paid them and put his hand on the ice and held it there for several minutes as his heart filled with fear and jubilation at the contact with mystery. Gabriel Garcia-Marquez, One hundred years of Solitude. XI But the question you were asked, Theaetetus, was not, what are the objects of knowledge, not yet how many sorts of knowledge there are. We did not want to count them, but to find out what the thing itself - knowledge - is. Socrates, Dialogs of Plato. x i i Chapter 1 Introduction Imagine that the topography information of the whole world were available with high accuracy; we would know the elevation of mountains and valleys, the flowing courses of water, the movement of glaciers, and the expansion of volcanoes before erupting. Earth science professionals would have essential information, and the people working on SAR interferometry would be - at least - satisfied of doing something of benefit for others. Synthetic aperture radar (SAR) interferometry or InSAR, is a technique for obtaining accurate elevation maps of the Earth from radar images. When a point on the Earth's surface is illuminated by two antennas located in different positions, the elevation of this point is geometrically related to the difference of its line-of-sight distances to the antennas. The way InSAR estimates the difference of these distances can be explained using the basic concept of velocity. The distance traveled by an object (in this case a signal) equals the product of the media velocity and the observed travel time. The difference in the signal travel times is estimated from the difference in the signal phase delays. The velocity of the media, on the other hand, is assumed constant. In this way, by estimating the difference in the signals phase delays an estimate of the point's elevation is obtained. An interferogram image is computed from the signal phase delay difference of two overlapping radar images. The interferogram is processed in several stages, and after adjusting certain parameters, the scene's elevation map is obtained. Practical applications include topographic and thematic mapping. 1 In order to obtain dynamic topographic information, SAR images pairs of the scene are taken at different times. Their phase is subtracted in order to isolate the topography's displacement. Displacements of few centimeters can be measured using this differential interferometry technique. Its applications include glaciology, land subsidence, volcanoes and earthquakes monitoring. Even though many theoretical and technological advances have been achieved in the last years, SAR interferometry is not a mature technique yet. There are several difficulties in obtaining good (low levels of noise) interferograms and estimating parameters such as the baseline accurately, which makes the process of obtaining elevation maps, thus far, not fully operational. This thesis addresses one specific problem in interferometry, which is the estimation of the fringe frequency in the interferogram signal. The concept of frequency is related to the rapidity at which a variable periodically changes. In the case of interferogram signals, this variable is the phase, and a complete variation cycle is described as a fringe. The frequency of the fringes is related to the local topographic slope of the scene, and it is required in several stages of the interferogram processing. In this thesis, a new algorithm to perform fringe frequency estimation is developed and evaluated. 1.1 Thesis Contribution Motivation for Research: Local fringe frequency estimates are needed in three stages of the interferogram processing chain: coherence estimation, phase noise filtering and phase unwrapping. • In the coherence estimation stage, the frequency estimates are used to correct the effect of the slope, as described in Section 2.3. • There are three general approaches to perform phase noise filtering (they are presented in Section 2.4). If the frequency estimation approach is chosen, local fringe frequency estimates are needed for either defining the center frequency of the band-pass filter or modeling the local phase in the slope-correction filter. 2 • Finally, fringe frequency estimation has been shown to facilitate the phase unwrapping process, making the process feasible even for very noisy interferograms (Section 2.5). Thesis Scope: Given the importance of local fringe frequency in interferogram processing and the limitations in resolution and accuracy of present methods, the objective of this thesis is to investigate algorithms for local frequency estimation in interferograms. The specific objectives are to: • review current fringe frequency estimation methods, • develop an algorithm to overcome the limitations of current methods, • evaluate the proposed algorithm on synthetic and real interferograms, and • compare the proposed method against current methods. Thesis Contribution: The contribution of this thesis is the development of a new algorithm for local frequency estimation. Under the circumstances of rough topography and moderately high coherence, the proposed algorithm has advantages over other schemes because it achieves a better space-frequency resolution and therefore tracks the topographic changes of the scene more accurately. 1.2 Organization of the thesis The thesis report is divided in eight chapters. After this introduction, basic concepts on SAR interferometry are reviewed in Chapter Two, where the importance of fringe frequency estimation in interferogram signals is further explained. In Chapter Three, existing approaches to perform fringe frequency estimation are reviewed. Several methods have been proposed in the literature. They vary in two aspects: the way they split the image in order to obtain the frequency estimates, and the way the frequency estimation itself is performed. The first aspect influences the spatial resolution of the estimation, while the 3 second aspect influences the estimation accuracy. Both aspects are analyzed for the existing methods, in order to design a method that overcomes the limitations of the existing ones. In Chapter Four the proposed method is described. The proposed method is based on the A M - F M (amplitude modulation - frequency Modulation) signal representation model, which is used in other areas of signal processing theory, like in image analysis and compression. The mathematical framework of the signal model is presented, and the application of the proposed method to one and two-dimensional signals is explained. One important part of the proposed method is the estimation of the instantaneous frequency. The selected algorithm for this stage is explained in Chapter Five. Chapter Six is devoted to the evaluation of the proposed method. Synthetic interferogram signals with different spatial distribution and levels of noise are used to evaluate the performance of the proposed method. One of the applications of local fringe frequency estimation is phase noise filtering in interferogram signals. A new method based on the estimated frequencies is presented and its performance is compared against two existing methods. Chapter Seven, on the other hand, is devoted to the application of the proposed frequency estimation and phase filtering methods to real interferogram signals. Finally, conclusions and suggestions for future work are made in Chapter Eight. In addition to the main text, four appendices explain specific topics, which are referred to the thesis. 4 Chapter 2 SAR Interferometry Background The subject of this thesis is the development of a new algorithm to accurately estimate the local fringe frequency in SAR interferograms. In this chapter, a review of SAR interferometry is made and the importance of local fringe frequency estimation is further explained. There are several stages to produce and process an interferogram in order to obtain topographic information. These stages form a processing chain, which is reviewed in Section 2.2. Among all the processing stages, there are three where fringe frequency estimation is relevant. These three stages are further explained in Sections 2.3 to 2.5. 2.1 SAR Interferometry Principles Synthetic aperture radar (SAR) images are obtained by processing the data collected by an aircraft or satellite with a side-looking antenna. The system transmits a stream of radar pulses and records the back-scattered signal from each pulse. A digital processor is used to focus the received radar data into a high-resolution image [37]. Two types of information are included in the processed SLC (Single Look Complex) image. The first type corresponds to the magnitude component | / , ( x ) | , which represents the backscatter of microwave energy of the ground, and depends on the properties of the surface such as slope, roughness and dielectric constant. The second type of information is its phase component 0i (x), which is related to the distance of closest approach of the antenna and the ground patch. In the absence of noise, a SAR image is represented by: 5 /,(*) = \I1(x)\txpU0l(x)) (1) where x = (xx,x2) are the image pixel coordinates, which are commonly referred as the range and azimuth directions (see Figure 2-1). An interferogram is obtained by multiplying the SAR image Ix(x) by the complex conjugate of another overlapping co-registered SAR image I2(x) obtained using a slightly different antenna position: I(x) = WxVtW = | A || /21 expCyC^Cx)-692(x))) (2) The relationship between the resulting interferogram phase and the scene's topography is explained as follows. If the same patch area A of the Earth's surface is observed from two different points SA, SB as illustrated in Figure 2-1, the phase shift between the images acquired by the antenna position SA and the antenna position SB is proportional to the difference in line-of-sight distances: Ajr <S> = —(RA-RB) (3) where <I> = 0] -62 is the interferogram phase of I(x) in (2), and RA,RB are the slant range distances in Figure 2-1. The phase shift in (3) can be approximated by: [1] . Anfbx zbHx\ A d d1 where: X is the radar wavelength, 6 b is the normal baseline between the two antennas, x is the ground range distance, H is the antenna altitude, and d = -JH2+x2 . The right part of (4) has two contributions to the phase variations. The first term b xl d does not depend on the elevation z and corresponds to the fringe pattern in the case of flat terrain, it is referred as the orbital fringes phase. The second term zbH xl d3 is proportional to z and contains the topographic information, it is referred as the topographic fringes phase. The estimated phase of the interferogram signal in (2) can be used to infer the scene's topography. This result contains the essence of the SAR interferometry technique or InSAR. Figure 2-2 illustrates the magnitude of the interferogram obtained from two SLC images acquired by the satellites ERS1 and ERS2 (tandem mission) in a region of the Chilcotin River's upper valley in B.C. (Canada). Figure 2-3 illustrates the interferogram phase. The almost parallel fringes in Figure 2-3 correspond to the orbital fringes in (4). The topographic fringes, on the other hand, are immersed in the regular pattern and cause local variations in the phase. The data acquisition parameters are listed in Table 2.1. Another type of information that can be extracted from two overlapping SAR images is its coherence magnitude. The coherence magnitude is a sensitive indicator of the similarity of the two SAR images. The coherence magnitude is used as a quality measure of the interferogram. Therefore, it can be used to estimate the accuracy of the obtained elevation map. It is also used in image landcover classification. The maximum likelihood estimate of the coherence assuming constant interferogram phase is given by: [3] N X VAm.n e X P ( ~ 9l.m. J ) (5) 7 Normal basel ine H l \ range C'Ve, C f ; o n Figure 2-1: Interferogram acquisition geometry Parameter Value Units A l t i t u d e 7 8 2 K m R a d a r w a v e l e n g t h 0 . 0 5 6 6 M e t e r s I n t e r f e r o g r a m B a s e l i n e 4 2 M e t e r s A c q u i s i t i o n d a t e S L C 1 S e p t . 1 0 / 9 5 -A c q u i s i t i o n d a t e S L C 2 S e p t . 1 1 / 9 5 -S c e n e ' s c e n t e r l a t i t u d e 5 1 ° 3 5 ' 2 0 " N -S c e n e ' s c e n t e r l o n g i t u d e 1 2 1 ° 5 9 ' 2 0 " W -S p a t i a l R e s o l u t i o n 1 2 . 5 x 1 2 . 5 M e t e r s Table 2.1 Example of interferogram acquisition parameters. 8 Input Interferogram - Magnitude <— Range direction Figure 2-2: Magnitude of ERS Tandem mission interferogram The interferogram corresponds to a 5 Km x 5 Km area in the upper valley of the Fraser River in B.C. (Canada). 9 Input Interferogram Figure 2-3: Phase of ERS Tandem mission interferogram The topographic phase information is immersed in the regular pattern and cause local variations in the phase. Note: The wrapped phase [-n,n] is normalized to [0,1] as shown by the gray scale 1 0 where: V is the coherence magnitude, y/ is the maximum likelihood (ML) estimate of the interferogram phase, /,•_,„_„ exp( j(Qinhn )) is the complex value of the pixel (m,n) in the averaging window for SLC image i, i = 1,2. In the next section, the processing steps to produce and process an interferogram are reviewed. 2.2 Interferogram Processing Chain When the two SLC images are acquired in different passes (repeat-pass interferometry), the Doppler centroid of the two images is not necessarily the same. This implies that the images sample different parts of the azimuth spectra and the non-overlapping or un-correlated parts generate noise in the interferogram. Azimuth pre-filtering to the common azimuth frequency bandwidth has been reported to increase the coherence magnitude [2] and constitutes the first step to obtain an interferogram in the processing chain of Figure 2-4. Range filtering is usually performed to extract the common range spectra as well. The second stage in the processing chain is image registration. Typically, one of the images is shifted and rotated to coincide with the other one. Accuracy in the registration is very important, so over-sampling is first performed to obtain registration errors less than a 1/10 of a pixel and to prevent aliasing when the image multiplication is performed. There are several methods for registering the two SAR images [3]. The effectiveness of each method depends upon the input data characteristics. After registration, the coherence magnitude is computed as in (5) based on the pixel values of both SLC images. To form the interferogram, both images are multiplied pixel by pixel as in (2). The resulting interferogram exhibits two contributions in the phase variations as given by (4): the flat-earth fringes which are produced by the acquisition geometry, they are the same as if the scene's topography were flat, and the topographic fringes which are phase variations produced by the real topography. 11 The average magnitude spectrum in the range direction of the interferogram in Figures 2-3 and 2-4 is illustrated in Figure 2-5. The maximum range frequency peak is located close to the frequency origin; this frequency corresponds to the frequency of the flat-earth fringes. According to the spectral shift theorem [4], layover areas in an interferogram correspond to negative frequencies. Therefore, they can be filtered out using a half-band filter. In Figure 2-5 the dashed line draws the ideal frequency response of the half-band filter. By multiplying the interferogram with the complex conjugate of a phase function that models the flat terrain, the flat-earth fringes are removed. Flattening is normally performed before any further processing because it moves the interferogram phase derivatives to lower values, which helps to simplify subsequent processing steps. After flattening, the signal is sub-sampled to its original resolution and phase filtering is performed to reduce the phase noise. The last stage in the processing chain is phase unwrapping. In the phase unwrapping stage, the phase of the processed interferogram is unwrapped to produce terrain elevations. Unwrapping means adding or subtracting integer multiples of 2n to produce unambiguous phase values whenever possible. Figure 2-6 illustrates the extracted wrapped phase from the transect line A in Figure 2-3. Notice that the phase values are bounded within the interval [—TC,Tt). Fringe frequency estimation is relevant to three stages in the processing chain: coherence estimation, phase noise filtering and phase unwrapping. These stages are shaded gray in the block diagram of Figure 2-4. The use of the frequency estimates in these stages is explained in the following sections. For the phase noise filtering and phase unwrapping stages, their different processing approaches are reviewed. 2.3 The Use of Frequency Estimates in Coherence Estimation The M L coherence estimator in (5) assumes that the phase is constant. Therefore, for rough topographic areas, the estimated coherence is distorted. By locally modeling the phase and correcting the effect of the slope, a better coherence estimation is achieved [4,5,6]. In [6] the phase behavior is assumed linear in a window around each pixel P in the image. Based on the local frequency estimates CQx,CQy, the local topography is modeled as: 12 I,(x) SAR IMAGES I 2(x) AZIMUTH PRE-FILTERING AZIMUTH PRE-FILTERING IMAGE REGISTRATION HALF-BAND FILTER PHASE NOISE FILTERING , i , PHASE UNWRAPPING Unwrapped phase Figure 2-4: Interferogram Processing Chain The stages where Fringe Frequency Estimation is relevant are shown shaded. 13 - 2 - 1 0 1 2 3 Frequency in Range direction —> rads/sample Figure 2-5: Interferogram average magnitude spectrum in the range direction The ideal frequency response of the half-band filter is drawn by the thick dashed line. Wrapped interferogram phase 20 40 60 80 100 120 140 160 180 200 Figure 2-6: ID slice of the phase in the transect line A of Figure 2-3 Notice that the phase values are bounded within the interval [-n,n). The horizontal axis is sample number in the interval shown by transect line A in Figure 2-3. 14 (j)top{m,n) = coxm + coyn (6) where m, n are the coordinates of the pixels in the window. Then, the terms 7 / > m > n exp(jdi ) with i = 1,2, in (5) are multiplied by exp(-j(/)top(m,n)) to correct the effect of the local slope. In this way, the interferogram is locally flattened before the coherence is estimated. 2.4 Phase Noise Filtering Approaches There are several sources of decorrelation noise in interferogram signals, namely: system noise, quality of focusing and registration, temporal decorrelation, and spectral or baseline decorrelation [7]. It has been shown that filtering the interferogram prior phase unwrapping stage greatly reduces the complexity of the process and leads to better results [6,8]. Different approaches to interferogram phase noise filtering have been proposed. They can be analyzed using the theoretical framework of signal estimation theory [9] and can be grouped into three classes: • constant phase estimation approach, • frequency estimation approach, and • post-filtering techniques. 2.4.1 Constant Phase Estimation Approach In this approach, the phase is assumed constant in a small region of the image. The noise and phase statistics are modeled and optimum filtering is performed seeking to optimize a criteria. In the presence of additive noise, the observed SLC image magnitude | 7, | is represented by a deterministic part a and random part £ , | 7, (x) \ = a£. The random part £ obeys a Raleigh distribution derived from the vector sum of two orthogonal normal distributions. On the other hand, the phase of 7, is the sum of a deterministic part 0X and a random part T]l, 15 argf/^x)} = 6l+rjl. The random part 77, can be assumed uniformly distributed in the interval [-Tt,?:). The observed signal is represented by: I(x) = /,(je)/2(jc)* = | / 1 | | / 2 | e x p [ 7 ( ( ^ - ^ ) + ( ^ , - ^ ) ) ] (8a) Making (j) = 9x-02 and T] - T]x-t]2, the deterministic and random parts of the interferogram phase respectively, I(x) = | / | e x p ( ; » = \I\exp[j(cp+r])] (8b) The probability density function (pdf) f^ity) is derived in [10] for single-look and multi-look phase distributions. The pdf is a function of the coherence of the interferogram among other parameters. If the interferometric phase is assumed constant along m samples of the image, the following model is used to estimate the phase: y/i = tp + r]i for i = l,...,m (9) where: tp is the deterministic phase which is assumed constant in the analysis window, and T]i is the random part of the phase as described by (8b). Assuming that all pixels in a square window have the same deterministic phase value, a maximum likelihood (ML) phase estimate is obtained by coherently averaging the complex values in the window [11]: cp = a r g { £ I(m,n)ejn,n'n)} (10) m,n 16 where I(m,n) e • 7 ¥ ( m ' " ) are the input interferogram samples in the window, and ip is the estimated phase. The phase estimator in (10) is named standard multi-look, and it is found to reduce the noise variance at the expense of losing spatial resolution. A better estimation is achieved by using only samples belonging to a directional window. The local statistics estimator proposed in [12] seeks to minimize the mean square error (MSE) between the estimated and the real phase signal. In this case, the phase estimation problem consists in estimating 0 from the samples y/i in (9), such as the mean square error (<p-(p)2 is minimized. The estimator is given by: where: y/i j is the locally unwrapped input phase, y/i j is the local mean unwrapped phase, and btj is given by: EK^j-w,,,)2] The filtering is performed by the following steps: i) The coherence map is computed. ii) Pixels in the moving window are unwrapped locally around the average value of the center 3x3 pixels. iii) Sixteen directional windows (in TC/16 steps) are used to determine the direction with minimum variance o~n2. It is assumed that in this direction the phase is constant. 17 iv) The phase is estimated using (11) where i,j are the coordinates of the pixels in the window that has lowest variance. To obtain the variance an2 the coherence computed in step i) is used in an analytical function that relates variance and coherence [10]. v) When applying equation (12), b can become negative. If so, it is set to zero to ensure the weight is always between 0 and 1. In this algorithm, low coherence areas are smoothed more heavily than high coherence areas, which is a desirable characteristic of the filter. At the same time, by filtering in the orientation of minimum variance, the resolution is not degraded if the neighborhood pixels along that direction have the same phase (same elevation). An alternative method to computing (11) without performing local phase unwrapping is reported in [13]. The authors reported a qualitative evaluation of the filter showing that the filter substantially reduces noise, facilitating the phase unwrapping process. In addition, the results were similar to that obtained by filtering the data with directional Gauss filters [14], but the preservation of the phase gradient and orientation of the fringes was not as good, especially in low-coherence areas. A refinement of the algorithm described in the previous section is proposed in [17]. In areas where the phase gradient is very high, narrow signal fringes are lost after spatial filtering. In those areas, a parametric estimation algorithm is used to remove the high frequency signal phase before filtering is carried out. The application of the modified filter to synthetic data shows that the combination of the parametric estimation with the refined local statistics filtering retained both high-density fringes and spatial resolution. 2.4.2 Frequency Estimation Approach In the frequency estimation approach, the phase is not assumed constant anymore. Instead, it is assumed to be linearly dependent on the spatial coordinates. This converts the estimation problem into a spectral estimation one, where the best complex sinusoid signal is estimated from noisy samples of the signal. In this approach the measurement model is a complex sinusoid signal of unknown amplitude A and frequency fa . In one dimension, the signal phase model is: 18 Vi = <Pi+rli for i = 1,..., m (13) where cpi varies linearly with i, i.e. = 2nf()i. The estimation problem is to estimate fa based on the m noisy observations xi: xt = Aexp(j( ft +r]i)) = A exp( j (In fai + 77,.)) (14) Several frequency estimation methods have been proposed for interferogram images [15,16,18]. In Chapter 3, some of these methods are reviewed. After the frequency has been estimated, a band-pass filter centered on the estimated frequency fa is applied to the data. As an alternative to band-pass filtering, slope-correction filtering has been proposed [6]. In this case, for every square window in the interferogram, the estimated frequency (fx,fy) is used to locally model the phase as in (6). By multiplying the input image with the 2D complex function exp(- j2n{fxm + f n)) the effect of the slope is corrected. After this local slope correction, the phase is assumed to be a stationary process that only varies because of the noise and the differences between the exact fringe pattern and the first order model. Then, standard complex averaging is carried out, reducing the noise variance. A detailed description of multi-look filtering is presented in Section D.2 of Appendix D. 2.4.3 Post-filtering Approach Non-linear image processing techniques have also been applied to phase noise filtering. In [19], the authors use three techniques: a modified median filter, a statistical mode filter, and a morphological filtering. Results on real data showed a reduction in the number of residues1. From the three approaches reviewed, only the frequency estimation approach - Section 2.4.2 - makes use of local fringe frequency estimates. 1 The concept of residue is explained in Appendix A, together with an algorithm to detect them. 19 2.5 Main Approaches to Phase Unwrapping In the absence of residues, a simple integration procedure in which phase gradient estimates (partial derivatives) are integrated along vertical or horizontal lines suffices to produce the unwrapped terrain elevation. In practice, however, there are noise and aliasing effects that cause phase inconsistencies in the interferogram. Residues reveal the existence of these phase inconsistencies. The presence of many residues in the interferogram indicates that the phase is too noisy or that the fringe pattern is not regular enough to allow straightforward phase unwrapping. Unwrapping difficulties result from two kinds of corruption: • Decorrelation noise that affects phase estimation. This noise is high when the interferogram coherence is low. In this case, fringe patterns exist but they are corrupted by noise. • Areas in which there are no fringes or the fringe continuity has disappeared. This happens in the presence of shadows and layover, or in the case of abrupt topography (such as cliffs) which the interferogram pair is unable to map. In both cases, the integration of simple phase gradients propagates and accumulates phase inconsistencies giving rise to large errors in the final result. The proposed approaches to overcome this situation and effectively perform 2D phase unwrapping can be separated in two groups: local and global approaches. Both benefit from prior fringe frequency estimates and are briefly described below. • Local approach: In this approach, discrete phase differences are integrated along selected paths. In order to avoid error propagation, connection of residues is performed to define barriers, which the integration paths are not allowed to cross [20]. When the noise creates many residues they can not be correctly connected. Therefore, the use of the local fringe frequencies to reduce the number of residues by adaptively filtering the local fringe pattern, help phase unwrapping when using this approach. • Global approach: This approach is based on the minimization of the difference between the phase gradient measured on the interferogram and the gradient of the unwrapped phase [21]. The use of discrete phase differences corrupted by noise leads to biased estimates of the gradient [8]. Thus, the obtained least-squares solution is 20 undervalued, resulting in an underestimated topographic slope. Local frequency estimates can replace the discrete phase gradients, leading to more accurate phase gradient estimates and therefore more accurate elevation maps. 21 Chapter 3 Local Fringe Frequency estimation - Existing methods In order to estimate the local fringe frequency of an interferogram signal, researches have used different approaches. In this chapter those approaches are reviewed in a structured way. The most straightforward way to find about the frequency content of a signal is by computing its Fourier transform. However, this transform assumes that the signal is stationary, which is not the case when the interest is to find the local fringe frequency of an interferogram signal. The short time Fourier transform (STFT) for non-stationary signals is considered in Section 3.1, where basic definitions are given with an analysis of its space-frequency resolution. The STFT is used to compute the periodogram, which is an estimate of the local power spectrum and serves as a tool to estimate the local frequency of a signal. In Section 3.1.3 results using this approach are reviewed with an implemented example. In an attempt to improve the accuracy of the local fringe frequency estimation, a parametric method is reviewed in Section 3.2, and a method that estimates the local fringe frequency at different resolution levels is analyzed in Section 3.3. Finally, the adequacy of each approach is discussed based on the assumptions made in each case. 22 3.1 Short Time Fourier Transform Because interferogram phase signals are non-stationary the use of the short time Fourier transform (STFT) 2 is considered instead of the traditional Fourier Transform. The STFT of a continuous function s(x) is defined as: oo STFT\x, co) = \s(x)g(x-z)e-j0*dx S,gzL2(R) (1) —oo This function can be described as the Fourier transform of the signal s(x) previously windowed (multiplied) by the function g(x) around the position coordinate X- The variable CO is continuous in the interval (-00,00) and corresponds to the angular frequency of the complex functions e~jm. As the window function is shifted in space over the whole signal and consecutive overlapped transforms are performed, a description of the evolution of the signal spectrum is achieved. This method assumes that the signal is stationary over the limited window g(x), which is a valid assumption if the window is short enough. This window has an important effect upon the frequency resolution of the space-frequency distribution of the signal. If a short window in space is used to compute (1) the transform will exhibit good resolution in space, but the frequency resolution -defined as the ability to discriminate complex sinusoids of different frequencies - will be low. Alternatively, if a longer window in space is used in order to increase the frequency resolution, the space resolution (localization) will be poor. In the case of interferogram phase signals, it is important to obtain a space-frequency representation as accurately as possible in both space and frequency in order to obtain accurate frequency estimates of the signal. In order to quantify the effect of the window function g(x) on the space and frequency resolution of the STFT, let the square root of the energy distribution be a measure of the effective spatial width of the signal. 2 The name short time Fourier transform is used even though space is the domain of interferogram signals. 23 \x2g(x)* g(x)dx (Ax)2 = ^ (2) jg(x)*g(x)dx . —oo where * denotes complex conjugate, and the effective bandwidth of its Fourier transform G(co) be also defined in terms of its energy distribution. oo jco2G(co)*G(co)dco (Aco)2 = (3) JG(co)*G(co)dco —oo It has been reported [22] that the product of (2) and (3) has a lower bound given by: AXAOJ > (47Z-)-1 (4) Equation (4) is called the uncertainty principle in analogy with the Heissenberg uncertainty principle in quantum mechanics. This equation implies that space and frequency resolution can not be made arbitrarily small. One resolution must be traded for the other. For the STFT, the resolution in space and frequency are constant values independent of (%,QJ). These resolutions only depend on the type and length of the window function g(x). This constant joint space-frequency resolution results in the STFT covering the space-frequency plane with a uniform array of cells. 3.1.1 The STFT as a Signal Decomposition Transform The intention of this section is to present the STFT in a broader range of signal decomposition transforms, where basis functions other than complex sinusoids could be used to more efficiently decompose the signal for frequency analysis purposes. 24 Equation (1) was analyzed as the Fourier transform of the windowed signal s(x)g(x-X), but for real-valued signals it is also valid to describe this function as the decomposition of the signal s(x) using the basis functions g(x — X)e~JCM • These basis functions form a complete set of functions that when combined in a weighted sum can be used to construct a given signal. In the case of the STFT these basis functions are complex sinusoids e~]m windowed by the function g{x) centered around X • Using this description it is possible to write a general equation for the STFT in terms of inner products of the basis functions K m(x) = g(x- %)e~jm and the signal s{x) . By mapping the signal onto these basis functions, a space-frequency description of the signal is generated. As an example, consider the one-dimension synthetic phase interferogram signal of Figure 3-l(b). The instantaneous frequency of this signal covers the interval [0,7i/4), as illustrated in Figure 3-l(d). In Figure 3-l(c) its discrete Fourier transform (DFT) has been computed. In this case, the signals are discrete with N samples so (5) becomes: where m,n are sample numbers of the discrete sequences s, g , and £2 denotes discrete angular frequency in rads/sample. Notice that STFT(-) has been changed to S$F (m, Q) for the discrete case and it is simply called a space-frequency representation. The result of computing (6) by using rectangular function windows of length N=32 is illustrated in Figure 3-2. In this case, the basis functions are discrete sequences of complex exponentials with frequencies: Q = 2nx{ 0 , 1 / N , l — l/N}. (5) SSF(m,Q) = ^s(n)g(n-m)e -jQn (6) 25 a) Phase signal 50 100 150 sample number • 200 250 c) Interferogram spectrum 2 4 frequency (rads/samp) -> A 0. | o . c § 0 . Cr b) Interferogram signal TTT 50 100 150 200 250 sample number -> d) Instantaneous frequency / / / • / / / / . . . / / 50 100 150 200 sample number -> 250 Figure 3-1: One-dimensional Synthetic interferogram signal a) Quadratic phase function of the signal Q(x)= Kx2 b) Real part of the Interferogram signal s(x)=eii2<x> c) Magnitude of signal spectrum. d) Instantaneous frequency in radians computed as the derivative of £2(x). 26 Space Frequency representation Figure 3-2: Space -Frequency representation of signal in Figure 3-l(b) This figure shows the evolution of the spectrum. A s the input sample number increases, the center frequency of the local spectrum increases. 27 Although the signal is locally monochromatic this space-frequency representation shows a locally band limited signal whose central frequency is increasing. Other windows can be considered in order to improve the frequency resolution, but they all have the restriction of the uncertainty principle (4). Alternatively, the window size can be adapted to the local frequency content of the signal (see Appendix B). It can be easily shown (like in [23]) that Gaussian functions used as windows achieve the lower bound of (4). These functions and their Fourier transform have the form: g(x) = e-*2'2*2 Gift)) = 0 -V2^~" a 2 ( y 2 / 2 (7) These functions are named Gabor filters [36] and they are used in the next chapter to effectively filter the interferogram signal. 3.1.2 Use of the Periodogram in Local Frequency Estimation The periodogram of a continuous signal x(t) is by definition the stochastic process: ST(co) = -^\XT(vf (8) where XT (co) is the Fourier Transform of the truncation of the signal x(t): T XT(co) = jx(t)e~ja"dt (9) -T For large values of T the periodogram ST(co)h an unbiased estimator of the power spectrum S(co) of the signal x(t). However, if S(co) is not nearly constant in an interval of the order of l/T, the periodogram is a biased estimate of S(co) [24]. 28 In order to reduce the bias the signal x(t) can be windowed by a function c{t) yielding the modified periodogram: SC(GJ) = — C IT jc(t)x(t)t - jcot dt\ (10) Alternatively, in order to reduce the variance and the bias of the estimation, performing convolution with a spectral window function W(oj) can smooth the periodogram. Sw(co) = —ST(co)*W(co) (11) In [24], windows W(co) are found to minimize either the variance or the bias of the estimation, which are two conflicting requirements. Once a good estimator is found, the frequency for which there is a peak in the power spectrum S(co) constitutes a maximum likelihood (ML) estimator of the signal local frequency. In other words, among all possible local frequencies cot the peak of the power spectrum S(co) of the signal x(t) indicates which frequency a\ is the most probable to occur given the observed signal samples x(t), which are used to compute the estimate of S(co). Figure 3-3 illustrates the result of the process when applied to the synthetic signal of Figure 3-l(b). Figure 3-3(a) illustrates the input signal and one example of an analysis window of 32 samples. The window is shifted over the whole signal to compute the local periodogram and to estimate the frequency by finding the maximum of the local periodogram. A circle in Figure 3-3(b) shows each frequency estimate, and the theoretical frequency - defined as the derivative of the unwrapped phase of the signal - is 29 illustrated in continuous line. Interpolation by eight samples has been applied to the computed periodogram in order to find the maximum more precisely. The accuracy of the estimation of the frequency C0k depends on the accuracy of the estimation of S(C0) . The error presented in such estimation is due to the fact that there is not only one frequency present in the analysis window, but a continuous range of increasing instantaneous frequencies of same amplitude. In order to increase the localization property of the estimation, the size of the window must be made small, but this decreases the resolution of the frequency representation since the frequency bins become wider. The same trade-off between localization and frequency resolution explained in Section 3.1 is involved in the choice of the windowing function c(t) in (10). Long analysis windows will behave better, but will have poor localization properties. On the other hand, notice that by producing a smoother version of the periodogram as in (11), its frequency resolution is reduced. 3.1.3 Reported Works Using this Approach In [25] the authors used overlapping windows to locally estimate the spectrum of the signal and adaptively filter the interferogram. The procedure is as follows: i) The interferogram image fSijJ is divided into windows of 32x32 pixels, with 50% overlap in the x and y direction. ii) The 2D FFT is computed on each window. iii) The following operator is applied to each pixel in the frequency domain: z,j = I -v,7 | a *,j (12) where a belongs to the interval [0,1]. 30 iv) On each window its inverse FFT is computed, and for the overlapping pixels a linear taper is used. In the work reported above no M L frequency estimation is performed. The signal is adaptively filtered according to the locally computed periodogram. The filtering depends on the value of a. As a is increased from zero, the difference in amplitude between low and high magnitude values is increased. For a =1 strong filtering is applied, the signal band-pass area is halved and the SNR squared. This filter increases the magnitude for the frequencies where there is significant signal content, and it reduces the magnitude for the frequencies where there is no significant signal (ideally only noise). However, this assumption is not valid for areas where there is rough topography or there is low coherence. Figure 3-4 illustrates this point, the filter in (12) with a =0.7 has been applied to a 32x32-pixel window of an interferogram image. In Figures 3-4(a) and (b) the input and output interferograms are shown, and in Figures 3-4(c) and (d) their corresponding 2D-DFT. Due to the low frequency resolution of the DFT, the dominant frequency in the window can not be easily estimated. It can be concluded that this filtering approach performs a simple and effective filtering. Its disadvantage is that the window-overlapping and linear taping schema does not follow fast changes in the signal spectra. The analysis window can be chosen wide enough to increase the frequency resolution but the small signal length required to maintain the stationary condition in interferogram phase signals, limits the resolution of this approach to estimate the local fringe frequency. The interferogram signal can be over-sampled before computing its spectra, in order to obtain estimations with higher resolutions in frequency, but this highly increases the computational load. 3.2 Two Dimensional MUSIC Algorithm Parametric spectral analysis has been broadly used in different applications to improve the frequency resolution of frequency estimation. The authors in [16] proposed 31 an approach that uses an efficient version of the two dimensional (2D) multiple signal classification (MUSIC) algorithm. The authors subdivided the image in square cells and assumed that the signal is locally monochromatic in each cell. The MUSIC algorithm estimates the frequency vector (cox,COy) for each cell. The MUSIC spectral estimator models the signal as a sum of K sine waves of different frequencies (CO*,0)J) i £ [1,...A'], with random amplitudes Ai and variances 2 cr corrupted with additive white noise. K *(n) = ^Akcos(cvkn + y/k) + v(n) n = 1,2,...,N (13) k=i The additive white noise is assumed non-correlated with the signal, and with variance Gv . The data length N of the received signal is assumed to be large compared to K. In one dimension, the frequency estimation problem is stated as follows: Given the samples x(l),x(2),...,x(N) of the received signal (x(n)J, estimate the frequencies CO], tQi,..., o\ of the individual complex sinusoids contained in the received signal. To solve this problem the correlation matrix whose elements are the values of the auto-correlation function rx(m) is approximately computed using N samples x(n) of the input signal. rxQ) rx<X) rx(0) rx(M-l) rx(M-2) rx(M-l) rx(M-2) rx(0) (14) 32 a) Interferogram phase signal and analysis window -1 \ I / ' \ / \ / 1 \ / 1. rt \ \ \ 1 n V 1 \ l l \ 1 \ / 1 \ (1 / 1 \ I I _ | _ L 1 _ _ l ; 1 \ 1 \ / w \l V V V v ' ' " I " 1 1 1 1 1 100 150 200 250 sample number -> b) Theoretical and estimated Frequencies 50 100 150 200 250 sample number -> Figure 3-3: M L Frequency estimator based on the local periodogram The periodogram is computed for the samples in the analysis window (the rectangular box shows an example). b) The M L frequency estimator finds the frequency at which the local periodogram has a maximum. 33 INPUT INTERFEROGRAM 32x32 window 5 10 15 20 25 30 INPUT 2D FFT 0 0 OUTPUT FILTERED INTERF. 5 10 15 20 25 30 OUTPUT 2D FFT x 107 2 v 0 0 Figure 3-4: Adaptive filtering using 32 pixels F F T a) Real part of 32 x 32 pixel window of a Radarsat interferogram b) Real part of filtered interferogram. - Notice that there is not only dominant direction in the window. c) 2D DFT magnitude of input window in a) d) 2D DFT magnitude of output window in b) Notes: For cells a) and b) the x and y axes correspond to range and azimuth directions. For cells c) and d) the x and y axes correspond to spatial frequency in the range and azimuth directions. 34 The dimension of this matrix is MxM where M is the number of lags used in the auto-correlation function. The value M is assumed larger than 2K i.e. twice the number of frequency components present in the signal. The eigen-vectors of the correlation matrix in (14) are shown to partition the RM space into two subspaces: a signal subspace that is spanned by the eigen-vectors associated with the 2K largest eigen-values of Rx, and a noise subspace that is spanned by the eigen-vectors associated with the smaller eigen-values of Rx, the multiplicity of which is M-2K [26]. In this way, the most representative frequency o\ of the sequence x(n) is estimated by finding the largest eigen-value of Rx. For interferogram signals in two dimensions, the frequency estimation problem is stated as: Given the samples of each square cell fx(m,n)J, estimate the angular frequency (0JX ,C0y) of the most representative complex sinusoid contained in each cell. The most representative frequency in the cell corresponds to the largest eigen-value of the matrix Rx. The computational load to estimate the eigen-vectors in two dimensions is very high, but the authors in [16] implemented an efficient way taking advantage of the block structure of the auto-correlation matrix. The authors reported that for small window sizes the noise reduces the dynamic range of the eigen-values, making the separation of the signal space from the noise space more difficult. This method has the advantage of high resolution in frequency, but the assumption that the interferogram complex signal for each cell constitutes a 2D one-frequency complex sinusoid is equivalent to say that all the pixels in the ground elevation fit an inclined plane for each cell. This assumption is not valid if the topographic slope of the scene varies within the cell. For scenes with rough topography it is very likely that the windowed portion of the image does not contain one single frequency for each cell as illustrated in the following example. 35 Figure 3-5(a) shows a raw interferogram of the Chilcotin area in BC (Canada) processed from two SLC images acquired by the satellites ERS1 and ERS2 , before flattening. For the 16x16-pixel cell named B in Figure 3-5(c) there will be at least two predominant frequencies, as illustrated by the background image that corresponds to the instantaneous frequency of the signal. After having selected one frequency as the most dominant, the 2D MUSIC algorithm does not take into account the others since there is no information about their location within the cell. If the image is divided in smaller windows, the problem of not having enough samples to estimate the frequency is faced. The authors in [16] proposed a scheme to adapt the window size based on the coherence map for the filtering stage, but not for the frequency estimation stage on which the filtering is based. 3.3 Multi-grid method In [27] a method was proposed to estimate the local fringe frequency based on different resolution cells. The advantage of using this approach is that the estimations over large windows can be used to reduce the effect of the phase slope on the estimation in smaller windows. The process is illustrated in Figure 3-5(d) when the reduction factor equals 2. The largest window w? - corresponding to the lowest resolution - is used to obtain an estimate of an average frequency / j over the whole window w?. The frequency / j is then used to flatten the phase over the next smaller window W2 - the next higher resolution level - by a complex multiplication by e ; 2 ^ 3 . The estimated frequency from window W2 is the difference frequency d2 = f2- f3, where fi is the frequency that would had been estimated from window W2 without the complex multiplication. The frequency/2 is calculated and then used to reduce the phase variation in the next higher resolution window wj on the original signal, giving the difference estimate: 9, = fx — f2 . 3 : T h e a c q u i s i t i o n p a r a m e t e r s a r e l i s t e d i n T a b l e 2 . 1 ( C h a p t e r 2 ) 36 In general, given N+l resolution levels with the 0 resolution being the highest resolution, the estimation result is a series of difference frequency estimates, dn, where n = 0,..., N and d N = fN. In the absence of aliasing, the sum: n=N fo = X a„ ( i s ) n=0 is the frequency estimate at the 0 t h level, corresponding to the smallest desired estimation window size. For any resolution level m, the estimated frequency is: n=N f,n = X 5, (16) Since each estimate is made up of a sum of difference frequencies - which are relatively small - the effect of slope on the estimation is approximately removed. To estimate the local frequency at each resolution level m the authors proposed to use the centroid of the local spectrum [33] after multiplying by eil7fin (local slope correction). The centroid of the local spectrum is given by the spatial frequencies in the x and y directions (fx,fy), where: . r = arg{ z(i,k)z(i + l,k) } fy = arg{ z\i,k)z(i,k + l) } (17) 37 For a window of one pixel size, the local frequencies fx= arg{z (i,k)z(i + l,k)} and fy = arg{z*(i,k)z(i,k + 1)} are simply the calculation of the wrapped phase differences in x and y. The advantage of the multi-grid method is that by reducing the effect of the phase slope iteratively, the phase variation is reduced and then it is possible to locally filter the interferogram by standard pixel averaging. The condition for the method to reduce the phase variation of the higher resolution levels is that the phase is monotonic in the analysis window. This condition is not met for scenes with rough topography. However, even in those cases the method reduces the overall phase variability. 3.4 Summary and Conclusions Three methods that perform local fringe frequency estimation on interferogram signals were described in this chapter: I. M L estimator based on local periodogram (Section 3.1.2) II. 2D MUSIC algorithm (Section 3.2) III. Multi-grid method (Section 3.3) These methods are different in two aspects: first, the way they split the image in order to obtain the set of samples {SyJ that are used to estimate the local fringe frequency; second, the way such local fringe frequency estimation is performed. In order to illustrate the first aspect, Figure 3-5 shows the window splitting performed by each method. The background in figures 3-5(b)-(d) corresponds to the magnitude of the locally computed instantaneous frequency for the interferogram signal in figure 3-5(a) using the method explained in next chapter. In the first method the signal is divided in windows of 32x32 pixels which are shifted 16 pixels in each direction, this is illustrated by the windows marked as A and A ' in Figure 3-5(b). 38 In the second method the window size can become smaller, but there is a limit when the noise brings together the eigen-values and it is not possible to distinguish the signal from the noise. In Figure 3-5(c) a 16x16 pixels window (marked as B) is illustrated as an example. Notice that even for this small window size the signal can not be assumed monochromatic for each cell in the image. The authors in [16] reported that the window can be as small as 9x9 pixels, but the practical limit actually depends on the topography of the scene and the noise present in the signal. In the third method, the image is divided into windows of multiple sizes. As an example, windows of 8x8, 16x16 and 32x32 pixels are shown in Figure 3-5(d) and are marked as wl,w2,w3. The advantage of using this method is that the estimations over the large 32x32 pixel windows can be used to reduce the effect of the phase slope on the estimation over the smaller 16x16 and 8x8 pixel windows. Its main disadvantage is that the spectral estimation algorithm is not as effective as in other methods. For all methods there is always a trade-off between the size of the analysis window and the precision of the estimation. Moreover, the restriction of the signal to be locally stationary forces the window to become smaller, which degrades the precision of the estimation. In this sense robustness of the method to small windows is highly desired. Based on the assumptions made for each algorithm, but without a numerical evaluation of each method, it can be concluded that methods II and III have advantages over method I because they can split the image in smaller windows. In the second aspect - the way the local fringe frequency estimation is performed -method III has the disadvantage of not using an optimal estimation algorithm. But on the other hand, it is the method that better adapts the window size to the topographic changes of the scene. 39 a) b) v • c c c) d) Figure 3-5: Window splitting of different methods a) Input interferogram signal (ERS tandem mission data before flattening). b) 32x32 pixels window splitting used by the frequency estimation method based on the local periodogram. c) 16x16 pixels window splitting used by the two-dimensional MUSIC algorithm. d) 32x32, 16x16, and 8x8 pixels windows splitting used by the multi-grid method. Note: the background image in b), c) and d) corresponds to the magnitude of the instantaneous frequency for the input interferogram in a), using the method explained in next chapter. Darker shades indicate higher frequency magnitude. 40 Chapter 4 Local Fringe Frequency estimation - Proposed method The main limitation of the reviewed methods in Chapter Three is that they require a large number of samples {SijJk to estimate the frequency of the window A;, and this implies a reduction in the spatial resolution (localization) of the estimator. In order to get good estimates it is required to keep the window size large, but at the same time only small windows can properly track the topographic changes of the scene. In contrast, the proposed scheme pre-filters the signal before the local frequency estimation is performed, allowing a simpler and more efficient frequency estimator that operates at the pixel level. Of course, it is not possible to estimate the instantaneous frequency (LF) based on only one sample of the signal, so 1x5 pixel windows are used, implying an efficient frequency estimation. Because of their time-frequency localization - more specifically the achievement of the lower bound of the uncertainty principle - the filters used in the pre-filter stage are modulated Gaussian functions. In order to better adapt this pre-filtering scheme to the characteristics of low frequency and high frequency events, the filters have variable spatial aperture and variable frequency band-width, similar to the variable space-frequency domain filters used in Wavelet analysis. 41 The proposed approach is presented in the framework of the A M - F M (amplitude modulation - frequency modulation) signal representation model, which is introduced in Section 4.1. In Section 4.2, the proposed method is explained through an example using a one-dimensional signal, and in Section 4.3 the criteria to choose the multi-band filters is addressed. Finally, the application of the proposed scheme to two-dimensional signals is presented in Section 4.4. 4.1 AM-FM Signal Representation Model Since the topographic slope computed at every pixel in the ground elevation model produces an interferogram with one uniquely defined instantaneous frequency at each sample, only one frequency for each sample is necessary to represent the topographic information of the interferogram. Exceptions are the areas with shadows, layover and phase aliasing which the InSAR technique is not able to properly map and need to be identified and processed separately. Therefore, in the absence of noise, interferogram signals can be represented by locally monochromatic signals, i.e. they can be represented by one amplitude and one frequency for each pixel in the image: s(x) = A( x) exp( j Q(x)) x = (xl,x2) (1) Where A(x) > Ois the amplitude of the interferogram and V Q ( x ) = ( u(x),v(x)) ( dQ(x) dQ(x) ^ d xx ' d x2 is the gradient of its phase Q, (x). The phase function Q (3c) can not be directly and efficiently estimated because it is ambiguously embedded in the interval [-n,n), instead its gradient V£2(x ) can be efficiently estimated. The instantaneous absolute frequency / and the instantaneous orientation 0are defined as: f(x) = 4u\x) + v\x) Q(x) = t a n " 1 ^ (2) u (x) 42 where/and 8 are. related to the magnitude and orientation of the plane's slope, tangent to each point in the ground elevation model of the scene. By taking the derivative of s(x) respect to x in (1) the frequency is computed from the continuous signal as follows: VQ(x) = R e ( ^ ^ ) (3) And the amplitude as: A(x) = \s(x)\ (4) Under the assumption that A(x),VQ(jt) are locally smooth, a complete representation of the signal s(x) is the pair of estimates {A(x),V£2(i)} for all points x. The previous model constitutes the AM-FM (amplitude modulation - frequency modulation) mono-component representation of the signal s(x) [28]. In this representation, the input is modeled as a signal modulated in both amplitude and frequency. In [35], the authors established the smoothness conditions on A ( i ) , V Q ( i ) under which narrow-band channels filters effectively capture the local frequency structure of the signal s(x) under the model (1). For interferogram signals in the absence of noise, the smoothness conditions on the signal s(x), imply smoothness conditions on the scene's topography. Therefore, very abrupt topographic changes and discontinuities may not be completely captured. In general, for no-monochromatic images the one component assumption is not valid and (1) is replaced by: s(x) = X A(*)exp(M(*)) (5) 43 The AM-FM multi-component representation of the signal is the set of estimates { A , . ( i ) , V S ; ( i ) } for all points x [29]. If the signals are locally narrow-band the different frequency components can be separated using a filter bank that covers all the frequencies of interest. Although in the InSAR case there is only one frequency component of interest at each sample, it will be shown that this scheme is adequate to separate the noise components n(x) from the received signal: sn(x) = exp( jQ(x))+ n(x) (6) The objective is to estimate the phase gradient function. V Q ( x ) from the noisy samples sn(x). In addition, the phase function Q(x) can be recovered by integrating V£2(x) and s(x) reconstructed as in (1). 4.2 Proposed Fringe Frequency Estimation Method A discrete version of (3) could be applied directly to the input signal in order to estimate the instantaneous frequency (LF) at every pixel in the interferogram, but in the presence of noise these estimates will be highly degraded. Therefore, it is necessary to filter the signal before estimating its instantaneous frequency. Standard linear filtering techniques for noise reduction do not work in this case because the signal spectrum overlaps with the noise spectrum. Instead, the fact that the signal has only one dominant local frequency can be used to reduce noise by multi-band filtering. As an example, consider the ID synthetic phase interferogram signal of Figure 4-1(b). In this case x = x in equations (1) to (6) of the last section. Figure 4-l(a) shows a quadratic phase function £l(x) based on which the signal is constructed by doing s(x) = ejaM . Figure 4-l(c) illustrates its magnitude spectrum computed through FFT and Figure 4-1(d) its instantaneous frequency computed as the derivative of the unwrapped phase. The instantaneous frequency is linear within the interval [0,n]. 44 In general, if noise is added to the interferometric signal as in (6), it can not be entirely removed by standard linear filtering. However, based on the fact that the signal has only one dominant frequency at each sample, it can be locally band-pass filtered around its instantaneous frequency. Since the frequency at each specific sample is unknown for an arbitrary input, the signal is passed through a bank of band-pass filters like the ones illustrated in Figure 4-2(b), and the output of the filter with the highest response is used to estimate the instantaneous frequency. The processing steps are illustrated by the blocks: Multi-band filtering, MAX and Frequency Estimator in the diagram of Figure 4-3. The pre-filtering stage is named Multi-band filtering because the interval [0,n] is covered by the multiple pass-bands of the filters. The outputs signals yix) of the multi-band filter are M channels containing filtered versions of the input signal. After selecting the maximum channel output in block MAX, a discrete version of (3) is computed, and an estimation of the instantaneous frequency is obtained for all samples. The filters' choice can be arbitrary, the only requirement is that they sample the frequency space densely so the filter including the local fringe frequency has a maximum response and the other filters receive mainly noise. Because of their time-frequency localization - more specifically the achievement of the lower bound of equation (4) in Chapter 3 - the filters used in this work are modulated Gaussian functions of the following form: 1 -x2 12a2 ;'&»; x e e ' with Fourier Transform: for i = 1,2,..., M. (7) 45 where CO, is the center frequency of each filter, and a its spatial aperture parameter. These modulated Gaussian functions are called Gabor filters in the literature [22]. Figure 4-2(a) illustrates the real part of M Gabor functions modulated at C0„ = n C O m a x / ( M - 1 ) , (n = 0 ,1 , . . . ,M-1) covering the interval [0,n] for M = 6 , CO m a x = 2.845 . Figure 4-2(b) illustrates their FFT. The input signal is convolved with each Gabor filter in the multi-band filter scheme of Figure 4-3. The main lobes of neighboring filters intersect at 50% of its maximum value, so the filters densely sample the frequency domain and only one filter output is high while the others receive noise. The selection of the filters' parameters is discussed in Section 4.3. Figure 4-4 illustrates the real part of the outputs of the filter bank when the input signal is the one plotted in Figure 4-1(b). The solid line shows the real part of the signal and the dash-dot line its magnitude. Notice that the filter that has the maximum magnitude output at the beginning of the sequence is the one with the lowest center frequency , and moves towards the one with the highest center frequency at the end of the sequence. The output of the filter with maximum response is used to estimate the instantaneous frequency, in this way the noise out of the band of the specific filter being used is filtered out. Figure 4-5(a) illustrates the estimated frequency (continuous line) using a discrete version of the operator described by (3). Notice that the application of the multi-band filtering has an edge effect at the beginning and at the end of the sequence. The theoretical frequency - defined as the discrete derivative of the unwrapped phase - is illustrated in Figure 4r5(a) by the dash-dot line, the difference between the theoretical and the estimated frequencies is very small. Among the stages of the interferogram processing chain of Figure 2-4 where frequency estimates are needed, phase noise filtering is selected here for illustration purposes. The steps of the proposed method illustrated in Figure 4-3 are summarized as follows: 4 6 Figure 4-1: One-dimensional Synthetic interferogram signal a) Quadratic phase function of the signal £2(x)= K x 2 b) Real part of the Interferogram signal s (x )=e j n ( x ) Note: The amplitude of the signal is constant for all samples, misperception is caused by the use of discrete samples. c) Magnitude of Interferogram spectrum. d) Instantaneous frequency in radians computed as the derivative of £2(x). 47 Real part of Filter no.1 Real part of Filter no.2 40 50 60 70 80 Real part of Filter no.3 40 50 60 70 80 Real part of Filter no.4 40 50 60 70 80 Real part of Filter no.5 40 50 60 70 80 Real part of Filter no.6 Horizontal axes are signal sample number Figure 4-2(a): Real part of Gabor filters in the space domain Dashed lines illustrate the filters' magnitude, which is constant. Their modulating frequencies vary. S p e c t r u m o f G a b o r F i l t e r s Figure 4-2(b): Magnitude Spectrum of Gabor filters A l l filters have the same shape but are centred at different frequencies in order to cover the interval [0, n). 48 i) The input interferogram signal s(x) is convolved with each filter gn and the magnitudes of the filter outputs are compared for each pixel in the image. The comparison process occurs in block M A X where the channel y(x) with the highest magnitude is selected. ii) The instantaneous frequency (LF) is estimated on y(x) of the M A X output, giving the frequency estimate V Q ( x ) for each pixel x in the image. (The implemented IF estimator is the DESA-1 algorithm given by equation (19) in Chapter 5). iii) The estimate V Q ( x ) is used to perform slope-adaptive or band-pass filtering on the input interferogram (block A in Figure 4-3) and then phase unwrapping. Instead of phase unwrapping, the frequency estimates VQ(x) can be used to obtain the estimated unwrapped phase Q,(x) by integration (indicated by block B in Figure 4-3) and the signal sR(x) = ej^x) obtained. Notice that sR(x) is a smooth version of s(x), since its phase is reconstructed from frequencies estimated on filtered samples of s(x) . Figure 4-5(b) illustrates the reconstruction of the signal, based on the estimated frequencies illustrated in Figure 4-5(a). The continuous line is the real part of the reconstructed signal sR(x), and the dotted line the real part of the original signal s(x). The difference between the two signals is small and it is due to errors in the estimated frequency V£2(x). These frequency estimation errors are further discussed in Chapter 5. The consecution of sR(x) is not as important as the unwrapped phase VQ(x ) estimation itself, which gives the topographic information. Although alternative B in Figure 4.3 offers a way to obtain such estimation, only alternative A is considered in this thesis. 49 Interferogram signal s(x) MULTIBAND FILTERING: 9i 92 9M-I M A X J(x) INSTANTANEOUS FREQUENCY ESTIMATOR VQ (X) Estimated frequency] VQ (X) B: FREQUENCY INTEGRATION T s(x) Interferogram signal A: ADAPTIVE FILTERING Q (X) Integrated phase 1 Filtered Interferogram PHASE UNWRAPPING T Unwrapped phase Figure 4-3: Block diagram of the proposed method. The input signal is convolved with the multi-band filters g,{.) Then, the instantaneous frequency is estimated on the channel with maximum response. The estimated frequency can be used to perform A : slope-correction filtering on the input signal or B: to reconstruct the signal through frequency integration. 50 Filter Outputs Figure 4-4: Outputs of the multi-band filters of Figure 4-2 for the input signal in 4.1(b). The solid line shows the real part of the channel output signal for each filter. The channel with maximum magnitude response at each sample is used to estimate the instantaneous frequency. The signals magnitude is always indicated by the dash-dot line in spite of their plotting appearance. 51 a) Es t ima ted F requenc ies in rads 50 100 150 200 250 -> s a m p l e number b) R e a l part of Recons t ruc ted and original S igna l s a m p l e number Figure 4-5: Constant bandwidth filters results a: Estimated frequency on the maximum output of the filtered channels in Figure 4-4. The continuous line shows the estimated instantaneous frequency at each sample on the channel with maximum magnitude in Figure 4-4. The dash-dot line shows the theoretical instantaneous frequency computed as the derivative of the unwrapped phase in Figure 4-1(a). b: Signal reconstruction. The continuous line shows the reconstructed signal based on the integration of the estimated frequency of Figure 4-5 (a). The dotted line shows the original signal. Note: The amplitude of the signal is constant for all samples, misperception is caused by the use of discrete samples. 52-4.3 Filters Parameters Selection The filter's objective is to densely sample the frequency space of the interferogram signal in order to detect in which band the instantaneous fringe frequency is located. In the example of the previous section, Gabor filters of constant band-width were used to cover the interval [0, n). Narrow-band filters are desired because they reject more noise, but they also imply a higher number of filters and their sampling functions are more spread in space reducing the space localization property. On the other hand, very space-localized filters are more spread in frequency, according to the uncertainty principle of equation (4) in Chapter 3. For analysis purpo*ses, let the interferogram be divided in two kinds of events: low-frequency events due to smooth changes in the topographic elevation of the scene, and high-frequency events due to abrupt changes in the topographic elevation. Low-frequency events have larger space periods than high-frequency events. Therefore, they can be sampled with spatially wider Gabor filters (cr high). On the other hand, high-frequency events will require very localized filters (cr low), since these events are more localized in space. The desired filtering scheme should therefore have adaptive window sizes. 4.3.1 Wavelet-like Filters Scheme In order to adapt the multi-band filtering scheme to the characteristics of low-frequency and high-frequency events, the filters must have variable aperture parameter a. In this case, (7) becomes: 1 - x 2 / 2 o y / ' c o , * e e ' with Fourier Transform: 53 G,.(co) = e - a ^ ~ ^ 1 2 for i = 1,2,... ,M. (8) In order to maintain the product of the space and frequency resolutions A x - A c o , constant along all filters, the filter parameters C7-,Aco. must satisfy the following conditions: a M-i = r ' A M (9) A(0M_,. = for r>l, i = l,...,M-1 (10) r so that, a , _ ! A G O - , , = csi A c o , = constant, for i = 2,...,M-1. (11) For Gabor filters given by (8), it can be easily found that the frequency dispersion parameter A d ) , is proportional to the center frequency CO-. Therefore, the filters' center frequencies are also related by: <*M-i=^f for r > l , i = l , . . . , M - l (12) r In this way, given the parameters for one filter the parameters for the other filters are set such they satisfy conditions (9) and (12). If the parameter r is set to 2 or higher, the narrowest filter bandwidth Aojj will be too low and distortion on the instantaneous frequency estimator will occur. Therefore, the parameter r must be set to a value lower than 2. Figure 4-6(b) illustrates the filters frequency distribution for r=1.8, M=6. Notice that all filters frequency responses intersect at the same magnitude value (75%) but their bandwidths are not constant. The value of this threshold is inversely related to the amount of noise rejected. 54 4.3.2 Advantages of the Wavelet-like Filters Scheme. The multi-resolution feature of wavelet-like filters can be illustrated in terms of its space-frequency partition. Figure 4-7 illustrates the space-frequency partition for the constant cr configuration of Figure 4-2, the space-frequency plane is divided in cells of constant area and shape G X BW . The [0,n]frequency interval is divided in regular cells of value BW = 0.5712 rads/sample. The spatial aperture parameter is also constant a = Ami. Figure 4-8 illustrates the wavelet-like configuration of Figure 4-6, in this case the space-frequency plane is divided in cells of constant area - except for filter number 1 -but different shape, with dimensions a , X BWi. The filters parameters are listed in Table 4.1. For the low frequency filters the space windows are larger, and for the high frequency filters the space windows are smaller. In Figures 4-7 and 4-8, the smaller the frequency partition the narrower the filter's band-pass is, and consequently it will reduce more noise. From this point of view, constant band-width filters will have uniform performance along the frequency axis when stationary signals are considered. But, interferogram signals are not stationary, if a real frequency variation of short duration occurs within a long window, it will be identified as a disturbance in the constant bandwidth scheme. In Figure 4-7 the real part of a test signal is drawn, the signal's frequency is constant for all samples except for a short interval where it is increased to a higher frequency. In such interval the filter whose magnitude is maximum will distort the information contained by the signal since the analysis window is not short enough. In contrast, the wavelet-like configuration of Figure 4-8 will capture the non-stationary characteristics of the test signal by adapting the length of the analysis window. The interval where the signal changes its frequency will be detected by filter number 6 which has a narrower analysis window. In this way, local changes in the signal will be detected by the proper filter. 55 A relationship between the proposed scheme and the STFT method (reviewed in Section 3.1) has been established. Such relationship is explained in appendix B and shows the advantages of using analysis windows of variable size. On the other hand, condition (11) is equivalent to say that all filters have constant quality factor (Q). In [28] it is shown that for the non-linear energy operator used to estimate the instantaneous frequency (equation 13 in Chapter 5) uniform worst-case performance across the spectrum can only be attained by using this constant-Q or wavelet-like filter bank. This further motivates the use of the wavelet-like filtering scheme. -An additional advantage of using wavelet-like filters is that in the statistical distribution of typical interferograms, low instantaneous frequencies are more likely to occur than high instantaneous frequencies. Since the noise reduction capability is higher in narrower band-pass filters - for a given number of filters M - there will be more noise reduction if the wavelet-like filter distribution is used instead of the constant band-width filter distribution. Filters Centers (u)i) in rads/sample Spatial Aperture (<Ti) Bandwidth (BWj) in rads/sample 0.083 14.23 0.157 0.234 17.85 0.137 0.421 9.91 0.245 0.757 5.50 0.442 1.364 3.06 0.773 2.454 1.7 1.387 Table 4.1 Filters Parameters of Figure 4-6 Notes: The ratio between consecutive filters centers is 1.8. The listed bandwidths have been computed on rounded numbers. The parameters CDj and CM are as in (8), and obey conditions (9) and (12). 56 Real part of Filter no.1 Real part of Filter no.2 0.2 0.1 0 -0.1 -0.2 40 50 60 70 80 90 Real part of Filter no.5 ; / A \ 1 \ V' ; / 0.2 0.1 0 -0.1 -0.2 40 50 60 70 80 90 Real part of Filter no.4 40 50 60 70 80 90 Real part of Filter no.6 40 50 60 70 80 90 40 50 60 70 80 90 Horizontal axes are signal sample number Figure 4-6(a): Real part of Gabor wavelet-like filters in the space domain. 57 Constant BW Gabor filters 3h 2.5h Figure 4-7: Constant bandwidth space-frequency partition This figure illustrates the space-frequency partition for the constant bandwidth filters scheme of Figure 4-2. A frequency change of short duration is identified as a disturbance because the analysis window is not short enough. Note: The signal plot is interpolated between samples to facilitate interpretation. 58 Wavelet-like Gabor filters 3 BW-j Q | I I I ! I I I I ° 1 I I L 0 2 4 6 8 10 12 14 16 18 20 Space (sawples) -> Figure 4-8: Wavelet-like filters space-frequency partition This figure illustrates the space-frequency partition for the wavelet-like filters scheme of Figure 4-6. The interval where the signal changes its frequency is sampled by filter number 6 which has the narrowest analysis window. Note: The signal plots are interpolated between samples to facilitate interpretation. 59 4.3.3 Filters Parameters Selection Procedure Once the wavelet-like filters have been chosen because of the advantages mentioned above, a simple procedure has been designed in order to set the filters parameters. The procedure is as follows: Given: i) The number M of filters to be used which is proportional to the number of operations to perform on the input image. ii) The threshold p at which the filters' main-lobes intersect in the frequency domain. (This parameter is inversely related to the amount of noise rejected). Find the parameter r in (12) such that: The M filters cover the interval [0,n) with one of the filters output being higher than the normalized magnitude p for any frequency in the interval [0, n). By setting the frequencies such as in (12), they will only cover the interval (0,n) for M —> °°. Therefore, a limit has to be imposed to the series and the rest of the interval be filled with a filter covering the lower frequencies, as it is illustrated in Figure 4-6(b). Notice in Table 4.1 that the bandwidth of the two lowest frequencies filters is close. 4.3.4 Application of the Wavelet-like Filter Scheme The synthetic interferogram signal of Figure 4-1(b) is passed through the set of filters in Figure 4-6 producing the outputs illustrated in Figure 4-9. The maximum channel response is found for each sample, and the instantaneous frequency is estimated on that channel. Figure 4-10(a) illustrates the estimated frequency in continuous line, and the theoretical frequency -defined as the discrete derivative of the unwrapped phase - in dash-dot line. The difference between the theoretical frequency and the estimated frequency is small and is further discussed in Chapter 5. Figure 4-10(b) illustrates the reconstructed signal by integrating the estimated frequencies (continuous line), and the original signal illustrated by the dotted line. The difference between the two signals is small and it is due to errors in the estimated frequency V£2(x). The 60 reconstruction process in Figure 4-10(b) is started at a sample where the edge effect of the filters is not present anymore. 4.4 Application of the proposed method to 2D interferogram signals For two-dimensional signals, the same filters as in the one-dimensional case are used, but the filters are rotated each TI/8 radians in order to cover the two-dimensional frequency plane. Figure 4-11 (a) illustrates this case when the one-dimensional wavelet-like filters of previous section are used. Before an interferogram is flattened, negative frequencies of the signal (left side of the 2D frequency plane) correspond to areas of layover, as explained in Section 2.2. Any event on the left side of the two-dimensional frequency plane in Figure 4-11(a) will be filtered out, and any event on the right side of the plane will lay on the main lobe of one of the filters in the multi-band scheme (except at very high frequencies). The application of these filters could be enough to reduce layover areas but, since the filters along the vertical axis do not have sharp edges, half-band filtering must be applied in advance to reduce layover areas properly. The total number of filters in Figure 4-11 (a) is 9M, M being the number of filters used to cover the interval [0,7t]'m one dimension (M=6 in this case). The implementation of the multi-band filtering stage is done through FFT computations, which are more computational efficient than direct convolution. The FFT of all filters is computed in advance, and only the input FFT is computed on-line. On the other hand, the 9 filters close to the origin can be substituted for one common DC filter for efficiency purposes, but the center frequencies must be re-arranged to properly cover the interval [0-n). This reduces the number of filters in the scheme to 9(M-1)+1. Figure 4-11(b) illustrates the filters center frequencies in this case. After the instantaneous frequencies have been estimated, the frequencies must be corrected to eliminate the orbital frequency component. The interferometric phase in the absence of noise is modeled as a contribution of the orbital and topographic fringes: O(x) = Qorb(x) + Q t ( x ) (13) 61 where Qorb (x) = ^ ( ^ ) and Qtopo (x) = ^ ( ) . (as in (4) of Chapter 2). A d A d After the instantaneous frequency estimates V<f>(x) are obtained, the gradient of the orbital fringes phase model VQorb(x) must be subtracted in order to obtain the topographic fringe frequency estimates VQtopo(x). Usually Q,orb(x) is modeled as a linear phase trend in the range direction, toorbW = °>orb *1 ( 1 4 ) so the desired instantaneous topographic fringe frequency estimates are obtained by: VQtopo(x) = V 6 ( S ) - coorb (15) The procedure above simplifies the process of correcting the estimated frequency to just subtracting the orbital fringes frequency. 4.4.1 Application to Flattened Interferograms If the interferogram is first flattened, the resulting interferogram spectrum is shifted to the left in the range direction by an amount equal to the orbital fringes frequency C00rb. In Appendix C, it is shown that in the absence of noise, equal performance of the instantaneous frequency estimator is achieved for all bands of the multi-band filtering stage. Therefore, shifting the spectra does not have any influence in the instantaneous frequency estimation, as long as the multi-band filters match the signal's spectrum location. However, in the presence of noise (which is the practical case), there is an advantage when the signal's dominant components are close to base-band. In this case, the low frequency filters of the multi-band scheme are more likely to contain the dominant frequency, and since they have narrower bandwidths, more noise is filtered out. This is partially fulfilled if the interferogram is first flattened and then symmetric filters are created to cover the negative frequencies interval. 62 Unfortunately, this implies the creation of new filters (and their FFT computation) for each new value of 0Jorb. The use of a pre-existing D E M (digital elevation model) has been shown to help the process of obtaining accurate elevation maps from interferograms. Particularly, if a pre-existing D E M is used to locally flatten the interferogram, its spectra will be centered at the origin in both range and azimuth directions [3]. Moreover, if the input D E M had no errors, all the energy in the bandwidth of the original interferogram would be concentrated at zero frequency. In the case of interferograms flattened with a pre-existing D E M , the filters' frequency response can be located within the frequency sub-plane - n/2 < (x)x < Tl/2, -Tl/2 < C0y < 7l/2 as illustrated in Figure4-12. The total number of filters is: 16(M-1)+1 =49 with M=4. The filters' distribution parameters of Figure 4-12 are listed in Table 4.2. 4.4.2 Practical Procedure In this work, the wavelet-like filter distribution was chosen because of the advantages mentioned in Section 4.3.2. Regarding the filters' parameters, they are not unique and can be optimized for a specific type of interferogram depending on the scene's topography statistics. However, such optimization is not addressed in this thesis. In [34], the filters' frequency distribution of Figure 4-11(b) was used to evaluate the performance of the proposed frequency estimation algorithm in synthetic interferograms before flattening, with satisfactory results. These filters sample all the range positive frequencies of the spectra. Interferograms before or after flattening can be demodulated in the range direction in order to locate the most significant part of the spectra at base-band. Then, the filters' frequency distribution of Figure 4-12 can be used in order to achieve higher noise reduction. The estimated frequencies have to be subsequently adjusted by subtracting the frequency value used to demodulate the signal. If the spectrum of a flattened interferogram (without DEM) happens to be concentrated at low frequencies, the use of the filters' frequency distribution of Figure 4-12 without demodulation produces very similar results. An example is shown in Section 7.1. 63 Should a D E M be available for flattening the phase, the filters' frequency distribution of Figure 4-12 is the one to be used because it matches the spectra and achieves higher noise reduction. On the other hand, the D E M used to flatten the interferogram also reduces the complexity of the phase unwrapping operation by reducing the variability in the interferogram phase. In conclusion, the filters' frequency distribution of Figure 4-11(b) is the one to use for any interferogram before flattening. The symmetric filters' frequency distribution of Figure 4-12 is preferred because is achieves higher noise rejection as explained in Section 4.4.1. For interferograms flattened with a D E M , the signal's spectrum is located at base-band so the symmetric filters' frequency distribution is used. When no D E M is available, demodulating the interferogram signal in the range direction provides a method to match the signal's spectra and the symmetric filters' distribution. The main interest in next chapters is to evaluate the performance of the proposed method, using the symmetric filters' distribution of Figure 4-12 and the parameters listed in Table 4.2. Filters Centers (cDj) in rads/sample Spatial Aperture (Oi) Bandwidth (BWj) in rads/sample 0 5.20 0.225 0.326 11.55 0.205 0.620 6.08 0.380 1.178 3.20 0.740 Table 4.2 Filters Parameters of Figure 4-12 Notes: The filters' magnitude responses intersect at 50% of its peak in the radial direction. The ratio between consecutive filters centers is 1.8. The listed bandwidths have been computed on rounded numbers. The 2D frequency filters are obtained by rotating the one-dimensional filters each 22.5°. 64 Filter Outputs Figure 4-9: Outputs of the multi-band filters of Figure 4-6 for the input signal in Figure 4-1 (b) The continuous line shows the real part of the channel output signal for each filter. The channel with maximum magnitude response at each sample is used to estimate the instantaneous frequency. The signals magnitude is always indicated by the dash-dot line in spite of their plotted appearance. 65 a) Es t ima ted F requenc ies in rads -> s a m p l e number b) R e a l part of Recons t ruc ted and or iginal S igna l s a m p l e number Figure 4-10: Wavelet-like filters results a: Estimated frequency on the maximum output of the filtered channels in Figure 4-9. The continuous line shows the estimated instantaneous frequency at each sample on the channel with maximum magnitude in Figure 4-9. The dash and dot shows the theoretical instantaneous frequency computed as the derivative of the unwrapped phase signal in Figure 4-l(a). b: Signal reconstruction. The continuous line shows the reconstructed signal based on the integration of the estimated frequency of Figure 4-10(a). The dotted line shows the original signal. Note: The amplitude of the signal is constant for all samples, misperception is caused by the use of discrete samples. 66 Filters center frequencies Oo, oo' o :oo 0 o o -3 -2 -1 0 1 2 3 4 5 x spatial frequency —> Figure 4-ll(a): Center frequencies of the filters for the two-dimensional case. Filters center frequencies - with common DC filter CL 10 >-1 -2 -3 c ?0 c J o 0 ° 0 o o o o o£ 0°o 0 o o c .0. 0 ° 0 0 0 0 1 2 x spatial frequency —> Figure 4-ll(b): Center frequencies of the filters for the two-dimensional case Substitution of the 9 filters closest to the origin in Figure 4-11(a) by one D C filter. Note: For both figures, the x and y axes correspond to frequencies in the range and azimuth directions, respectively. 67 GABOR FILTER'S CENTER FREQUENCIES 1.5h c o o o o c 3 o o o o o o o o o o o ^ o o o o o O 0 O o o o o o o o o o o o o o o o o ° o ° o o o 6 0.5 •3 -0.5 -1.5 h--1.5 -1 -0.5 0 0.5 1 X spatial frequency > rads/sample 1.5 Figure 4-12: Symmetric filters distribution for interferograms flattened with a pre-existing D E M . Note: The x and y axes correspond to frequencies in the range and azimuth directions respectively. 68 Chapter 5 Instantaneous Frequency Estimation In the scheme proposed in Chapter Four, it was required to estimate the instantaneous frequency of the interferogram signal after filtering it and selecting the channel with maximum magnitude. The assumption made is that the signal power is strong enough compared to the noise power, so that the multi-band filtering stage is able to properly find in which band the instantaneous frequency is located. By selecting that specific channel, the noise out of that band is filtered out and then, the instantaneous frequency (IF) of the filtered signal is estimated. The objective of this chapter is to present the IF estimation algorithm to be used in the proposed scheme. The definition of instantaneous frequency is given in Section 5.1, and the presentation of the selected instantaneous frequency estimator - the DESA-1 algorithm - is made in Section 5.2. Both Sections 5.1 and 5.2 constitute prior work. An intensive evaluation of the DESA-1 algorithm has not been reported. In order to investigate the adequacy of the algorithm for the SAR interferometry application, an evaluation of the DESA-1 algorithm was performed under different conditions of the input signal. The details are documented in Appendix C and only the conclusions are cited in this chapter. 5.1 Instantaneous Frequency Estimation The instantaneous frequency (IF) of a monochromatic continuous signal s(x) = AejQ{x) is defined as the derivative of its phase: 69 = 1 dQ(x) ( 1 ) 2TT dx If the signal is discrete, a discrete finite impulse response (FIR) differentiator can be used to achieve the differentiation operator. Then, the discrete IF can be defined as: / , (n) = ^-Q(n)*d(n) (2) 2,71 where, d(n) is the impulse response of a FIR differentiating filter, and * denotes convolution. It has been found that the FTR filter d(n) exaggerates the effects of high frequency noise [30], so an approximation to the differentiation operator is used instead. i p For polynomial phase signals of any order p, Q (n) = ^ akn whose IF is: k=0 2n k = ] A q-th order generalized phase difference estimator can be defined as: 1 i / 2 / ( » ) = T ~ £ b^in + l) (4) 271 l=~q/2 with q an even integer. The coefficients bt are found so that f^ri) = f(ri), i.e. 1 P 1 9 / 2 — l k a k n k - 1 = — X * / Q ( » + 0 (5) Z7T /=_ 9 /2 70 For example for p=2, the phase difference estimator is the so-called central finite difference operator. 1 (Q(n + 1 ) -Q(n -1)) (6) 471 This estimator is unbiased for linear F M signals, but in the presence of noise the signal is not purely monochromatic and the phase signal can not be extracted without distortion. Although the noise present in the signal is reduced by the multi-band filtering stage, there will be some noise within the selected band. In the presence of additive noise the signal model is: Where v(n) corresponds to the additive part of the noise. In the presence of noise, the extracted phase will have noise contributions and the operator given by (6) will exhibit very high variance. In order to reduce the variance of the estimator, the phase difference can be smoothed, but it will produce biased estimates. There are many IF estimators that could be used for interferogram signals. In [30] a review was made where the following IF estimators among others were analyzed: the Estimator Phase difference given by (6), the Peak of the Short Time Fourier Transform (STFT) - explained in Chapter 3, Section 3.1.2 - , and the Peak of X W V D (Cross Wigner-Ville Distribution). Among all IF estimators mentioned above, the peak of the X V W D has the best performance when considering additive noise to a F M (Frequency Modulated) signal. The Cross Wigner-Ville Distribution (XVWD) between a reference signal s(n) and an observing signal r(n) is defined as: s(n) = eKQ{n)) + v{n) (7) M 12 Wsr{n,k) = X s(n + m) r*(n-m) e j2nmkl M (8) m=-M 12 71 In this method the reference signal is estimated from the observed signal by obtaining an initial IF, any estimator can be used to estimate the starting IPs. The X W V D procedure is as follows: i) Form a unit amplitude reference signal from the estimated IF. ii) Form the X W V D between the reference and the observed signal and estimate the new IF from the peak of the X W V D . Repeat steps i) and ii) until the difference of IF estimates from successive iterations is less than a specified amount. The author in [30] evaluated the methods using a linear F M signal, with IF varying from 10 to 90 Hz, sampled at 200 Hz. The IF at each sample was calculated using 128 samples window, which is extremely large for the InSAR application. The X W V D method gave the best performance among the methods that meet the Cramer-Rao (CR) bound [30]. However, the high computational complexity of the X W V D method constitutes its main disadvantage. In order to implement the scheme proposed in Chapter 4, a simple and efficient IF estimator is considered here based on the energy signal analysis ESA operators introduced in [31]. 5.2 Energy Separation Algorithm The energy separation algorithm (ESA) was introduced to analyze speech signals from the point of view of the energy required to generate them. Considering a linear un-damped oscillator consisting of a mass m and a spring of constant k, its displacement x(t) is governed by the motion d2x equation m—=- + kx — 0. Whose general solution has the form: xii) = Acos(coot + 0) (9) The instantaneous Energy E0 of the system is the sum of the kinetic and potential energies and for the un-damped case it is constant: E0 = 2X •T ( A a ,- ) 72 The Energy tracking operator l//c is defined by: (10) which when applied to equation (9) produces: yrc[ Acos(cvot + 0)] = A2co02 = E. o (11) m / 2 Therefore, the output of y/c is equal the energy (per half unit mass) of the source producing the oscillation. If x(n) is a sampled version of the continuous-time signal x(t), and the derivative is replaced with the 2-sample backward difference: (x(ri) - x(n -1)). Then, the discrete version of the energy operator lf/c is: When (12) is applied to a discrete signal with constant A amplitude and frequency Q,c, it information about the frequency and amplitude is contained in the expression. The constant frequency and amplitude can be found from the following equations [31]: y/[x(n)] = x2(n) - x(n +1) x(n -1) (12) produces: lf/[Acos(Qcn + 0)] = A s i n 2 ( Q ( . ) , which is similar to (11) in the sense that the (13) 73 Considering now the application of the operator (12) to an F M signal of the form: x(n) = A cos [Qcn+ 0 +Q,m\q{m)dm] (14) o n where (j)(n) = Q,cn + 6 + Q m jq(m) dm is the phase of the signal. Its instantaneous frequency o is defined by: = Qt(n) = Qe + Qmq(n) (15) dn with the condition that the instantaneous frequency Q ; (n )G [0,7T), which is equivalent to 0 < Q c - Q m < Q( in) < Qc + Q m < n, and implies | q(rri) |< 1 and Qm < Qc. The result of applying the operator (12) to the F M signal in (14) is: Wix(n)] = A2 y/[cos(tp(n))]. (16) For slow-varying Q,-(n) compared to Q c , y/[cos(ip(n))] ~ sin 2 (£2,. (n)) , and expression (16) becomes: y/[x{n)] « A2sin2(Q(-(n)) (17) On the other hand, the application of the operator (12) to the backward difference is: y/[x(n) - x(n -1)] - 4a2(n) s i n 2 ( ^ i M ) s i n 2 ( Q . ( n ) ) . From which the instantaneous frequency can be obtained by: Q | ( N ) , c o s - { l - ^ U ( n ) - ^ - 1 ) ] } (18) 74 A further improvement results if the term y/[x(ri) - x(n -1)] is made symmetrical by averaging it with y/[x(n +1) - x(n)], thus repeating the procedure yields another algorithm given by the expression: Q ( ( B ) . e o .-.{i_rt*W-M»+i)]} 1 m m ] ' This is called DESA-1 (discrete energy separation algorithm), and it is chosen to implement the IF estimation in the proposed method scheme of Figure 4-3. The DESA-1 algorithm has a low computational complexity, only five neighboring samples are required to compute (19). The DESA-1 algorithm can estimate frequencies up to one half of the sampling frequency of the signal. However, if the noise present in the signal makes the argument of cos _ 1(v) such that v £ [0,7t), the IF is set to the value of the previous sample - assuming that the IF varies slowly. On the other hand, if l//[x(n0)] = 0 and assuming Q(no) ^ 0 the frequency is also set to the previous value. As an example, consider the application of the DESA-1 algorithm to the F M signal of Figure 5-l(a). The input sequence is: x(n) = exp(j<p(n)) = exp( j[ £lc n + Af s i n ( ^ m o d n ) ]) (20) Its instantaneous frequency is defined by: ^ = Q,(n) = R ; + A / ^ m o d c o s ( Q m o d n ) (21) dn With parameters £ 2 m o d =7r/100, Q C =3^ /8 and Af-7l/S, the instantaneous frequency is bounded to the interval [ ^ / 4 , ^ / 2 ] . Figure 5-l(b) shows the estimated IF and Figure 5-l(c) shows the absolute value of the normalized estimation error defined as: 75 8(n) ^ Q ' ( n ) - Q ' ( n ) x l O O % (22) Q,.(n) where Q.(n) is the estimated IF and Q-(n) is the theoretical IF derived from (21). The maximum absolute value of the error is 0.62% which is very small and can be neglected in practical terms. 5.3 DESA-1 Algorithm Performance In Appendix C, the performance of the DESA-1 IF estimator given by (19) is obtained by using the F M test signal (20) with varying parameters and different levels of noise. The results show that the RMS estimation error is considerably lower when the signal's IF changes slowly, and it is moderately high (23%) when the IF changes fast. The performance curves in Appendix C, show that the DESA-1 algorithm performance is adequate for IF estimation in interferogram signals, when used in conjunction with the multi-band filtering scheme. The disadvantage of the DESA-1 algorithm, of not following very fast changes of the IF can be solved by over-sampling the input data, but this increases the computational complexity of the algorithm. The proposed method is not, restricted to the use of the DESA-1 estimator. In general, any IF estimator will benefit by receiving a less noisy input signal. Therefore, there is a need for research on alternative IF estimators, which can then be used in the proposed scheme. In the next chapter the performance of the scheme proposed in Chapter Four using the DESA-1 algorithm given by (19), is evaluated. 76 a) Real part of Input F M signal i i i i i i i I 50 100 150 200 250 300 350 400 Sample number —> Figure 5-1: Example of IF estimation by the DESA-1 algorithm a) Real part of the signal: x(n) = exp(;[Q c« + Af sin(Omod/z)]) with parameters Q m o d = n/W0, Q ( ; = 3TZ78-b) Estimated frequency using D E S A - 1 , Q.e[nlA,nl2] c) Normalized error. The maximum error achieved ( 0.62% ) is negligible. 77 Chapter 6 Results on Synthetic Interferograms In this chapter, results are presented using the local fringe frequency estimation method proposed in Chapter Four. Synthetic interferograms with different spatial distribution and coherence are used to simulate real interferograms. The model used to generate them is presented in Section 6.1, followed by the error performance evaluation. In order to compare the proposed method against other window-based frequency estimation methods, a hypothetical ideal estimator is introduced in Section 6.1.2. Its performance is obtained and conclusions are made based on their comparison. One of the applications of frequency estimation in interferograms is phase noise filtering. In Section 6.2 two existing methods to perform such filtering are presented together with a proposed method based on the IF estimates. A comparison between them is done using synthetic interferograms. 6.1 Frequency Estimation Results In order to illustrate and evaluate the use of the proposed method in the estimation of instantaneous frequency (IF) for interferogram signals, consider a "Gaussian hill" synthetic phase signal (p(x, y): 78 with x, y = 1,...,N. The synthetic interferogram formed by two SLC (Single Look Complex) images 5j, s2 is obtained using the following model: [7] s(x,y) = sls2 where nx,n2 are independent complex random variables with normal distribution, zero mean and variance <Jn . The correlated part common to both signals c, is modeled by a 2 complex random variable with normal distribution, zero mean and variance cr. . Using this model the interferogram coherence is given by: [7] The phase pattern is plotted in Figure 6-1(a) for a=70 and the magnitude of its theoretical IF is given in Figure 6-1(b). The real part of the interferogram signal is shown coherence value of 0.44. Figure 6-l(d) illustrates the 2D FFT magnitude of the interferogram signal. After convolving the input signal with each one of the Gabor filters illustrated in Figure 4-12, the filter output with maximum magnitude is selected. As an example of the filtering properties of the proposed scheme, Figures 6-2(a) to (c) give the real part of three of the Gabor filters' outputs. Figure 6-2(a) illustrates the output of the DC Gabor filter. Figure 6-2(b) illustrates the output of the Gabor filter centered at the second radial frequency in the horizontal direction and Figure 6-2(c) of the Gabor filter centered at the third radial frequency in the horizontal direction as well. The filter response is of high 5, (x, y) = ce1<p + ft, and s2 ( x, y ) = c + n2 (2) (3) in Figure 6-l(c) where the noise variance is o~n2=1.25 and ac2=l, which implies a 79 a) Syn. Interferogram Phase b) Magnitude of Theoretical Freq. Figure 6-1: Synthetic interferogram example. This figure illustrates the effects of noise and frequency slope on the estimation accuracy. a) Phase pattern for a = 70 b) Magnitude of theoretical IF c) Real part of the interferogram signal - coherence value = 0.44. d) 2D FFT magnitude of the interferogram signal. Note: In plot 6-l(b) the frequency has been multiplied by 12 for better visualization. 80 contrast (high magnitude) in the areas where the frequency and orientation of the input signal correspond to that of the filter. Notice that the outputs are smooth in the areas where the magnitude is high, since the out-of-band noise has been filtered out in these regions. Comparing Figures 6-2(a), (b) and (c), it can be observed that they complement each other. For instance, at the center of the image, the response of the filter in 6-2(a) is high and it is low in the other two filter outputs. For higher frequencies, the filter response is low for the DC filter and it is high for the second and third radial frequencies - Figures 6-2(b) and(c). Figure 6-2: Outputs of three of the Gabor filters. The center frequency of the channel with the maximum response produces a piece-wise approximation of the interferogram IF at every pixel, as illustrated in Figures 6-3(a) and (b). The estimated IPs on the channel with maximum response using the DESA-1 algorithm of equation (19) in Chapter Five are illustrated in Figures 6-3(c) and (d) for the range and azimuth directions respectively. In order to illustrate the accuracy of the IF estimation in the presence of noise, Figures 6-4(a)-(d) show one-dimensional (ID) slices of the IPs. The ID slices are done by the center of the image in the range and azimuth 81 directions of Figures 6-3(c) and (d). The estimated IF is drawn with a continuous line and the theoretical JJF with a dash-dot line. For the purposes of illustrating the specific channel used at every pixel, the center frequency of the filter with maximum output used to estimate the IF is drawn with a dotted line. Note that the centre frequencies of the filters are heavily quantized, but the estimated frequencies are continuous. The range and azimuth frequency estimations are made independently and they are stored in a complex number as: IF(x,y) = IFrange(x,y) + j IFazimuth(x,y) (4) In this way, the magnitude | IF(x, y) | and angle tan - 1 I F ^ m M ^ x , y > > 0 f the IF -1Fran8e(X>y) which are geometrically related to the topographic slope magnitude and angle - can be easily computed. 6.1.1 IF Estimation Error Performance There are two sources of errors in the IF estimation of the synthetic interferograms: • The noise added to the synthetic signal in the model given by (2). • The distortion error inherent to the DESA-1 algorithm. This error is present even for noise-free signals and it is due to the inaccuracy of the IF estimator. It is particularly high when the signal's IF changes very fast (see Appendix C). To separate both sources of error, consider the IF estimation results on the signal without noise. In Figure 6-5, ID slices of the estimated IF are drawn. Notice that there is ripple in the estimation of the IF of Figure 6-5(a) and (d), and there are jumps in the estimated frequency of Figures 6-5(b) y (c) in the areas where the filter center frequency 82 a) Filter Center Range Frequencies b) Filter Center Azimuth Frequencies c) Estimated IF Range Frequencies d) Estimated IF Azimuth Frequencies Figure 6-3: Output of the frequency estimation process: a) - b) Filter center frequency of filter with largest response, c) - d) Estimated frequency of each pixel in the interferogram. 83 Figure 6-4: ID slices of estimated local frequency The estimated IF is given by the continuous line and the theoretical IF by the dash-dot line. Filter center frequencies are given by the dotted line. 84 Figure 6-5: ID slices of estimated IF when no noise, is present in the signal The estimated IF is given by the continuous line and the theoretical IF by the dash-dot line. Filter center frequencies are given by the dotted line. 85 changes. This distortion error will increase with higher variations of the input IF which can be simulated by increasing the slope of the Gaussian hill, i.e. decreasing am (1). 2 The interferogram coherence is related to the noise variance by (3), <7C is kept constant and equal to 1.0. Figure 6-6 illustrates the M S E (Mean Square Error) of the estimation vs. interferogram coherence, when the values of a in equation (1) are set between 20 (simulating steeper topography) and 90 (flatter topography). The M S E error in Figure 6-6 is defined by: £MSE= ~ ^realin) f & where Q,is the estimated IF, Q r e a / i s the real IF and TV the number of samples in the image. In Figure 6.6 the estimation error is low for coherence values higher than 0.5. For slow IF variations (high a) the error is lower than for fast variations, because it is easier for the DESA-1 algorithm to track slower frequency changes as it was found in Chapter Five. When the coherence of the simulated interferogram drops below 0.57, the M S E error increases dramatically for values of a higher than 20. At this value of coherence, the multiband filtering and the highest magnitude filter output selection are not able to locate the local frequency due to the larger phase noise. Based on these simulations it can be concluded that the proposed IF estimation scheme works well for interferograms with coherence higher than 0.57. 6.1.2 Comparison Against Window-based Methods. A complete comparison between the proposed method and other window-based frequency estimation methods requires their implementation, in order to obtain 86 IF Estimation performance 0.3 0.4 0.5 0.6 0.7 0.8 0.9 .1 1.1 Interferogram Cohe rence —> Figure 6-6: Estimation error of the proposed method for different values of a. Lower values of a simulate steeper topography and higher values flatter topography. 87 a) Syn. Interferogram Phase b) Magnitude of Theoretical Freq. Figure 6-7: Synthetic interferogram example with rough topography (a=20) This figure illustrates the localization properties of the proposed method on a synthetic interferogram with rough topography (a=20). Note: In plot 6-7(b) the frequency in the vertical axis has been multiplied by 12 for better visualization. 88 performance curves such as in Figure 6-8. Instead, the introduction of a hypothetical estimator in this section serves as an ideal model for window-based frequency estimators. The main advantage of the proposed method is its high space localization. To illustrate the localization properties of the proposed IF estimator consider the input synthetic interferogram illustrated in Figure 6-7 for a =20 and no noise. In this case, the signal has rapid spatial variations and the content of the signal is spread in many more filters than in the previous case of a =70 - see the signal's 2D FFT in Figure 6-7(d). A ID slice of the estimated IF is illustrated in Figure 6-8 by the continuous line and the theoretical IF by the dash-dot line. In order to compare the performance of the proposed method against other window-based IF estimators, consider a hypothetical IF estimator that exactly computes the average of the noise-free signal's IF in cells of A ; 2 pixels. To implement such an estimator, the theoretical frequency is averaged in cells of A j 2 pixels and the resulting average is assigned as the IF of each one of the pixels in the window. Such a hypothetical IF estimator is drawn in Figure 6-8 by the dashed lines for A=15. Comparing the proposed IF estimator against the hypothetical one in Figure 6-8 for the noise-free case, it is observed that the proposed method achieves a better performance due to its localization properties. Figures 6-9(a) to (d) show the estimation error of the hypothetical estimator in dashed lines, superimposed onto the performance curves of Figure 6-6. Different values of A -namely 5,10 and 20 pixels - have been used to compute the average IF in the hypothetical estimator. For each experiment (value of a), the IF estimation error of the proposed method (drawn in continuous line) is always less than that obtained by the hypothetical estimator when the averaging window of the hypothetical estimator is 10 or 20 pixels, and the coherence higher than 0.6. The estimation error is approximately equal for both estimators when Ai equals 5 pixels and the coherence is 0.66. It is important to keep in mind that the number of pixels used in the IF estimation by the DESA-1 algorithm is also 5 pixels. 89 Fast IF variation 20 40 60 80 100 120 sample number —> Figure 6-8: ID slice of estimated IF when no noise is added to the signal and a = 20 The estimated IF is given by the continuous line and the theoretical IF by the dash-dot line. The dashed line shows the IF estimated by a hypothetical estimator, which computes the average of the theoretical frequency in non-overlapping windows of 15 pixels. 90 1 0 A 1 0 " I L U CO 1 0 " 1 0 " A 1 0 ' I L U CO y= io-3 1 0 " a) a=30 1 0 " A 2 0 pixels A 1 0 " I \ 1 0 pixels L U CO y=io 5 pjxejs _ O 1 0 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1 Coherence -> c) a=70 1 0 \ 2 0 pixels \ 1 0 pixels IF MSE -> o o j k — 5 pixels b) a=50 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1 Coherence -> d) a=90 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1 Coherence -> 0 . 5 0 . 6 0 . 7 0 . 8 ' 0 . 9 1 Coherence -> Figure 6-9: Performance of proposed method vs. Hypothetical Estimator The performance of the proposed method (in solid line) is shown for different values of or (topography roughness). The performance of the hypothetical estimator is shown in dashed lines for three different window sizes: 10, 20 and 30 pixels. Note: The plotting scales in all cells are constant. The vertical interval plotted is varied for better visualization. 91 Since the estimated IF of the hypothetical method is computed in the absence of noise, a real non-overlapping window-based IF estimator will have estimation errors higher than that plotted in the dashed lines in Figure 6-9. Therefore, it can be concluded that the error of the IF estimation using the proposed method is less than that obtained by non-overlapping window-based methods for coherence values higher than 0.6. For overlapping window-based methods using 10 pixels or more with 50% overlapping, the previous conclusion is held for coherence values higher than 0.66. 6.2 Phase Noise Filtering Results Phase noise filtering is one of the important applications of IF estimation in interferogram processing (see Section 1.1). In this section, the performance of two existing methods for phase noise filtering is compared against the performance of a new method that uses the estimated IFs. The main difference between the three methods is the assumption they make about the local topography, as explained below. (More detailed descriptions of these filtering methods are found in Appendix D). • Standard multilook averaging: This method performs standard samples averaging to estimate the phase. It assumes that the phase is locally constant. This method reduces the phase noise variance at the expenses of losing resolution [7]. • Slope correction filter: This method models the phase linearly to correct the effect of the topographic slope. Then, it performs standard multilook averaging, reducing the noise variance with less distortion [6]. • Non-linear phase modeling filter (new method): This method makes use of the IF estimates given by the proposed IF estimation method. It builds a local non-linear model of the phase by integrating the estimated IPs, then performs multilook averaging after correcting the effect of the local topography in the phase. 92 6.2.1 Comparison between Phase Noise Filtering Methods Using the synthetic interferogram model of Section 6.1, interferograms of different coherence and phase spatial distribution are generated in order to evaluate the proposed phase noise filtering method. The advantages of the non-linear phase modeling filter over standard multi-look averaging are pronounced when the input IFs vary rapidly. Consider, for example, the noise-free interferogram of Figure 6-10(a). Figure 6-10(b) shows the result of applying the non-linear phase modeling filter. In addition, Figure 6-10(c) shows the result of applying the slope-correction filter (linear phase modeling) and Figure 6-10(d) shows the result of applying standard multi-look averaging. A l l filters used window sizes of 11x11 pixels. Visually, the image in 6-10(d) is the most distorted. To further visualize their differences, the filtered phases have been unwrapped in one dimension across the center of the image in the range direction. They are illustrated in Figures 6-11(a). The plots show the unwrapped phase after slope-correction filtering in dash-dot lines. In addition, the unwrapped phase after non-linear phase modeling filter is shown in continuous lines and the unwrapped phase after standard multi-look averaging in dashed lines. The phase difference errors are shown in Figure 6-11(b), in which the distortion caused by the standard multi-look averaging is much higher than that caused by the other two filters. The distortion is caused by averaging pixels that do not have the same original phase. The phase error differences for the other two methods are close. However, for the non-linear phase model filter the IFs used were estimated on the noisy input, in contrast with the noise-free hypothetical frequency estimates used in the slope correction filter. This hypothetical estimator - as explained in Section 6.1 - computes the window average of the theoretical frequency on noise-free data. Figure 6-12(a) illustrates the real part (range direction) of the input theoretical frequency for the interferogram illustrated in Figure 6-10(a). Figure 6-12(b) illustrates the real part of the average frequency used to build the linear phase model in the slope-correction filter, which corresponds to the output of the hypothetical estimator for 11x11 pixel windows. Figure 93 6-12(c) shows the real part of the estimated frequency used to build the non-linear phase model in the proposed method. Figures 6-13(a)-(d) illustrate the case when the interferogram coherence is 0.5. The image in Figure 6-13(d) - corresponding to standard multi-look averaging - is the most distorted. The real part of the estimated frequency used to build the non-linear phase model is illustrated in Figure 6-12(d). These frequencies are estimated on the noisy data, whereas the frequencies used to build the linear phase model are the window averages of the theoretical frequency shown in Figure 6-12(b). The unwrapped phases are illustrated in Figure 6-14(a) and the phase error differences are shown in Figure 6-14(b). Again, the distortion caused by multi-look averaging is much higher than the caused by the slope-correction filter or the non-linear phase modeling filter. In order to illustrate how the phase estimation error behaves for different coherence values, the M S E of the differences between the unwrapped filtered phase and the theoretical one have been computed for the synthetic interferogram of Figure 6-13(a). Figure 6-15(a) and (b) illustrate the M S E of the phase difference error when 5 and 11 pixel windows are used respectively. The interferogram coherence is varied between 0.36 and 1.0 using the model presented in Section 6.1. In Figure 6-15(a) (where five pixel windows have been used), all filters behave approximately equal for high coherence values (>0.8), achieving small errors (<102 rads.). When the coherence is lower than 0.8, the standard multi-look filter error increases as depicted by the square symbols. The slope correction filter (star markers) and the non-linear modeling filter (circle markers) continue behaving similarly, but when the coherence is lower than 0.5, the threshold where the EF estimation method fails is reached. In conclusion, for five pixel windows the error achieved for all methods is similar if the interferogram coherence is higher than 0.8. This interval is indicated by the letter A in Figure 6-15(a). Because the multi-look filter is the simplest filter, it is the one to use in such cases. For interferogram coherence between 0.5 and 0.8, both slope correction 94 filtering and non-linear phase modeling filtering produce low errors. For low interferogram coherence (<0.5), the use of the proposed filter is not adequate. On the other hand, when the window size is increased to eleven pixels - Figure 6-15(b) - the error obtained by standard multi-look averaging (square markers) and slope correction filtering (star markers) is higher than that obtained when five pixel windows are used for the same coherence values. The reason for this is that the averaging window size is larger than the spatial variation of the input phase pattern. In contrast, for the non-linear phase modeling filter (circle markers), the error performance using eleven pixels is lower than the error achieved using five pixel windows, when the interferogram coherence is less than 0.8. The reason for this error reduction, is that non-linear phase modeling better tracks the phase input changes and by averaging more samples it reduces more noise variance. Moreover, the non-linear phase modeling filter achieves the lowest error of all three methods when the coherence is within the interval [0.5, 0.8] and 11 pixel windows are used. This interval is indicated by the letter B in Figure 6-15(b). For interferogram coherence less than 0.5, the proposed IF estimation method fails and slope-correction filtering achieves the lowest error. The estimated frequencies used in the slope-correction filter method are given by the hypothetical ideal estimator using 5x5 and 11x11 pixel averages respectively, whose behavior does not change in the presence of noise (see Figures 6-9(a)-(d)). Therefore, the performance curves obtained by the slope correction filtering (shown by star markers in Figure 6-15) are actually lower bounds of the performance obtained by the method when real frequency estimators are used. 6.3 Conclusions IF Estimation Results: In Section 6.1, simulations showed that the estimation error achieved by the proposed method is lower than that achieved by non-overlapping window-based method for coherence values higher than 0.6. This performance improvement is achieved because of 95 i) the use of highly localized filters that adapt their spatial aperture to local frequency contents of the signal, and ii) the IF is estimated on each pixel of the image. These results are based on the model used to generate the synthetic interferograms (2), (3) and the phase spatial distribution of (1). Although, they cover different ranges of coherence and spatial variation, the use of models closer to real topographic scenes is suggested, especially when dealing with abrupt topographic changes. The last implies the design of a procedure to effectively detect the non-interferometric areas like the areas with shadows, layover and phase aliasing. Phase Noise Filtering Results: Based on the performance obtained in Section 6.2 (Figure 6-15), it can be concluded that for high spatial phase variations and moderately high coherence (>0.5), the proposed filtering method achieves a better performance than slope-correction filtering because non-linear phase modeling follows the phase changes better than linear modeling. In addition, the advantages of the proposed method against standard multi-look averaging are substantial. The examples showed that standard multi-look averaging highly distorts the input phase signal, which results in large height estimation errors in comparison with the proposed method. Another advantage of non-linear phase modeling is that the averaging window can be extended in order to reduce more noise variance. Filtering with larger windows does not distort the input phase pattern as much as standard multi-look averaging or slope-correction filtering. However, the errors in the estimated IF are integrated (summed) when the non-linear phase model is being (equation (5) in Section D.3), which limits the size of the averaging window in the proposed method. An extensive evaluation of the proposed filtering method needs to be addressed. Such evaluation must use 2D phase unwrapping instead of ID phase unwrapping, because the second dimension could resolve or add ambiguities in the phase. 96 a) Input b) Non-linear phase modeling 50 100 150 200 range samples -> c) Linear phase modeling 50 100 150 200 range samples -> 50 100 150 200 range samples -> d) Standard averaging 50 100 150 200 range samples -> Figure 6-10: Phase distortion by different filters using 11 pixel window. a) Input phase interferogram - coherence=1.0. b) Filtered interferogram using the non- linear phase modeling filter (proposed method). c) Filtered interferogram using the slope-correction filter (linear phase modeling). d) Filtered interferogram using standard multi-look averaging. 97 a) Unwrapped phase Figure 6-11: Unwrapped Filtered Phase -1.0 coherence. a) ID Unwrapped phase through the center of images in Figure 6-10 along the range direction. b) Difference error between unwrapped filtered phase and original phase. - Continuous line: using the non- linear phase modeling filter. - Dash-dot line: using slope-correction filtering. - Dashed line: using standard multi-look averaging. 98 a) Theoretical IF b) Hypothetical Estimator Range samples -> Range samples -> c) IFs on noise-free data d) IFs on noisy data Range samples -> Range samples -> Figure 6-12: Range Frequency estimates. a) Theoretical Frequency of interferogram shown in Figure 6-10(a). b) Estimated Frequency by an ideal estimator, which takes 1 l x l 1 averages of the theoretical frequencies. c) Estimated Frequency by the proposed filter on the noise-free data of Figure 6-10(a). d) Estimated Frequency by the proposed filter on the noisy data of Figure 6-13(a). 99 a) Input 50 100 150 200 range samples -> c) Linear phase modeling 50 100 150 200 range samples -> b) Non-linear phase modeling range samples -> d) Standard averaging 50 100 150 200 range samples -> Figure 6-13: Phase noise filtering using 11 pixel windows. a) Input phase interferogram - 0.5 coherence. b) Filtered interferogram using the non-linear phase modeling filter. c) Filtered interferogram using the slope-correction filter (linear phase modeling). d) Filtered interferogram using standard multi-look averaging. 100 a) Unwrapped phase 20 40 60 80 100 120 140 160 180 200 220 range s a m p l e s -> b) Phase error v : I : 1 v \ > / ^ J \ \ / / / I 20 40 60 80 100 120 140 160 180 200 220 range s a m p l e s -> Figure 6-14: Unwrapped Filtered Phase - 0.5 coherence, 11 pixel window. a) ID Unwrapped phase through the center of images in Figure 6-13 along the range direction. b) Difference error between unwrapped filtered phase and original phase. - Continuous line: Filtered using the non-linear phase modeling filter. - Dash-dot line: Filtered using slope-correction filtering. - Dashed line: Filtered using standard multi-look averaging. 101 a) 5 pixels window 8- i V* B- ° -• a - - ^ - _ "* s s — _ 1 0.4 0.5 0.6 0.7 coherence -> 0.8 0.9 b) 11 pixels window o- -© • • • • o B 0.3 0.4 0.5 0.6 0.7 coherence ->• 0.8 0.9 Figure 6-15: M S E of the phase difference error for 5 and 11 pixel windows - Circ le markers: Non- linear phase modeling filter. - Star markers: Slope-correction filtering (linear phase modeling). - Square markers: Standard multi-look averaging. 102 Chapter 7 Results on Real Interferograms In this chapter, results with real SAR data are presented. First, the IP estimation method proposed in Chapters Four and Five is applied to real interferograms. Then, phase noise filtering is performed based on the estimated frequencies using the procedure presented in Section 6.2.3. In Section 7.1 results on the ERS Tandem Mission data are presented. The selected SLC (single look complex) images are taken with one-day difference. Consequently, the resulting interferogram is of high coherence (-0.7). Results on Radarsat data are presented in Section 7.2. Usually Radarsat data are of low coherence due to its longer repeating time period (24 days), but in this case a scene from a region with scarce vegetation was chosen which has moderately high coherence (-0.6). The interesting characteristic of this interferogram is that the IP spatial variations are very high, being a good example to show how the phase filtering based on IF estimates behaves versus standard methods. Results from both data sets are analyzed and the limitation of noise phase filtering based on the estimated IPs is discussed. 7.1 Results on ERS Tandem Mission Data Figure 7-1 illustrates the interferogram magnitude of a topographically rough area in British Columbia (Canada). The interferogram was obtained from two SLC (Single Look Complex) images acquired by the satellites ERS 1 and ERS2. The acquisition and processing parameters are listed in Tables 7.1 and 7.2. The interferogram has been half-band filtered to reduce layover 103 areas as explained in Section 2.2. Sub-sampling by 4 in the azimuth direction has been performed in order to have pixels of equal size in range and azimuth. Figure 7-2 illustrates the interferogram phase after flattening. The areas where the fringes are closer correspond to areas with steeper topographic slope, like the areas near the river in the lower right part of the image. The average magnitude spectrums in the range and azimuth directions are illustrated in Figures 7-3(a) and (b) respectively. Notice that the spectrum in the range direction is asymmetrical due to the application of half-band filtering before flattening. The spectrum tends to be concentrated in the lower frequencies. For example, outside the range frequency interval [-0.5, +1.5] rads/sample, energy is more than 10 dBs below the peak. This motivates the use of narrow-band Gabor filters concentrated at lower frequencies. In Section 4.4 the procedure suggested for this kind of interferograms was to perform range-modulation in order to match the interferogram's spectrum location and the filter's distribution of Figure 4.12. Given that the signal's spectrum is highly concentrated at lower frequencies, a simplification of the process is achieved by directly using the interferogram without range-modulation. The frequency response of the Gabor filters of Figure 4.12 (whose parameters are listed in Table 4.2) is duplicated in Figure 7-4. Figure 7-5(a) illustrates the real part of the filters unit pulse response in one-dimension. The frequency response of the filters illustrated in Figure 7-5(a) is illustrated in Figure 7-5(b), notice that the filters frequency responses intersect at -3 dB of its maximum magnitude. Using these filters, any frequency in the 2D circle with ratio n centered at the origin of the frequency plane will have a filter's response higher than -3 dBs. After convolving the input data with all filters, the IF is estimated on the channel with maximum magnitude as it was illustrated in the block diagram of Figure 4-3. The DESA-1 algorithm, given by equation (19) in Chapter Five, is used to estimate the IF. The resulting IF in the range direction is illustrated in Figure 7-6, and the IF in the azimuth direction in Figure 7-7. These two figures show some characteristics of the topography since the IF is geometrically related to the topographic slope in either the range or the azimuth direction. 104 Some areas in Figures 7-6 and 7-7 show abrupt changes between negative frequencies (black shades) and positive frequencies (white shades). In such areas the IF may not be properly estimated. These areas correspond to areas where the coherence is low. They are highly correlated with the bright areas of the interferogram magnitude (Figure 7-1) and correspond to steep slopes. In order to check the validity of the EF estimation, a very detailed topographic slope map is required. In this example, the pixel spacing is 12.5x12.5 meters; therefore, in order to check the estimated IPs at each pixel a topographic slope map of that resolution is required. Because such high-accuracy maps are usually not available, an alternative is to compare the resulting phase (height) after frequency integration between two points against in situ elevation measurements like the ones obtained with a GPS. This method implies that all the parameters in the height to phase transformation are already solved. Based on the previous IF estimates the phase is filtered using the non-linear phase modeling filter presented in Section 6.2. Figure 7-8 illustrates the results when seven pixel windows are used. Notice the reduction in noise compared to Figure 7-2. The standard multi-look average filter is used to obtain the results shown in Figure 7-9. For most areas the results in Figures 7-8 and 7-9 are similar, but many details are missing in Figure 7-9. Notice for example the loss of fringes in the lower right part of the image near the river. Process Method Parameter Registration Maximum Fringe Frequency -Range Half-Band filtering Band-Pass Filtering [0,71,] Flattening Demodulation Peak range freq. Down-sampling in azimuth Low-Pass filtering and decimation By 4 Coherence Mean whole scene -0.7 Table 7.1 ERS Interferogram Processing parameters 105 Parameter Value Units Incidence angle: near swath 20 Degrees Incidence angle: far swath 26 Degrees Near Slant range 833 Km Far Slant range 873 Km Swath Width: ground range 100 Km Swath Width: slant range 40 Km Altitude 782 Km Radar wavelength 0.0566 Meters Interferogram Baseline 42 Meters Acquisition date SLC1 Sept. 10/95 -Acquisition date SLC2 Sept. 11/95 -Scene's center latitude 51°35'20"N -Scene's center longitude 121°59'20"W -Pixel Spacing 12.5x12.5 Meters Table 7.2 ERS Interferogram Acquisition parameters 106 Input Interferogram - Magnitude A 13 E N < <— Range direction Figure 7-1: Magnitude of ERS Tandem mission interferogram This interferogram is used to illustrate the IF estimation method and phase filtering. The area corresponds to a 5 Km x 5 Km region in the upper valley of the Fraser River in B.C. (Canada). 107 Input Interferogram <— R a n g e direction Figure 7-2: Phase of ERS Tandem mission interferogram The areas where the fringes are closer correspond to steeper areas l ike near the river in the lower right part of the image. Note: In all interferogram phase plots the wrapped phase [-71,71] is normalised to [0,1] as shown by the grayscale. a) -3 -2 -1 0 1 2 3 Range spatial frequency —> rads/sample b) - 3 - 2 - 1 0 1 2 3 Azimuth spatial frequency —> rads/sample Figure 7-3: Magnitude Frequency Spectrum of ERS Tandem mission interferogram. a) Average Magnitude Spectrum in the range direction b) Average Magnitude Spectrum in the azimuth direction The Gabor filters used in the multiband filtering stage are superimposed in dashed lines. 109 GABOR FILTER'S CENTER FREQUENCIES c o 5o o o c 3 o o o o o o o o o o o o O o o o o O 0 O o o o o o © o O : o o o o o o o o ° o ° o o o 0 -1.5 -0 .5 0 X spatial frequency -0.5 -> rads/sample 1.5 Figure 7-4: Gabor Filters Frequency centers in 2D. 110 a) Filters unit pulse response Figure 7-5: Gabor Filters in ID. a) 4 Gabor filters unit impulse response hj(n). In dashed line its magnitude \h,{n)\ and in continuous line its real part. b) Frequency response of the 4 filters shown in a). I l l The simulations results in Chapter Six have already shown that standard multi-look average filtering distorts the input phase pattern more in comparison with slope-correction methods, in this example, such distortion is present in the areas of high phase variation. In conclusion, for this data set the EF estimation tracks the continuous IF changes of the input phase except in areas of low coherence. Filtering the interferogram with the non-linear phase modeling method presented in Section 6.2, reduces the noise variance without significantly reducing the spatial resolution in areas of close fringes. 7.1.1 Limitation of Noise Phase Filtering Method in the Presence of Residues. When the non-linear phase model is being built, the path illustrated in Figure D-5 (Appendix D) is followed to integrate the frequency. If its corresponding phase field is locally conservative the integrated phase does not depend on the chosen path and the phase model is unique. One sufficient condition for path-independence is that the integral around every single closed path is zero. When this integral is not zero there is a residue (see Appendix A). Residues reveal the existence of phase discontinuities. They are caused by noise or by non-interferometric sources like shadows, layover or real discontinuities in the topography (like in cliffs), which Insar is unable to map. After multi-band filtering, residues caused by noise will likely disappear if the coherence is high enough, but residues due to the topography will remain. The frequency estimates used to build the phase model and filter the interferogram, will fail in the vicinity of these non-interferometric occurring residues. This constitutes the main limitation of phase filtering based on a non-linear phase model built from the estimated EFs. Figure 7-10 and 7-11 show the residues of the raw interferogram (Figure 7-2) and the filtered one (Figure 7-8). These maps show in darker shades the areas where residues with negative charge exist {-27t), and in lighter shades the positive charge residues (+27t). The number of residues have lowered from 2.31% to 0.18%, as a result of filtering the interferogram. Filtering has eliminated the phase residues that were caused by noise. The areas where phase residues persist in the filtered interferogram are located mostly where high IFs were estimated and correspond to shadows or steep slope areas. 112 IF in the Range direction <— Range direction Figure 7-6: IF in the range direction - Estimated on the interferogram of Figure 7-2. Note: Frequencies higher than 0.15 rads/sample are shown in white, and frequencies lower than -0.15 rads/sample are shown in black. 113 IF in the Azimuth direction A I I C o T3 E N < 0.15 Ho.1 H0.05 -0.05 -0.1 <— Range direction Figure 7-7: IF in the azimuth direction estimated on the interferogram of Figure 7-2. Note: Frequencies higher than 0.15 rads/sample are shown in white, and frequencies lower than -0.15 rads/sample are shown in black. 114 Filtered Interferogram <— Range direction Figure 7-8: Filtered interferogram using the IF estimates of Figure 7-5 and 7-6 Note: The averaging window used is 7x7 pixels. 115 Filtered Interferogram <— Range direction Figure 7-9: Filtered interferogram using the standard multilook averaging filter Note: The averaging window used is 7x7 pixels. 116 a) Input Interferogram Figure 7-10: Residues of input interferogram (Figure 7-2). 117 b) Fi l tered Interferogram <— R a n g e direction Figure 7-11: Residues of filtered interferogram (Figure 7-8). 118 7.2 Results on Radarsat data Figure 7-12 illustrates the magnitude of the interferogram obtained from two SLC images acquired by the satellite Radarsat in area of the Death Valley in California (U.S.A.). Figure 7-13 illustrates the interferogram phase. Notice that the fringes are very dense, which shows that this scene is a good example to test the spatial localization of the IPs estimation and phase filtering. The acquisition and processing parameters are listed in Tables 7.3 and 7.4. The average magnitude spectrum in the range and azimuth directions are illustrated in Figures 7-14(a) and (b) respectively. The effective spectrums is highly concentrated at the lower frequencies, because the interferogram was flattened using the topographic information of a pre-existing D E M . The Gabor filters distribution of Figure 7-4 is used in this case as proposed in Section 4.4. In Figure 7-14, the filters' frequency response is superimposed in dashed lines. The estimated IPs are illustrated in Figure 7-15 and 7-16 for the range and azimuth directions respectively. The estimated IPs do not seem to continuously follow the topographic changes in steeper areas, they present patches rather than continuous values. These areas are highly correlated with the strong back-scattered areas of interferogram magnitude (Figure 7-12) and show high phase variations in the interferogram phase. In these areas the IP estimation is governed by the central frequency of the dominant filter output in the multiband filtering stage. Figure 7-17 illustrates the filtered phase using seven pixel windows in the non-linear phase modeling method presented in Section 6.2 where the IF estimates are used to model the phase. In contrast, Figure 7-18 shows the result of using the standard multilook averaging filter with seven pixel windows. The difference between both results is substantial, most of the topographic details seen in Figure 7-17 are missing in Figure 7-18. In conclusion, for this data set the IF estimation does not track the continuous IF changes of the input phase, but the central frequency of the dominant filter in the multiband filtering stage approximately samples the topographic slope. In this way, a fairly good filtering is performed in comparison with standard filtering. Finally, Figures 7-19 and 7-20 illustrate the residue maps of the input and filtered interferograms. The number of residues have lowered from 2.39% to 0.76%, as a result of filtering the interferogram. Filtering has eliminated the phase residues that were caused by noise. However, many residues still remain due to the topography complexity. 119 Parameter Value Units Incidence angle: near swath 39.6 Degrees Incidence angle: far swath 41.8 Degrees Near Slant range -1036 Km Far Slant range -1070 Km Swath Width: ground range -53.5 Km Swath Width: slant range -34.8 Km Altitude 782 Km Radar wavelength 0.056 Meters Interferogram Baseline -360 Meters Acquisition date SLC1 Sept. 24 /96 -Acquisition date SLC2 Oct. 18/96 -Scene's center latitude 36°29 N -Scene's center longitude 116°43'W -Pixel Spacing 6.0(r) x 8.9(a) Meters Table 7.3 Radarsat Acquisition parameters 120 Process Method Parameter Registration Maximum Correlation -Range Half-Band filtering Band-Pass filtering [-71,71] Flattening Using coarse D.E.M. -Down-sampling in azimuth Low-Pass filtering and decimation By 2 Coherence Entire scene mode Value -0.6 Table 7.4 Radarsat Processing parameters 121 Input Interferogram - Magnitude A I I C o E N < <— Range direction Figure 7-12: Magnitude of Radarsat interferogram This interferogram is used to illustrate the IF estimation method and phase filtering. The area corresponds to a 3.0 Km in range by 4.45 Km in azimuth region of the Death Val ley, Cal i fornia (U.S.A.) 122 Figure 7-13: Phase of Radarsat interferogram The areas where the fringes are dense correspond to very steep areas. 1 2 3 a) - 3 - 2 - 1 0 1 2 3 A z i m u t h S p a t i a l F r e q u e n c y — > r a d s / s a m p l e Figure 7-14: Magnitude Frequency Spectrum of Radarsat interferogram. a) Average Magnitude Spectrum in the range direction b) Average Magnitude Spectrum in the azimuth direction The Gabor filters used in the multiband filtering stage are superimposed in dashed lines. 124 IF in the Range direction Figure 7-15: IF in the range direction - Estimated on the interferogram of Figure 7-13. Note: Estimated IF grayscale is [-7t/4,7t/2] rads/sample. 125 IF in the Azimuth direction <— R a n g e direction Figure 7-16: IF in the azimuth direction - Estimated on the interferogram of Figure 7-13. Note: Estimated IF grayscale is [-7i/2,7i/2] rads/sample. 126 Filtered Interferogram <— R a n g e direction Figure 7-17: Filtered interferogram using the IF estimates of Figure 7-15 and 7-16 Note: The averaging window used is 7x7 pixels. 127 Filtered Interferogram A I I C q o "5 E N < 0.9 0.8 0.7 H0.6 0.5 0.4 0.3 0.2 0.1 <— R a n g e direction Figure 7-18: Filtered interferogram using the standard multilook averaging filter Note: The averaging window used is 7x7 pixels. 128 P h a s e Res idues <— R a n g e direction Figure 7-19: Residual phase map of input interferogram (Figure 7-13) P h a s e Res idues o 0 <— R a n g e direction Figure 7-20: Residual phase map of filtered interferogram (Figure 7-17) 130 Chapter 8 Conclusions 8.1 Summary The main difference between the proposed and existing methods is that existing methods try to estimate the local frequency of the signal by dividing it in windows. They assume that the signal is stationary within each window, which is not always true. The frequency estimators in these methods define a fixed window size that tries to balance the trade-off between spatial resolution and accuracy of the estimation. In contrast, in the proposed method the signal is adaptively filtered before the local frequency estimation is performed, allowing to have a simpler and more efficient frequency estimator that operates on fewer pixels. The proposed approach consists of two stages: multiband filtering and instantaneous frequency (IF) estimation. For the first stage, Gabor filters with variable space-frequency domain were chosen because of their minimum uncertainty characteristics. On the other hand, the DESA-1 algorithm was chosen for the instantaneous frequency estimation stage because of its low complexity and acceptable performance. The adequacy of the DESA-1 algorithm for instantaneous frequency was investigated using synthetic signals immersed in noise. Finally, the complete method was evaluated on synthetic and real interferograms with variable coherence and topography. 131 8.2 Conclusions Accurate frequency estimation is a difficult task when noise is present in the signal. At the same time, spatial localization is highly desired in order to follow abrupt topographic changes of the interferogram in rough areas. In order to combine these two requirements, the proposed method filters the data before estimating the frequency and estimates the frequency with a simple highly localized instantaneous frequency algorithm. In this way, the complexity of the frequency estimation algorithm is moved from the frequency estimation stage to the multi-band pre-filtering stage. The assumption made is that - in general - this results in better frequency estimates. Simulations showed that this assumption is valid when the coherence of the interferogram is moderately high (>0.6). The performance achieved by the proposed method is better than that achieved by non-overlapping window-based methods in that coherence interval. The comparison has been made against a hypothetical estimator that takes window averages on the noise-free signal frequencies. Therefore, the performance curves achieved by the hypothetical estimator are actually lower bounds of the performance curves achieved by real frequency estimators. The improvement achieved by the proposed method is due to: i) the use of highly localized filters that adapt their spatial aperture to the local frequency contents of the signal, and ii) the IF is effectively estimated using only five samples in each direction of the image (range and azimuth). The estimation of the local frequency in interferograms has relevance in three stages of the processing chain: coherence estimation, phase noise filtering and phase unwrapping. An application of the frequency estimation results was developed for the phase noise filtering stage. In this application, the phase was locally modeled using a non-linear model based on the estimated IPs. Then, after flattening the local topography the complex samples were averaged. Simulations showed that this alternative method for phase filtering, achieves a better performance than other two reported methods when the interferogram coherence is higher than 0.5. However, the method does not solve the problem of phase discontinuities in the presence of topographically induced residues. 132 8.3 Practical Application of the Proposed Method The proposed method shows advantages in interferograms with moderately high coherence (>0.6). These include ERS tandem mission interferograms with one or three days repeating time period, and Radarsat interferograms (26 days repeating time period) in regions with very scarce vegetation. However, many other interferograms have coherence less than 0.6 in large regions of the image. The possible use of the proposed method (using narrower bandwidth filters) for these interferograms, needs to be addressed. On the other hand, the interferogram data acquired by the S R T M mission constitutes an excellent scenario to test the method. In January of the year 2000, the S R T M mission acquired simultaneous pairs of SLC images using two antennas, eliminating the temporal decorrelation which is the main source of incoherence. A l l continental areas between latitudes 60° N and 60° S were mapped with high coherence [35]. In order to give an idea of the processing time in the proposed method, a 500x500 pixels interferogram is processed in less than 10 minutes using M A T L A B on a SUN SPARC 20 workstation. At the same time, unwrapping the phase of the same interferogram can last several hours depending on the phase variability. Therefore, the application of the proposed method to filter the phase is justified because it simplifies the phase unwrapping process by reducing the number of residues without excessively reducing the topographic information of the signal. Less accurate elevation maps can be obtained in less processing time by filtering the data with standard multilook averaging. The application of the proposed method is then oriented to the consecution of accurate elevation maps based on interferograms of moderately high coherence. On the other hand, the parameters in the phase-to-height transformation (such as the baseline) are assumed to be accurately estimated. 8.4 Future work The frequency estimation performance obtained were M S E error vs. coherence curves for different spatial distributions of the interferometric phase. The validity of the coherence model and the type of synthetic signals used were not questioned. This issue must be reviewed because it influences the decision of which algorithm is better. A more realistic model should be used for 133 the simulations, which includes non-interferometric features. In addition, a full analysis of the method's advantages must include testing the height accuracy after 2D phase unwrapping. Since the frequency content of the interferogram signal is spatially correlated, the implementation efficiency of the filters in the multi-band filtering stage can be improved by the use of Kalman filters [29]. Alternatives algorithms for IF estimation must be investigated to overcome the DESA-1 algorithm limitations found in Appendix C. Although phase noise filtering was not the main objective of this thesis, it constitutes an important application of local fringe frequency estimation. In areas where residues are present, phase modeling does not have an unique solution and phase filtering based on the estimated IFs fails. An alternative to this problem is to circumvent those areas by detecting the residues first and then growing a mask around them. If all residues are masked and each mask encloses an equal number of positive and negative residues, the remaining phase field is conservative and phase filtering without distortion is achieved. The consecution of accurate frequency estimates by the proposed method opens the possibility of obtaining the elevation map by integrating the frequency estimates rather than by unwrapping the phase, which is usually a difficult task. The hypothesis that the IF estimates fully represents the scene's topography is based on the complete representation of a monochromatic signal, that can be achieved by the mono-component A M - F M model [28]. The main advantage of this alternative is that the IF estimates are less noisy than the filtered phase, creating an alternative to unwrap noisy interferograms. The field of estimated frequencies is also non-conservative, therefore, algorithms similar to the used in phase unwrapping could be used. The above implies that further investigation is required in order to achieve better results for noisier interferograms (coherence less than 0.66). This can be achieved by reducing the bandwidth of the filters in the multi-band filtering stage in order to filter out more noise at the expenses of losing spatial resolution. The automatic selection of the filters parameters based on the input characteristics of the interferogram should also be addressed to eliminate user interaction. 134 In addition, the conditions under which the mono-component A M - F M signal model represents the interferogram signal need to be studied for each acquisition geometry (ERS, Radarsat, SRTM, etc.), in order to obtain quantitative smoothness conditions on the scene's topography. 135 Appendix A Interferogram Residues A . l Definition The closed contour integral of phase differences along any valid path in a properly unwrapped interferogram phase should result in a net value of zero residual phase. Zero residual phase from a closed path integral corresponds to returning to the same point in the scene's topography. However, discontinuities in the image phase yield anomalous paths for which the residual phase difference is not zero. As stated by the residue theorem [32], if there are inconsistencies in the phase the sum of unwrapped phase differences in a closed contour is not zero but a multiple of 2n. | Vlf/(x)-dx = 2TTX( sum of enclosed phase residues ) (1) where l//(x) is the phase function and x the vector coordinates. Residues reveal the existence of phase discontinuities. They are caused by decorrelation noise or by non-interferometric features like shadows, layover or real discontinuities in the topography (such in cliffs), which the interferogram pair is unable to map. 136 Positive residues charges (+1) are assigned for positive residual, while negative residues charges (-1) are assigned for negative residual. When residues charges are balanced in a region, the integral around that region is zero for any simple path chosen. Thus, consistent phase unwrapping is possible only if all integration paths do not encircle unbalanced residue charges. A.2 Detection Algorithm Residues can be located by adding the wrapped phase along the smallest closed path (2x2) in the discrete two-dimensional array [32]. If y/(m,n) is the phase of the interferogram, for each pixel m,n the following wrapped phase differences are computed: Ax(m,n) = W{y/(m,n + \)-y/(m,n)}, A 2 (ra,n) = W{y/(m + l,n + l)-if/(m,n + l)}, A3(m,n) = W{y/(m + \,ri)-y/(m + \,n + \)}, A4(m,n) = W{y/(m,n)-y/(m + l,n)} (2) where W{-\ denotes the phase wrapping operator. Then, the phase residue is: 4 q(m,n) = £ \(.m,n) (3) /=i As an example, Figure A-1(a) shows a phase function described by: m(f) = 6(7), -n<6<n (4) where 9(r) is the angle between the vector r = [x, y] and the horizontal axis. The center of the image corresponds to the origin of coordinates. The topography associated 137 to the phase (p(r) in (4) can be visualized as a conic surface with its vertex located at the origin. y Transforming (4) to x, y coordinates using (p(x, y) = arctan(—), the phase gradients x are given by: <V = -y d £ = x ~v 2 2 ' "v 2 2 ^ ' dx x + y dy x + y The phase is discontinuous at the origin and this causes the gradient to be undefined. Figures A-1(b) and (c) show the phase gradient computed by taking discrete differences of the wrapped phase in the vertical and horizontal directions. Near the origin the phase gradient estimates behave like the electrical field of a dipole. Figure A-1(d) illustrates the phase residues computed using equations (2) and (3). There is only one residue found at the origin and its value is +2n. The situation will be similar in non-interferometric occurring residues like in areas of shadows, layover or cliffs, but the residues will appear along the phase discontinuities. 138 a) Input phase b) Vertical Phase Gradient Figure A - l : Example illustrating naturally occurring residues a) : Phase function b) and c): Phase gradients on the vertical and horizontal directions, d): Phase residues 139 Appendix B Relation of the Proposed Method and the STFT The use of the Short Time Fourier Transform (STFT) in the estimation of local frequency in interferograms has been reported in [5, 15]. The STFT of a signal x(n) is defined by: The squared magnitude of the STFT is an estimate of the periodogram of the signal. The frequency of the peak in the periodogram constitutes a Maximum Likelihood (ML) estimate of the local frequency of the signal (see Chapter 3). The spatial extent of the window w(n) has an important effect upon the frequency resolution. If a short window is used to compute (1), the transform will exhibit good resolution in space, but the frequency resolution will be low. Alternatively, if a longer window is used in order to increase the frequency resolution, the spatial resolution (localization) will be poor. The alternative proposed in this work is clearly seen when (1) is expressed as a convolution of the frequency-shifted signal x(m)e~j0"n with w(m). At a specific sample, (1) can be re-SSTFT(n,w) = ^ x(a) w(n-a) e (1) a written as: SSTFT(n,co) = x(n)e-jmn*w(n) = {x{n)*w(n)eJ(on} e - jcon s STFT (n,CO) = {xin)*^ (n)} (2) 140 If the frequency is held constant CO- C00, then (2) is the magnitude response of a filter hm to an input x(n). The frequency response of hm is (co) = W(co -CO0), which is the same as the window function frequency response shifted to the frequency C0o. The analysis above shows that the proposed approach alleviates the problem of poor space-frequency resolution of the STFT, because the window is locally optimized. For instance, a high frequency event is analyzed with a small window and a low frequency event is analyzed with a large window. Thus, by varying the window size and approximately match filtering the signal, the proposed approach applies an appropriate window to the signal at all samples. The task of finding which window fits the local characteristics of the input best, is simplified to finding which filter hm(n) maximizes the output in (2). Such maximization could be performed over all frequencies CO in the frequency plane. Instead, a much smaller number of filters are used and the Instantaneous Frequency (IF) is estimated at the pixel level on the output x(n) * hw (n) using the DESA-1 algorithm presented in Chapter Five. 141 Appendix C Evaluation of the DESA-1 IF Estimation Algorithm The purpose of this appendix is to evaluate the adequacy of the DESA-1 IF estimation algorithm for the InSAR application. The proposed methodology can also be used to compare other IF estimators against DESA-1. However, such comparison is not done. C l Frequency Distortion of the DESA-1 Estimator There are several characteristics that need to be considered in order to evaluate the adequacy of an IF estimator. The most important ones for the InSAR application are: • How well the estimator works estimating the IF of a complex signal in the interval • How well the IF estimator works when the IF varies rapidly? • How well the'estimator wOrks estimating the IF in the presence of noise? • How the performance of the IF estimator is affected by the multi-band filtering stage? In order to evaluate the DESA-1 IF estimator a test signal is used for which several parameters control the excursion of the IF, the velocity of the IF variation and consequently its effective bandwidth. An adequate test signal is a F M signal, where the modulation parameters Q m o d , Af, Q c are changed to obtain different IF excursions and bandwidths. 142 The test F M signal is described by: x(n) = exp(j(/)(n)) = exp( j[Qcn + Af s i n ( Q m o d n ) ] ) (1) Its theoretical instantaneous frequency is: d(j)(ri) dn = Q,.(n) = Q c + A / Q m o d c o s ( Q i n o d n ) (2) and the absolute value of the normalized estimation error is defined as: £(n) = Q,.(w)-fl,.(w) .Q,.(n) x l 0 0 % (3) where is the estimated IF and Q,.(n) is the theoretical IF. First of all, the performance of the IF estimation in the range [0,7r] is studied for the noise-free case. The term: distortion of the IF estimation is used in reference to the error between the real and the estimated IF in this case. In the example of Chapter Five - Figure 5-1 - it was shown that the DESA-1 algorithm estimates the IF with a maximum error of 0.62% when the excursion of the IF is within [71 / 4,71 / 2] rads/sample and the frequency at which IF varies is £ 2 m o d =7z7100. However, if the excursion of the input signal IF is changed from [7114,71/2] to [0.05,7z74] rads/sample as in Figure C - l , the error increases dramatically, implying that the IF estimator does not work well at low frequencies. This behavior of the DESA-1 Algorithm at low frequencies can be fixed by modulating the input signal to a higher frequency. Because the band at which the signal belongs has been previously detected in the multi-band filtering stage, the low frequency channels can be modulated to a higher frequency COH , where the IF can be more accurately estimated and then COH be subtracted from the estimated IF to obtain the real one. Figure C-2 illustrates this procedure where the COH used was 0.46 rads/sample. The maximum error is reduced from 100% to 2.67%. 143 The procedure is generalized to obtain the best performance of the DESA-1 by splitting the [0,7l] range in three sections, and shifting the signal center frequency (by modulation) to a frequency closed to 7112 where the estimator works better. In this way, an approximately constant performance is obtained for all frequencies. C.2 Distortion of DESA-1 with Fast IF Changes When the EF of the input interferogram changes fast there will be more error in the approximation (17) of Chapter 5. This can be noticed in the example of Figure C-2 where the error is higher for the areas at low frequencies where the derivative of the IF is high, i.e. around samples number 70 and 130. To evaluate the Distortion of DESA-1 with fast changes of IF, the F M modulating frequency ^ m o d * s v a r i e d in the F M test signal of (1), which results in faster or slower changes of the input IF. The RMS error of the estimation is defined as: ^^IW^^)xm% (4> where Q,(n) is the estimated IF, £2,(n) is the theoretical IF derived from (2), and N is the number of samples of the test signal. This kind of error is used because it produces a percentage figure that is illustrative of the accuracy of the estimator. Figure C-3 illustrates the RMS value of the error as defined in (4) for £ 2 m o d varying between 0 and 3 rads/sample. Q c = 3^"/8 and Q, 6 [ ^ /4 , 7T / 2 ] rads/sample. In the absence of noise (diamond markers) the error increases up to = 28% when the F M modulating frequency Q m o d reaches = 2.4 rads/sample. The IF of an interferogram phase signal varies accordingly to the changes in the topographic slope, but the rapidity of these changes is not known in advance. For the example above the frequency £ 2 m o d at which significant phase aliasing occurs is 2Q,modA = 71 => QmodA =7112. 144 rads/sample. Notice that the estimation error is moderately high (= 23%) before phase aliasing occurs, implying that the algorithm is less able to track very fast changes of the IF even if the signal is properly sampled. This effect can be fixed by over-sampling the signal, which lowers the signal's frequency contents but increases the computational load of the algorithm. To analyze the case when noise is present in the input signal, a simulation has been done assuming the noise is additive and normally distributed with zero mean. Figure C-3 illustrates the case when the variance of the noise is 0.015 (star markers) and 0.03 (circle markers). At low F M modulating frequencies the difference in error between the signal with and without noise is substantial, but at modulating frequencies higher than 0.5 rads/sample the different is small. Therefore, for high modulating frequencies most of the error is due to the incapacity of the algorithm to follow fast changes of the IF, and the behavior of the curves are close to the case of no noise. C.3 Noise Reduction by Pre-filtering In this section the efficacy of pre-filtering the signal to improve the performance of the IF estimator is studied through simulation experiments. The results of these experiments are RMS error vs. SNR curves as indicatives of the performance of the proposed method. Experiment 1: Consider the test signal of (1) with £ 2 m o d = 71 /100 and other parameters as in last section. For the no-filtering case the RMS estimation error - defined by (4) - is illustrated by the dashed line in Figure C-4. By filtering the signal with a Gabor filter centered at the signal's carrier frequency, the IF RMS estimation error is lowered from 34% to 16%. In Figure C-4 the continuous line illustrates the case when filtering is applied with a Gabor filter of a-3, and the dash and dot line when filtering is applied when a-1.5 (Remember that higher values of or imply a narrower filter). Experiment 2: As seen in last experiment if the filter's bandwidth is narrower more noise is filtered out, but the filter's bandwidth can not be excessively narrow because the signal's IF contents can be altered. This is illustrated in Figure C-5 where the same experiment has been done to a test signal whose £ 2 m o d =71/16. The error obtained in this case, is higher for the filtered signal with a=3 (continuous line) than for the non-filtered signal (dashed line) for SNR higher than approximately 11 dB. The above example illustrates that the distortion error caused 145 by the DESA-1 IF estimator can be higher than the effect of noise when the local bandwidth of the signal is large compared with the bandwidth of the filter. To better illustrate this distortion error, consider Figure C-6(a) where the Magnitude of the Spectrum of the F M test signal used to obtain Figure C-4 is shown. If the input signal is convolved with a Gabor filter centered at 2>n IA rads/sample, there will be a distortion in the estimated IF because the response of the filter is not flat in the pass-band. In Figure C-6(a) the frequency response of the Gabor filter with o=1.0 has been superimposed to the spectrum of the input signal, and in Figure C-4(b) the spectrum of the filtered signal is illustrated. Notice that the distortion is very low. If the bandwidth of the Gabor filter is narrower there will be more attenuation in the frequencies at the edges of the pass-band. Figures C-6(c) and (d) illustrate the case for o=3.0. Notice in this case that the filter reduces stronger components of the spectrum. If £ 2 m o d is increased, the bandwidth of the signal increases and there will be more attenuation due to the filter. Figures C-7(a)-(d) illustrate the case when £ 2 m o d = 711A . Notice that for a Gabor filter with <j=3.0 - Figure C-7(d) - , only the carrier frequency component remains after filtering, which highly increases the distortion error. The general case is illustrated in Figure C-8. The DESA-1 algorithm has been applied to the test signal (1) in the absence of noise after filtering it with a Gabor filter whose cr varies from 0 to 6 in steps of 0.5. The different values of £ 2 m o d (labeled as Q in Figure C-8) are expressed in rads./sample. The excursion of the IF is kept constant within \7T I A,7l 12\. For each curve, the RMS estimation error is computed for different values'of a. For o=0 (no filtering) the error achieves the minimum since no noise is added to the signal. For low values of cr the RMS error remains approximately constant until a value of <y=acut. At this value the effective band-width of the signal (BWS) is higher than the filter's band-width (BWg). The values of a=acM at which BWs=BWg, vary from curve to curve since BWS is a function of £ 2 m o d . 146 a) Real part of Input F M signal 50 100 150 200 250 300 350 400 b) Estimated Instantaneous Frequencies - in rads/sample 50 100 150 200 250 300 350 400 Sample number —> Figure C - l : DESA-1 IF estimation for IF within [n/3,2n/3] a) Real part of the signal: x(n) = exp(j[£2cn + Af sin(£2modn)]) with parameters Q m o d = 7z7100,Qc =0.42,A7 =23.4 b) Estimated frequency using D E S A - 1 , Q. s [0.05,TZ74] c) Normalised error. The error is very high at low frequencies reaching up to 100%. 147 a) Real part of Input F M signal 50 100 150 200 250 300 b) Estimated Instantaneous Frequencies - in rads/sample 350 400 50 100 150 200 250 c) Normalized error 150 200 250 Sample number —> 300 350 400 400 Figure C-2: DESA-1 IF estimation with modulation a) Real part of the signal: x(n) = exp(j[Q.cn + Af sin(Qm o d«)]) with parameters Q m o d = n 1100, Qc = 0.42, Af = 23.4 • b) Estimated frequency using DESA-1 and modulation at 0.46 rads/sample c) Normalised error - the maximum error is reduced from 100% to 2.67% by modulating the signal to a higher frequency where the estimator behaves better. 148 DESA-1 Performance vs. Modulating Frequency 35 30 25 20 M 15 10 0 * \ noise var.: \ \ noise var.=0.015 A 0.03 •. / / / noise var=0 >>1. -_ o- • / — » . • -o-• <> 0.5 1 1.5 2 Modulating Frequency in rads/sample —> 2.5 Figure C-3: DESA-1 R M S estimation error when IF varies fast The test signal is: x(n) = e\p(j[Cicn + Af sin(Qm o d/i)]) with parameters: Q.c = 3 ^ / 8 and variable Q m o d , A r The modulating F M frequency Qc is varied from 0 to ;r rads/sample, Notice the error increase at faster variations of the IF. 149 Figure C-4: DESA-1 Performance for Q m o d = n/100 The test signal is: *(#i) = exp(./[r2cn + Af sin(Qmodn)]) with parameters: Qc=3TT/& and variable Q m o d =^-/i00. The SNR is ioiog(-V)-150 Figure C-5: DESA-1 Performance for Q m o d = a/16 The test signal is: x(n) = exp(j[Qcn+ Af sin(Qmoin)]) with parameters: Qc=3x/8 and variable Q.moi -nl\6-The SNR is ioiog(-^)-151 For values of o> aCMi strong components of the signal spectrum are reduced and the RMS error increases asymptotically for values of £ 2 m o d < TT /16, and reaches a constant for values of £ 2 m o d >n 116. The reason for this behavior is explained in the sequel. For Q m o d > ^ / 1 6 the signal spectrum mainly consists of two components: the carrier frequency Qc and two peaks at Q c ± £ 2 m o d - as illustrated in the example of Figure C-l . After strong filtering only the carrier frequency remains as illustrated in Figure C-7(d), and further filtering does not affect the estimation error. On the other hand, for £ 2 m o d <Til 16 the signal spectra consists of components densely distributed within [£2C. + Q m j n , £lc + £ 2 m a x ] - as illustrated in Figure C-6(c) - , and after filtering the frequency components are gradually reduced when the values of cr increases. The distortion error in Figure C-8 means that rapid variations of the input signal's IF are smoothed producing a narrower band signal at the output of the multiband filter. The smoother estimated IFs are closer to center frequency of the filter as illustrated in Figure C-9(b) for the filtered signal of Figure C-9(a). The spectrum of this signal was illustrated in Figure C-7(b) after filtering it with a Gabor filter of cr=1.0. Figure C-10 illustrates the case when the signal was filtered with a Gabor filter of cr=3.0. Its spectrum is illustrated in C-7(d). For this case the estimated IF becomes close to the carrier frequency of the signal Qc, since other components of the signal are filtered out. C.4 DESA-1 Performance After Multi-band Filtering As seen in last section pre-filtering the signal with Gabor filters greatly reduces the error caused by noise in the IF estimation. The filter must be centered close to the local center frequency of the signal, and its bandwidth not be narrower than the local signal bandwidth. In reality, the local center frequency of an interferogram signal is unknown a-priori so the use of the multi-band filtering stage before IF estimation is well justified. The application of the proposed multi-band filtering scheme before IF estimation improves the performance of the IF estimation, as illustrated in Figure C-4 by the line with square markers for the experiment 1 of last section. The performance is better at all SNR, but its improvement is more significant at low SNR. The filters used are the ones illustrated in Figure 4-6 of Chapter Four 152 a) Input and Filter with sigma=1.0 c) Input and Filter with sigma=3.0 Figure C-6: Effect of the Gabor filtering for Q m o d = n 1100. a) Spectrum of the input signal (dotted line) and Gabor filter frequency response for sigma=1.0 (continuous line). b) Spectrum of the filtered signal. c) Spectrum of the input signal (dotted line) and Gabor filter frequency response for sigma=3.0 (continuous line). d) Spectrum of the filtered signal. Notice that strong signal components are removed. 153 Figure C-7: Effect of the Gabor filtering for Qmoi = nl 4 . a) Spectrum of the input signal (dotted line) and Gabor filter response for sigma=1.0 continuous line. b) Spectrum of the filtered signal. c) Spectrum of the input signal (dotted line) and Gabor filter response for sigma=3.0 continuous line. d) Spectrum of the filtered signal. Notice that only the carrier frequency component remains. 154 It is important to highlight that using one Gabor filter with a narrower bandwidth will not improve the performance, since it will filter out signal components increasing the error at high SNR. To illustrate this point consider Figure C-5 where the experiment 2 of last section has been done on a signal with faster IF variations ( £ 2 m o d = n 116 rads/sample). The continuous line curve with cross markers corresponds to the application of a stationary Gabor filter with o=3 whose performance is poorer than the case where no filtering is applied (dash line) at high SNR. In the experiment of Figure C-5 the use of the multi-band filtering is illustrated by the continuous line curve with square markers. In this case the performance is better than the other curves at low SNR and only slightly poorer than the no-filtering case at high SNR. The reason why in Figure C-5 the performance is better for the no-filtering case at high SNR, is that the signal's bandwidth is higher than the filter's bandwidth. Therefore, there is a distortion error produced by the filter, which is higher than the distortion caused by the noise. C.5 Summary and Conclusions Because its low computational complexity, the DESA-1 algorithm was chosen for the purpose of IF estimation in interferogram signals in the method proposed in Chapter Four. The computational complexity is an important criteria due to the high amount of data to be processed. The performance of the DESA-1 IF estimator was obtained by using a specific test signal corrupted with additive noise. The results showed that the RMS estimation error is considerably low when the signal's IF changes slowly, and it is moderately high (23%) when the IF changes very rapidly. The performance results obtained in this chapter - namely Figures C-3 to C-5 and C-8 - show that the DESA-1 algorithm performance is acceptable for IF estimation in interferogram signals in the proposed scheme. One disadvantage of the DESA-1 algorithm is its incapacity to follow very fast changes of the IF. This can be solved by over-sampling the input data, but it also increases the computational complexity of the algorithm. The proposed method is not restricted to the use of the DESA-1 estimator. In general, any IF estimator will benefit by receiving a less noisy input signal. Therefore there is room to search for another IF estimator and be plugged into the proposed scheme. 155 30 r 25 A ' 20f . o L_ CD | 151 re I *• 05 ^ LU R M S Estimation Error for different Gabor Filter Apertures -o-_ _ _Q O ©- 0 j — 9 & 0 Q=K/4 / 0 _ - 0 C/) 10h 0 0 0 Q=JI/8 0 0 6- - - 0 0-4) 0 / >v. . 0 " / / / 0 P _ -o-- 0 ' _ _ e _ _ © - - 0 ' - -0-£2=71/32 0 - -or n=7t/64 0 0 " Q=7t/10C^) - - © " 2 3 4 Gabor filter s igma —> Figure C-8: Influence of the Gabor filtering on the IF Estimation error in the absence of noise. Each curve illustrates the RMS estimation error vs. the aperture of the Gabor filter (sigma) for different F M modulating frequencies ^ m o d . The FM modulating frequency controls how fast the IF of the signal changes. 156 a) Filtered Signal when sigma=1.0 j i i i i i i i i l 5 10 15 20 25 30 35 40 45 50 c) IF Estimation error j 100 1 1 i i i i i i i i i 5 10 15 20 25 30 35 40 45 50 sample number —> Figure C-9: Effect of Gabor filtering for o=l, Q m o d =n/4-a) Filtered signal. b) Estimated and theoretical EF in rads/sample. c) Absolute value of the Estimation error in percentage. 157 a) Filtered Signal when sigma=1.0 20 25 30 sample number —> Figure C-10: Effect of Gabor filtering for o=2>, Q, = nl 4 • a) Filtered signal. b) Estimated and theoretical IF in rads/sample. c) Absolute value of the Estimation error in percentage 158 Appendix D Phase Noise Filtering Methods This appendix presents the phase noise filtering methods used in the performance evaluation of Section 6.2.1 (Chapter 6). The first two methods constitute prior work and are briefly presented in Sections D.1 and D.2. The third method, on the other hand, is a new method proposed in this thesis, which uses the output of the IF estimation method -presented in Chapter 4 - to model the phase. A detailed description of the new method is done in Section D.3. D.1 Standard Multi-Look Averaging In this method, the phase noise variance is reduced by combining returns from several interferometric pairs of equal mean phase [7]. In the absence of more than one interferometric pair, adjacent pixels can be averaged if they correspond to the same phase (height) - see equation (9) in Section 2.4.1. The filtered phase is obtained by averaging the complex interferogram values of adjacent pixels. 0ave(i,j) = arg{ Y,c(i + mJ + n ) e m ^ n ) } (1) m,n where c(i, j)ej^',j) is the input interferogram, and </>ave(i,j) is the filtered phase. 159 Standard multi-look filtering reduces the noise variance at the expense of losing spatial resolution, which lowers the contrast of high frequency fringes. D.2 Slope-Correction Filter The slope-correction filter proposed in [6] models the phase linearly to correct the effect of the topographic slope, and then performs standard multi-look phase filtering to reduce the noise variance. For each pixel P in the input signal, the topographic linear phase model based on the local frequency estimations is: where Q)x, Wy are the range and azimuth local estimated frequencies of a window centred on pixel P, and m, n are the coordinates of the pixels in that window. Although the phase behavior is not necessarily linear, it can be assumed so if the instantaneous frequencies do not vary much within the window. To illustrate the topographic phase model construction, consider the set of 16x16 pixels in window A of Figure D-l(a). Figure D-l(a) shows a synthetic interferogram corresponding to a=50 in equation (1) of Chapter 6, and 0.8 coherence. The input phase for the noise-free case is drawn in Figure D-2(a) and its topographic model in Figure D-2(b). The input phase is wrapped within [-7T,7T] but the modeled one is not. Considering only the phase term, the product between the input interferogram window ej^m'n) and the complex conjugate of the interferogram generated by the linear phase model e j< t>"" ,{m-" ) is: (ptop(m,n) = OJxm + COyn (2) (3) 160 It has a phase difference 6 which it is close to zero as shown in Figure D-2(c). In Figures D-3(a)-(d), I D slices of Figures D-2(a)-(d) have been drawn. In Figure D -3(a) the continuous line with circle markers illustrates the input phase and the line with cross markers illustrates the linearly modeled phase. Figure D-3(b) draws the phase difference between the model and the input phase. Although noise has not been introduced to the signal, the phase difference 9 is not exactly zero because the signal phase is not strictly linear. In other words, the IFs QX,Q are not constant within the window as illustrated in Figures D-3(c) and (d). a) I nput 20 40 6 0 80 100 20 40 60 80 100 Range sample -> Range sample -> Figure D - l : Input and filtered Synthetic interferogram a=50, 0.8 coherence. Figures D-4(a)-(d) corresponds to the noisy case showed in Figure D - l (a). The difference phase 0 in Figure D-4(b) has the effects of decorrelations in the interferogram. Its variance is highly reduced when neighboring samples are averaged by standard multi-look filtering. Figure D-4(b) illustrates only a I D slice of the phase 6 for 1x16 pixels after slope correction, but all 16x16 samples in the window wi l l be averaged to estimate the phase of pixel P. 161 Phase Phase (rads) a> (rads) Figure D-2: Linear Phase Model a) Input for noise-free case of Figure D- l (a) b) Topographic phase linear model. c) Phase difference between input and model d) Magnitude of input IF = [Q.x,Q.y]. 162 b) CO 5 10 Range sample -> 15 5 10 Range sample -> A I tn •a co co jo 03 T3 CO -0.06 -0.08 -0.12 -0.14 5 10 Range sample -> 15 5 10 Range sample -> Figure D-3: ID slices on Figure D-2 a) Input phase drawn with circle markers. Modeled phase with cross markers b) Phase difference between input and model c) Range direction IF Q,x .Theoretical in dash-dot line, estimated in continuous line d) Azimuth direction IF Q.y .Theoretical in dash-dot line, estimated in continuous line. 163 -0.14L 5 10 Range sample -> 15 -0.141 5 10 Range sample -> 15 Figure D-4: Noisy case a) Input phase drawn with circle markers. Modeled phase with cross markers b) Phase difference between input and model c) Range direction IF £lx .Theoretical in dash-dot line, estimated in continuous line d) Azimuth direction IF Q.y .Theoretical in dash-dot l ine, estimated in continuous line. 164 After the linear phase model is obtained - (/> in (2), standard multi-look filtering is performed on the locally corrected interferogram. The filtered phase is: where c(i, j)ej^''j) is the input interferogram, and cpfilt(i, j) is the phase of the filtered one. Figure D-l(b) illustrates the filtered interferogram of Figure D-l(a) using linear phase modeling and 7x7 pixels averaging window in (4). In this example, the spatial variations of the instantaneous frequency were low. Thus, the linear phase modeling was appropriate since the noise contributions were much higher than the discrepancies between the model and the real phase. However, the linear model will not be valid when much larger IF spatial variations are present in the input interferogram. D.3 Non-Linear Phase Modeling Filter (Proposed Method) When the phase behavior is non-linear, the difference phase 0 in (3) will have not only the effects of noise but also real phase information that will be filtered after averaging. To overcome this situation the EFs are estimated on a pixel basis rather than on a window basis, using the scheme proposed in Chapter Four. Then, a more accurate model is built by integrating the EFs. For each pixel P in the input signal, the topographic non-linear phase model is: </>fdt(i,j) = a r g { £ c(i + m,j + n)e j<p(i+m,j+n) -j<P,0p(m<n') (4) m,n n y/top{l,n) = 2^(1 .7) for n = l,...,d. and 165 m i//top(m,n) = l//top(l,n) + ^ &x(i,n) form = 2,...,d. (5) i=2 where Q ; c (m,n),Q ) ,(m,n) are the range and azimuth estimated IFs and m,n are the coordinates of the d2 pixels in a window. This non-linear phase model is used to locally correct the effect of the topography in the phase. Then, a standard multi-look phase filtering is performed. Equation (5) is equivalent to performing a two-dimensional discrete integration through two one-dimensional sums. Figure D-5(a) illustrates the paths followed by the ID summations. Initially, the IPs are summed along the azimuth direction in the first row. Then, the IPs are summed in the range direction for all columns. In the implemented algorithm, the model is built starting from the window center rather than from the window corner. Figure D-5(b) illustrates the followed paths in this case. The path indicated by the continuous line is first followed and then the path indicated by the dashed lines. The value of y/top at the window center is initialized to zero. As seen in Section 4.2, the integration of accurate IP estimates presents an alternative to obtain the unwrapped phase. However, the interest here is only to locally model the phase accurately. There are three main differences, which make the problem simpler: • There is not need of adding integer multiples of 2n to adjust the phase discontinuities like in the phase unwrapping case. • The estimated IPs are smoother than the real phase differences. • The integration is only performed locally in small windows. As an example of the topographic non-linear model construction, consider the set of 1x21 pixels in Figure D-6(a). The input phase for the noise-free case is shown by the continuous line, the linearly modeled phase (j) (2) by circle markers, and the non-linear 166 modeled y/top (5) by star markers. Notice that the modeled phases are zero at the center of the window (sample 11 th). For illustration purposes, the theoretical IFs Q(n) are used instead of the estimated ones Q(n) to compute (5). These are drawn in Figure D-6(b) with star markers. The average theoretical EF in the window is used to build the linear model in (2), and it is computed as: CO = i£ (6) " n=l This average frequency CO is illustrated by the dashed line in Figure D-6(b). The x and y sub-indices have been omitted because the same operation is performed in both dimensions. n=1 m=2 m=5 n=5 n=-d/2 • T • • azimuth direction a) m=-d/2 m=d/2 n=d/2 A A A i ! | — i. i 1 j s i ! T • T T azimuth direction b) CD Figure D-5: Integration paths a) Summation paths starting from the upper-left corner for d=5 b) Summation paths starting from the window centre. In Figure D-6(c) the slope-corrected phase is drawn with circle markers for the linear model and with star markers for the non-linear model. The difference between the real 167 and modeled phases reaches a value of almost 2 radians for the linear case, but it is exactly zero for the non-linear case because there is no error in the phase model since the theoretical IFs were used. Figures D-7 illustrate the case when the estimated IPs Q(n) are used to model the phase y/top in (5), instead of the theoretical IFs. These estimated IPs are drawn in Figure D-7(b) with star markers. Notice that they are closer to the real frequency (continuous line) than the average theoretical IF CO , which is illustrated by the dashed line. The phase difference in Figure D-7(c) is not exactly zero for the non-linear case as in Figure D-6(c) because of the incapacity of the DESA-1 algorithm to exactly follow the input phase changes. However, the difference between the modeled and the input phases is much higher for the linear model (circle markers). After the non-linear phase model is obtained - y/t in (5), standard multi-look filtering is performed on the locally corrected interferogram. The filtered phase is: VjuiUJ) = arg{ J cii + mij + tieM^+'V""™ } (7) m,n where c(i, j)eJ^l'j) is the input interferogram and if/filt (i,j) is the phase of the filtered one. Application examples of the proposed filter (7) on synthetic and real interferograms are found in Sections 6.2, 7.1, and 7.2. 168 Figure D-6: Linear vs. non-linear phase modeling - Theoretical IFs. a) Input phase drawn with continuous line, linearly modeled phase drawn with circle markers, and non-linearly modeled phase with star markers b) Range direction IF Q,x . Theoretical in continuous line, window average with dashed line and estimated IF in dash-dot line with star markers. c) Modeled phases: linear model - circle markers, non-linear model - star markers. Input phase value at the window center - dashed line. 169 -0.2 1 10 15 samples -> 20 Figure D-7: Linear vs. non-linear phase modeling - Estimated IFs. a) Input phase drawn with continuous line, linearly modeled phase drawn with circle markers, and non-linearly modeled phase with star markers b) Range direction IF Q.x . Theoretical in continuous line, window average with dashed line and estimated IF in dash-dot line with star markers. c) Modeled phases: linear model - circle markers, non-linear model - star marker. Input phase value at the window center - dashed line. Notice that the non-linear phase model is closer to the real phase. 170 References [1] Q. Lin, J. Veseky, and H . Zebker, "Topography estimation with interferometric aperture radar using fringe detection," in Proc. IGARSS'91, Helsinki, Finland, pp. 2173-2176, May 1991. [2] M . Schwabisch and D. Geudtner, "Improvement of phase and coherence map quality using azimuth pre-filtering: Examples from ERS-1 and X - S A R , " in Proc. IGARSS'95, Florence, Italy, pp. 205-207, 1995. [3] M . Seymour, Refining low-quality digital elevation models using synthetic aperture radar interferometry. Ph.D. thesis, University of British Columbia, 1999. [4] F. Gatelli, C. Prati, F. Rocca and others, "The wavenumber shift in SAR interferometry," IEEE Trans, on Geosci. Remote Sensing, vol. 32, pp. 855-865, July 1994. [5] D. Ghiglia, "IFSAR correlation improvement through local slope correction", in Proc. IGARSS'98, Seattle, U.S.A., pp. 2445-2447, 1998. [6] E. Trouve, J. Nicolas, H. Maitre, "Improving phase unwrapping techniques by the use of local frequency estimates", IEEE Trans. Geosci. Remote Sensing, vol. 36, pp. 1963-1972, November 1998. [7] H . Zebker and John Villasenor, "Decorrelation in interferometric radar echoes", IEEE Trans. Geosci. Remote Sensing, vol. 30, pp. 950-959, September 1992. [8] R. Bamler, N . Adam, G. Davidson, and D. Just, "Noise-induced slope distortion in 2D phase unwrapping by linear estimators with application to SAR interferometry" IEEE Trans. Geosci. Remote Sensing, vol. 36, pp. 913-921, May 1998. [9] M . Schwartz, L. Shaw, Discrete Spectral Analysis, Detection and Estimation. McGraw-Hill, 1988. 171 [10] J. Lee and others, "Intensity and phase statistics of multilook polarimetric and interferometric imagery", IEEE Trans. Geosci. Remote Sensing, vol. 32, pp. 1017-1028, September 1995. [11] E. Rodriguez and J. Martin, "Theory and design of interferometric synthetic aperture radar" Proc. Inst. Elect. Eng. F., vol. 139, no.2, pp.147-159. 1992. [12] J. Lee, T. Ainsworth and others, "Noise filtering of SAR interferometric phase images", in Proc. SPEE European Symp., Rome, Italy, pp. 735-742, Sept. 1994. [13] J. Lee, T. Ainsworth and others, " A new technique for noise filtering ,of SAR interferometric phase images" , IEEE Trans. Geosci. Remote Sensing, vol. 36, pp. 1456-1465, September 1998. [14] D. Geudtner, M . Schwabisch and R. Winter, "SAR Interferometry with ERS-1 data", in Proc. PIERS'94, pp. 821-825, 1994. [15] M . Hubig, " A maximum likelihood a priori filter for interferometric phase" in Proc. IGARSS'99, Hamburg, Germany, June 1999. [16] E. Trouve, M . Caramma and H . Maitre, "Fringe detection in noisy complex interferograms", Appl. Optics, vol. 35, pp. 3799-3806, July 1996. [17] I. Lim, T. Yeo, and others, "Phase Noise Filter for Interferometric SAR", in Proc. IGARSS' 97, Singapore, August 1997, pp. 445-447. [18] U . Spagnolini, "2-D phase unwrapping and instantaneous frequency estimation", IEEE Trans. Geosci. Remote Sensing, vol. 33, pp. 579-589, May 1995. [19] A. Candeias, J. Mura, J. Moreira and others, "Interferogram phase noise reduction using morphological and modified median filters" in Proc. IGARSS'95, Florence, Italy, pp. 166-168, 1995. [20] R. Goldstein, H. Zebker, and L . Werner, "Satellite radar interferometry: Two-dimensional phase unwrapping," Radio Sci., vol. 23, pp. 713-720, July 1998. [21] D. Ghiglia and L . Romero, "Robust two-dimensional weighted and un-weighted phase unwrapping that uses fast transform and iterative methods," J. Opt. Soc. Amer. A. , vol.11, no. 1, pp. 107-117, 1994. [22] J. Daugman, "Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by 2D visual cortical filters", J. Opt. Soc. Am., vol. 2 no 7, pp. 1160-1169, July 1985. 1 7 2 [23] H. Hsu, Fourier analysis. Simon & Schuster, New York, 1970. [24] A. Papoulis, Probability, Random Variables, and Stochastic Processes. McGraw-Hill, Inc. 1991. [25] R. Goldstein and C. Werner, "Radar Ice motion interferometry", in Proc. 3rd ERS symp., Florence, Italy, March 1997. [26] S. Haykin, Modern filters. MacMillan, Pub. Co., 1990. [27] G. Davidson and R. Bamler, "A multi-resolution approach to improve phase unwrapping", in Proc. IGARSS'96, Lincoln, U.S.A., pp. 2050-2053, 1996. [28] A. Bovik, P. Maragos, T. Quatieri, "AM-FM energy detection and separation in noise using multiband energy operators", IEEE Trans, on Signal Processing, vol. 41, No. 10, pp. 3245-3265, Dec. 1993. [29] J. Havlicek, A. Bovik, "Multi-component AM-FM image models and wavelet-based demodulation with component tracking", in Proc. IEEE Int. Conf. Image Processing, Austin, TX, Nov. 13-16, 1994, pp. 141-145. [30] B. Boashash, "Estimating and interpreting the instantaneous frequency of a signal - Part 2: Algorithms and applications", Proc. of IEEE, vol. 80, no. 4, pp. 540-568, April 1992. [31] P. Maragos, J. Kaiser and T. Quatieri, "Energy separation in signals modulations with application to speech analysis", IEEE Trans, on Signal Processing, vol. 41, no. 10, pp. 3024-3051, Oct. 1993. [32] D. Ghiglia and M. Pritt, Two-dimensional phase unwrapping. Theory, algorithms and software. John Wiley & Sons, Inc., 1998. [33] S. Madsen, "Estimating the Doppler centroid of SAR data", IEEE Trans, on aerospace and elect. Systems, vol. 25, no. 2, March 1989. [34] D. Perea-Vega and I. Cumming "Local frequency estimation in interferograms using a multi-band pre-filtering approach" in FRINGE'99 conference, Liege, Belgium, November 1999. [35] A. Bovik, A. Restrepo-Palacios, N. Gopal, and T. Emmoth "Localized measurement of emergent image frequencies by Gabor wavelets", IEEE Trans. on Information Theory, vol. 38, No. 2, pp. 691-712, March. 1992. 1 , 7 3 [36] D. Gabor, "Theory of comrnunication", J. Inst. Electr. Eng., vol. 93, pp. 429-457, 1946. [37] J. Curlander and R. Donough, Synthetic Aperture Radar: A signal processing approach. Kluwer, Boston, 1996. [38] S R T M mission JPL web site: http://www.jpl.nasa.gov/srtm/ 174
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Local fringe frequency estimation in synthetic aperture...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Local fringe frequency estimation in synthetic aperture radar interferograms using a multiband pre-filtering… Perea-Vega, Diego E. 2000
pdf
Page Metadata
Item Metadata
Title | Local fringe frequency estimation in synthetic aperture radar interferograms using a multiband pre-filtering approach |
Creator |
Perea-Vega, Diego E. |
Date Issued | 2000 |
Description | Synthetic aperture radar (SAR) interferometry is a technique for obtaining accurate elevation maps of the Earth from radar images. Local fringe frequency estimates are needed in several stages of the interferometry process. They are used to correct the effect of the topographic slope in the estimation of the interferogram coherence. They are needed to define the frequency center for adaptive band-pass filtering or to model the local phase in slope-correction filtering. Finally, accurate fringe frequency estimates facilitate the phase unwrapping process. In this work, I propose a new algorithm for local fringe frequency estimation in which the SAR interferogram signal is pre-filtered before the local frequency estimation is performed. This allows the use of a simpler and more efficient frequency estimator that operates at the pixel level. The proposed scheme shows advantages over other schemes because it achieves a better spacefrequency resolution and therefore tracks the topographic changes of the scene more accurately. The filters used in this work are modulated Gaussian functions with variable spatial aperture and bandwidth. In this way, the analysis window is adapted to the local characteristics of the signal at all samples. The variable-aperture filters are similar to the variable space-frequency domain filters used in wavelet analysis. I present results for synthetic and real SAR interferograms, as well as the performance of the proposed algorithm. An application of the frequency estimation method is developed for the noise filtering stage. A non-linear phase model is built to locally flatten the phase allowing the averaging of a higher number of samples without significant distortion. Simulations show that this alternative method achieves a better performance than two other reported methods when the interferogram coherence is moderately high. However, the alternative method does not solve the problem of phase discontinuities in the presence of topographically induced residues. |
Extent | 22756096 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-13 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065357 |
URI | http://hdl.handle.net/2429/10745 |
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 | 2000-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2000-0521.pdf [ 21.7MB ]
- Metadata
- JSON: 831-1.0065357.json
- JSON-LD: 831-1.0065357-ld.json
- RDF/XML (Pretty): 831-1.0065357-rdf.xml
- RDF/JSON: 831-1.0065357-rdf.json
- Turtle: 831-1.0065357-turtle.txt
- N-Triples: 831-1.0065357-rdf-ntriples.txt
- Original Record: 831-1.0065357-source.json
- Full Text
- 831-1.0065357-fulltext.txt
- Citation
- 831-1.0065357.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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0065357/manifest