A UNIFIED APPROACH TO THE GEOMETRIC RECTIFICATION OF REMOTELY SENSED IMAGERY by FRANK HAY-CHEE WONG B.Eng., McGill University, 1969 M.Sc., Queen's University, 1971 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY IN THE FACULTY OF GRADUATE STUDIES (Department of Computer Science) We accept this thesis as conforming to the required standard. THE UNIVERSITY OF BRITISH COLUBBIA APRIL 1984 fc) Frank Hay-Chee Wong, 1984 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f the requirements f o r an advanced degree a t the U n i v e r s i t y o f B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e copying o f t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the head o f my department o r by h i s o r her r e p r e s e n t a t i v e s . I t i s understood t h a t copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be allowed without my w r i t t e n p e r m i s s i o n . Department o f (^(^JH^J^ £d,<2A~C-J2. The U n i v e r s i t y o f B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date 27 fynlli ABSTRACT Many applications of remotely sensed d i g i t a l imagery require images that are corrected for geometric distortions. It i s often desirable to rectify different types of sa t e l l i t e imagery to a common data base. A high throughput rectification system is required for commercial application. High geometric and radiometric precision must be maintained. The thesis has accomplished the following tasks: 1. The sensors used to obtain remotely sensed imagery have been investigated and the associated geometric distortions inherent with each sensor are identified. 2. The transformation between image coordinates and datum coordinates has been determined and the values of the parameters in the transformation are estimated. 3. A unified rectification approach has been developed, for a l l types of remotely sensed d i g i t a l imagery, which yields a high system throughput. 4. Use of d i g i t a l terrain jnodels in the rectification process to correct for relief displacement has been incorporated. 5. An efficient image interpolation algorithm has been developed. i This algorithm takes into account the fact that imagery does not always correspond to sampling on a uniform grid. 6. The applications of rectified imagery such as mosaicking and multisensor integration have been studied. 7. Extension of the rectification algorithm to a future planetary mission has been investigated. The sensors studied include TIROS-N, Landsat-1, -2 and -3 multispectral scanners, Seasat synthetic aperture radar, Landsat-4 thematic mapper and SPOT linear array detectors. Imagery from the last two sensors i s simulated. i v TABLE OF CONTENTS PAGE ABSTRACT i i TABLE OF CONTENTS i v LIST OF TABLES x LIST OF FIGURES x i ACKNOWLEDGEMENT x v i 1. INTRODUCTION 1 1.1 Thesis Overview 1 1.2 Problem Statement 3 1.3 Background of Remote Sensing S a t e l l i t e s 6 1.4 Thesis Outline 11 2. INVESTIGATION OF AN APPROACH 13 2.1 Introduction 13 2.2 Geometric Transformation / Navigation 16 2.2.1 Weighted Arithmetic Mean 17 2.2.2 Polynomial Transformation 18 2.2.3 Spline Functions 19 2.2.4 General Time Series Model 20 2.2.5 Time Series Model with Recursive Estimator 21 2.3 Comparison of Geometric Transformation / Navigation Methods 23 2.3.1 D i s t r i b u t i o n of GCPs 23 2.3.2 Transformation Accuracy 25 2.3.3 Unique Advantages of Time Series Approach 26 2.3.4 Unique Advantages of a Recursive Estimator 27 2.3.5 Summary 27 2.4 Remapping and Interpolation 28 2.4.1 Background 28 2.4.2 B i l i n e a r Remapping 29 2.4.3 Two-Pass Processing 31 2.4.4 Three-Pass Processing 32 2.5 In t e r p o l a t i o n Kernel 38 2.5.1 Outline of Interpolation Methods 38 2.5.1.1 Nearest Neighbour and Linear Interpolations 38 2.5.1.2 Cubic S p l i n e 39 2.5.1.3 B-Spline 42 2.5.1.4 The Sampling Theorem 43 2.5.1.5 Cubic Convolution 44 2.5.1.6 Comparison of Interpolation Kernels 45 2.5.2 Preprocessing 47 2.5.2.1 Tabulation of Sine Function 47 2.5.2.2 Kernel S i z e Determination 51 2.5.2.3 Windowing 52 2.5.3 Nonuniformly Spaced Samples 55 2.6 System Data Flow 59 2.6.1 A u x i l i a r y Data Processing 61 2.6.2 R e l i e f Displacement Correction 62 2.7 Chapter Summary 65 3. DIGITAL TERRAIN MODELS 67 3.1 Introduction 67 v i 3.2 Ac q u i s i t i o n of DTMs 67 3.3 Types and Interpolation of DTMs 69 3.3.1 Irregular DTMs 70 3.3.2 Regular DTMs 71 3.3.3 Triangulated Irregular Network 72 3.4 DTM Accuracy 73 3.4.1 Accuracy Requirement i n Interpolated Coordinates 73 3.4.2 DTM Interpolation i n S c r o l l i n g Buffer 78 4. SCANNER IMAGERY 83 f- 4.1 Introduction 83 4.2 Overview of Geometric E f f e c t s 88 4.2.1 Common Dist o r t i o n s 88 4.2.2 Landsat-A TM Imagery S p e c i f i c D i s t o r t i o n s 92 4.2.3 SPOT Imagery S p e c i f i c D i s t o r t i o n s 92 4.2.4 Re l i e f Displacement 93 4.3 R e c t i f i c a t i o n Transformation 95 4.3.1 Closed Form Solution 95 4.3.2 Block Size 97 4.4 Recursive Estimator i n Navigation 98 4.4.1 State Vector Parameters 98 4.4.2 Extension to TIROS-N Imagery 101 4.4.2.1 Objective 101 4.4.2.2 TIROS-N Ch a r a c t e r i s t i c s 102 4.4.2.3 Orbit Accuracy Analysis 105 4.4.2.4 Attitude Accuracy Analysis 106 4.4.2.5 Experimental Results and Discussions 107 v i i 4.5 Experiments Interpolation i n the Presence of Scan Gap 114 4.5.1 Objective 114 4.5.2 Theory of Proposed Method 117 4.5.3 Implementation Considerations 122 4.5.4 Experiments with Sine Waves 124 4.5.5 Scan Gap Simulation i n Imagery 129 4.5.6 Experiments with Sine Wave Imagery 134 4.5.7 Synthesized TM Imagery 148 4.6 Experiments R e l i e f Displacement 154 4.6.1 Objective 154 4.6.2 Gaussian T e r r a i n with Landsat-4 TM Sensor Parameters 155 4.6.3 Real T e r r a i n with Landsat-4 TM Sensor Parameters 162 4.6.4 Real T e r r a i n with SPOT PLA Sensor Parameters 166 4.7 Experiments with Landsat-2 Imagery 170 5. SYNTHETIC APERTURE RADAR IMAGERY 181 5.1 Introduction 181 5.2 R e c t i f i c a t i o n Transformation 183 5.3 T e r r a i n Induced D i s t o r t i o n s 184 5.3.1 Rel i e f Displacement 184 5.3.2 Importance of Compression to Zero Doppler 189 5.3.3 Target Defocusing 191 5.3.4 Intensity Migration 191 5.4 Experiments with Seasat Imagery 194 5.4.1 Image Location Accuracy Study 194 5.4.2 Orbit Refinement 198 v i i i 5.4.3 R e c t i f i c a t i o n with DTM 201 5.5 Navigation by Automatic Focusing 208 6. APPLICATIONS 213 6.1 Introduction 213 6.2 Multisensor Integration 213 6.3 Image Mosaicking 227 6.3.1 Linear Interpolation 227 6.3.2 Feathering 229 6.3.3 F i l t e r i n g 231 6.3.4 Experiments 233 6.4 A Future Planetary Mission 237 6.4.1 Mission Description 237 6.4.2 Navigation 237 6.4.3 Venus Atmospheric Refraction 239 6.4.4 Use of Altimeter Data 242 6.4.5 D i g i t a l T e r r a i n Model 243 7. CONCLUSIONS AND RECOMMENDATIONS 245 7.1 Conclusions 245 7.2 Recommendations for Future Work 247 BIBLIOGRAPHY 253 APPENDIX A LANDSAT-4 SENSORS 259 APPENDIX B SPOT SENSORS 273 APPENDIX C RECURSIVE ESTIMATOR FILTER 277 APPENDIX D BILINEAR AND AFFINE TRANSFORMATIONS 282 APPENDIX E CLOSED FORM TRANSFORMATION FOR SCANNER IMAGERY 286 APPENDIX F THEORY OF SAR AND SAR DATA PROCESSING APPENDIX G CLOSED FORM RECTIFICATION TRANSFORMATION APPENDIX H REFRACTION OF ELECTROMAGNETIC WAVE LIST OF TABLES TABLE PAGE 1.1 Values of Re l i e f Displacement 5 3.1 Resolution and Worst R e l i e f Displacement Factor 75 3.2 DTM Resampling i n S c r o l l i n g Buffer (I) 81 3.3 DTM Resampling i n S c r o l l i n g Buffer (II) 82 4.1 TIROS-N C h a r a c t e r i s t i c s 103 4.2 Erro r Analysis of TIROS-N Navigation 115 4.3 RMS Error Versus Scan Gap Size Quantization 125 4.4 Attitude at 9 Points over the Image 178 5.1 Along Track and Across Track Errors 196 5.2 S e n s i t i v i t y of FM Rate to Parameter Variations i n Seasat Imagery 210 6.1 S e n s i t i v i t y of FM Rate to Parameter Variations i n VOIR Imagery 238 x i LIST OF FIGURES FIGURE PAGE 1.1 T y p i c a l S a t e l l i t e Orbit and Ground Track 7 1.2 An Evol u t i o n Cycle of S a t e l l i t e Sensor D e f i n i t i o n 10 2.1 Recursive Estimation 24 2.2 Two One-Dimensional Remapping and In t e r p o l a t i o n Processes 33 2.3 Three One-Dimensional Remapping and Int e r p o l a t i o n Processes 36 2.4 Coupling of Passes 37 2.5 Nearest Neighbour and Linear I n t e r p o l a t i o n Kernels and Frequency Responses 40 2.6 RMS E r r o r and Frequency Response for Cubic Convolution 46 2.7 RMS E r r o r and Frequency Response for Sine Function I n t e r p o l a t i o n 48 2.8 Approximation of Sine Function by Step Function 50 2.9 Examples of Autocorrelation Function 53 2.10 RMS Int e r p o l a t i o n Error Versus Kernel S i z e 54 2.11 Four-Point Sine Function and Its Frequency Response 56 2.12 Four-Point Sine Function with Kaiser Window and I t s Frequency Response 57 2.13 Sixteen-Point Sine Function with Kaiser Window and Its Frequency Response 58 2.14 System Data Flow 60 2.15 F i r s t Remapping and Interpolation Pass 64 3.1 DTM Accuracy Requirement for Landsat-4 TM Sensor 76 3.2 DTM Accuracy Reauirement for SPOT Sensors, Nadir Looking 76 3.3 DTM Accuracy Requirement for SPOT Sensors, Off Nadir Looking Angle = 26° 77 x i i 3.4 DTM Accuracy Requirement for Seasat SAR Sensor 77 4.1 Scanner Sensor System 84 4.2 T y p i c a l Mechanical Scanning Sensor 85 4.3 T y p i c a l E l e c t r o n i c Scanning Sensor 86 4.4 Earth Rotation E f f e c t 91 4.5 R e l i e f Displacement i n Scanner Imagery 94 4.6 Scan Vector Geometry 96 4.7 Accuracy Versus Block Size for Landsat Sensors 99 4.8 Accuracy Versus Block S i z e for SPOT Sensors 100 4.9 Resolution Versus Distance from Centre f o r TIROS-N Imagery 104 4.10 Uncertainties and Residual E r r o r Versus Number of GCPs Processed, V e r i f i c a t i o n Test 1 109 4.11 V e r i f i c a t i o n Test 2 111 4.12 Operational Run 111 4.13 Interpolation Kernel of Yen's Method 119 4.14 RMS Erro r Versus In t e r p o l a t i o n P o s i t i o n 127 4.15 Comparision of Methods — — RMS E r r o r Versus Scan Gap Siz e 130 4.16 E f f e c t of Overlap and Underlap on a Roadway 132 4.17 J i t t e r Model 133 4.18 Sine Wave Image 135 4.19 Distorted Sine Wave Image with Overlap 136 4.20 R e c t i f i e d Sine Wave Image (Overlap) with Sine I n t e r p o l a t i o n 137 4.21 R e c t i f i e d Sine Wave Image (Overlap) with Proposed Method 138 4.22 Radiometric Error i n R e c t i f i e d Sine Wave Image (Overlap) with Sine Interpolation 139 4.23 Radiometric Error i n R e c t i f i e d Sine Wave Image (Overlap) with Proposed Method 140 x i i i 4.24 Distorted Sine Wave Image with Underlap 142 4.25 R e c t i f i e d Sine Wave Image (Underlap) with Sine Interpolation 143 4.26 R e c t i f i e d Sine Wave Image (Underlap) with Proposed Method 144 4.27 Radiometric Er r o r i n R e c t i f i e d Sine Wave Image (Underlap) with Sine Interpolation 145 4.28 Radiometric Er r o r i n R e c t i f i e d Sine Wave Image (Underlap) with Proposed Method 146 4.29 E r r o r Magnitude 147 4.30 MSS Image Used to Synthesize TM Image 149 4.31 Radiometric Error i n R e c t i f i e d Synthesized TM Image with Sine Interpolation 150 4.32 Radiometric Er r o r i n R e c t i f i e d Synthesized TM Image with Proposed Method 151 4.33 RMS Error Versus Scan Gap Size i n Simulated TM Image 153 4.34 Gaussian T e r r a i n 156 4.35 Orthographic Image (Gaussian Terrain) 157 4.36 Distorted Image (Gaussian Terrain) 158 4.37 R e c t i f i e d Image (Gaussian Terrain) 159 4.38 Cross-Sectional Plot of a H i l l 160 4.39 Error Image (Gaussian Terrain) 161 4.40 DTM of S t . Mary Lake Area 163 4.41 Landsat-4 TM Orthographic Image 164 4.42 Landsat-4 TM R e l i e f Distorted Image 165 4.43 Landsat-4 TM R e c t i f i e d Image 167 4.44 Geometric R e c t f i c a t i o n Accuracy (Landsat-4 TM) 168 x i v 4.45 E r r o r Image (Landsat-4 TM) l b 9 4.46 SPOT Orthographic Image 171 4.47 SPOT R e l i e f Distorted Image 172 4.48 SPOT R e c t i f i e d Image 173 4.49 Geometric R e c t f i c a t i o n Accuracy (SPOT) 174 4.50 E r r o r Image (SPOT) 175 4.51 Landsat-2 MSS Raw Image 176 4.52 Residual E r r o r Versus Number of GCPs 178 4.53 R e c t i f i e d Landsat-2 MSS Image 180 5.1 Geometric Relation Between Slant Range and Ground Range 185 5.2 Plot of Ag/As Versus Range 186 5.3 R e l i e f Displacement i n Radar Imagery 187 5.4 Re l i e f Displacement and C r i t i c a l Angle 190 5.5 Predicted and Actual Target T r a j e c t o r i e s of an Elevated Target 192 5.6 Intensity Migration i n Radar Imagery 193 5.7 GCP E r r o r Vectors After Navigation 200 5.8 Block Size Versus Remapping Accuracy 202 5.9 Raw Seasat Image 204 5.10 R e c t i f i e d Seasat Image with DTM 205 5.11 R e c t i f i e d Seasat Image without DTM 206 5.12 E f f e c t of Azimuth Misfocusing 209 5.13 Navigation Loop Using Automatic Focusing 212 6.1 Seasat Component i n Landsat/Seasat Composite 215 6.2 Landsat Component i n Landsat/Seasat Composite 216 6.3 Landsat/Seasat Composite (12 Km x 12 Km) of the Ottawa Area 217 6.4 Histogram of Product of two Random Variables 219 6.5 Regions of Interest i n the Composite 221 6.6 Landsat/Seasat Composite (40 Km x 40 Km) of the Ottawa Area 223 6.7 Seasat Image of Vancouver, NTS Mapsheet 92G/3 224 6.8 Landsat/Seasat Composite of Vancouver, NTS Mapsheet 92G/3 225 6.9 NTS Mapsheet, 92G/3 226 6.10 Radiometric Mosaicking by Linear Interpolation 228 6.11 F i l t e r i n g Method of Image Mosaicking 232 6.12 Mosaic with no Overlap 234 6.13 Mosaic by Linear Interpolation 235 6.14 Mosaic by F i l t e r i n g 236 6.15 Refraction Geometry 240 6.16 Index of Refraction Versus A l t i t u d e 240 x v i ACKNOWLEDGEMENT I wish to thank my supervisor, Dr. Robert J . Woodham, for his valuable support, guidance, advice and encouragement throughout my years at the University of B r i t i s h Columbia. Dr. Woodham provided me with the d i g i t a l t e r r a i n model of the S t . Mary Lake area. He also imparted to me h i s knowledge of the synthesis of s a t e l l i t e borne imagery. Thanks are also due to Dr. Alan K. Mackworth, also of the Univ e r s i t y of B r i t i s h Columbia, f o r introducing me to Moire contourography and for widening my bird's eye view of a r t i f i c i a l i n t e l l i g e n c e . Manual d i g i t i z a t i o n of the d i g i t a l t e r r a i n models was performed by Mr. James L i t t l e and Mr. Ezio Cantanzariti of UBC using software o r i g i n a l l y supplied by Dr. Thomas Peucker of Simon Fraser U n i v e r s i t y . Mr. L i t t l e also provided me with a d i g i t a l t e r r a i n model of the Shawnigan Lake area. Sincere appreciation i s extended to the following colleagues at MacDonald, Dettwiler and Associates L t d . (MDA). I am extremely indebted to Dr. Robert Orth, without whose encouragement and permission I would never have undertaken t h i s research program. He has guided me t e c h n i c a l l y throughout t h i s research work. In addition he has con s i s t e n t l y shown deep concern f o r my career development at MDA and has taught me many other aspects of l i f e . Dr. Orth i s more than a friend and a colleague to me. I am equally indebted to Dr. John Bennett f o r his te c h n i c a l guidance i n the Synthetic Aperture Radar (SAR) area. He has continuously provided me, even though I am not formally working i n his section, with SAR projects which are relevant to the thesis t o p i c . I am also g r a t e f u l to the members of his SAR team, e s p e c i a l l y to Dr. Ian Cumming, Mr. Ron Fielden and Mr. Peter McConnell for t h e i r valuable discussions on the processing of SAR data and for providing me with SAR imagery. Thanks are also due to Mr. Daniel Friedmann for h i s discussion on i n t e r p o l a t i o n , and to Mrs. P a t r i c i a Palmer for developing part of the software used to correct TM imagery. Last, but c e r t a i n l y not l e a s t , I wish to acknowledge the Graduate Research and Engineering Technology award received from the Science Council of B r i t i s h Columbia and to thank MDA for allowing me the free and extensive use of i t s research f a c i l i t i e s . Consequently thanks are due to Dr. John MacDonald, president of MDA, for giving me the free p r i v i l e g e and for allowing me to pursue my academic int e r e s t s while employed at MDA. 1 CHAPTER I INTRODUCTION 1.1 Thesis Overview This thesis is concerned with the geometric rectification of remotely sensed d i g i t a l imagery including the incorporation of terrain data from d i g i t a l terrain models. The imagery studied includes that obtained from existing satellites in the Landsat series, the NOAA series and Seasat. It also includes that to be obtained from the proposed SPOT s a t e l l i t e . Landsat imagery and NOAA are obtained by mechanical scanning sensors while synthetic aperture radar (SAR) imagery is obtained by a radar sensor. SPOT imagery w i l l be obtained by a linear array sensor. The thesis develops a framework for an integrated image processing system with the f l e x i b i l i t y to deal with different sensors, to yield high throughput, to retain high geometric and radiometric f i d e l i t y and to be operationally viable. The thesis combines a review and an evaluation of previous work, extending i t when applicable, with the development of novel techniques to deal with problems previously not encountered. Specifically, the following is dealt with in this thesis: (a) The sensors used to obtain image data are investigated and the associated geometric distortions inherent with each sensor are identified. Particular attention is given to Landsat-4, SPOT and Seasat SAR. The geometric distortions in Landsat-4 Thematic Mapper (TM) imagery and in SPOT panchromatic linear array (PLA) and multispectral linear array (MLA) imagery are shown to be d i f f e r e n t from those i n Landsat-1, -2 and -3 multispectral (MSS) imagery. The geometric d i s t o r t i o n s i n Seasat SAR imagery are also d i f f e r e n t from those i n MSS and l i n e a r array imagery. These d i s t o r t i o n s depend on the method of processing the Seasat SAR data. (b) The transformation between image coordinates and datum coordinates i s determined. The values of the parameters i n the transformation are estimated. The o r b i t and attitude parameters are of p a r t i c u l a r importance. Their values can be determined with the aid of ground co n t r o l points (GCPs). In SAR imagery, only the o r b i t parameters are important. A method to determine t h e i r values by automatic focusing i s discussed. (c) A u n i f i e d r e c t i f i c a t i o n approach i s developed, for a l l types of remotely sensed d i g i t a l imagery, which can be implemented i n a conventional image processing system and which y i e l d s a high system throughput. High throughput i s achieved with one-dimensional processing of imagery which i s represented by a two-dimensional array of d i g i t a l values. (d) D i g i t a l t e r r a i n models (DTMs) are incorporated into the r e c t i f i c a t i o n process to correct f or r e l i e f displacement. Correction f or r e l i e f displacement i s also achieved by one-dimensional processing. (e) An e f f i c i e n t image i n t e r p o l a t i o n algorithm i s developed. This algorithm takes i n t o account the fact that imagery does not always correspond to sampling on a uniform g r i d . This i s required f o r Landsat-4 Thematic Mapper (TM) imagery where scan gaps are present. 3 (f) Two applications of rectified imagery are presented: image mosaicking and multisensor integration. In image mosaicking, two images bordering on the same area are merged into a single image. In multisensor integration, two images from different sensors, Seasat SAR and Landsat-2 MSS, are combined into a single image. (g) Extensions to the re c t i f i c a t i o n approach are considered for a future planetary mission. 1.2 Problem Statement Geometric rectification i s the warping of raw image data to a predefined datum such as a particular map projection. Many applications of remotely sensed imagery require rectification to a common datum. One example is the production of Seasat/Landsat composites (to be shown in Chapter 6) which combine the high resolution of Seasat with the spectral content of Landsat. Rectification is d i f f i c u l t since different s a t e l l i t e sensors have different orbit parameters and geometric attributes. Rectification of s a t e l l i t e imagery involves two steps: (a) Geometric transformation determination of transformation between input imagery and datum. The transformation depends on the following parameters: s a t e l l i t e platform and orbit data, terrain height, earth curvature and rotation rate. The transformation i s normally established with the aid of GCPs. (b) Interpolation sampling of image data using a specific 4 kernel; e.g., sp l i n e or sine function. The method of i n t e r p o l a t i o n depends upon the geometric transformation algorithm. R e c t i f i c a t i o n has been dealt with by many other authors [BAKE75, EBNE76, HORN79, KONE76, krat72, KRAU78, LARS 80, LECK80, ORTH78a, ORTH78b, OTEP78 and SIM075a]. Previous work c o l l e c t i v e l y shares the following d e f i c i e n c i e s when considered f o r the general case: (a) R e l i e f displacement i s not accounted f o r . Table 1.1 shows the p i x e l s h i f t s that would appear due to r e l i e f displacement f o r various s a t e l l i t e imagery at maximum scan angle. While r e l i e f displacement i s not s i g n i f i c a n t i n Landsat-1, -2, and -3 imagery, i t i s s i g n i f i c a n t i n Seasat SAR imagery and Landsat-4 TM imagery and w i l l be s i g n i f i c a n t i n SPOT imagery. T e r r a i n information i s required f o r geometric transformation. (b) The transformation and i n t e r p o l a t i o n algorithms presented do not address the issue of system throughput. Throughput i s of prime importance i n TM and SPOT imagery since each scene contains or w i l l contain as much as 6000 l i n e s by 6000 pix e l s per l i n e compared to approximately 3000 l i n e s by 3000 p i x e l s per l i n e i n an MSS scene. (c) The r e c t i f i c a t i o n algorithms are sensor and s a t e l l i t e s p e c i f i c . Imagery from d i f f e r e n t sensors cannot be r e c t i f i e d i n a s i n g l e shared system. (d) The geometric transformation algorithm i s often of generalized form (e.g., polynomial) and does not model the imaging geometry TABLE 1.1 VALUES OF RELIEF DISPLACEMENT 5 S a t e l l i t e s e n s o r A l t i t u d e Maximum scan a n g l e Nominal p i x e l s p a c i n g T a r g e t e l e v a t i o n f o r 1 p i x e l d i s p l a c e m e n t L a n d s a t - 1 , -2 -3 MSS 902 km 5.78° 58 m 507 m L a n d s a t - 4 MSS 705 km 7.47° 58 m 406 m La n d s a t - 4 TM (Bands 1, 2, 3, 4, 5 and 7) 705 km 7.47° 30 m 210 m La n d s a t - 4 TM (Band 6) 705 km 7.47° 120 m 840 m SPOT MLA ( n a d i r mode) 820 km 2.12° 20 m 472 m SPOT MLA ( s i d e l o o k i n g mode) 820 km 28.12° 26 m 44 m SPOT PLA ( n a d i r mode) 820 km 2.12° 10 m 236 in SPOT PLA ( s i d e l o o k i n g mode) 820 km 28.12° 13 m 22 m S e a s a t SAR 820 km 20.00° 25 m 11 m 6 ex p l i c i t l y . Subpixel rectification accuracy requires high quality GCPs which are d i f f i c u l t to establish. Algorithms which require a small number of GCPs are preferable, (e) The problem of scan gap in TM imagery is new. This thesis addresses these deficiencies e x p l i c i t l y . Each step is considered and the effect of terrain, i f any, is investigated. 1.3 Background of Remote Sensing Satellites S a t e l l i t e imagery i s obtained from a platform orbiting around the earth as shown in Figure 1.1. The combined effect of earth rotation and sa t e l l i t e velocity allows the s a t e l l i t e to give repetitive coverage of the entire globe. This produces a large volume of remotely sensed data in a reasonably short time, and above a l l , at a low cost. The f i r s t generation of remote sensing sat e l l i t e s began in the 1960s with the TIROS and NIMBUS series [NASA76a] which offered a spatial resolution of about 1 km. Despite the low spatial resolution, the imagery obtained demonstrated the potential of viewing the earth from space. The second generation of satellites was the Landsat [NASA76b] series launched in the 1970s. Each s a t e l l i t e recorded imagery in four spectral bands with resolution 1 greater than 100 metres. Currently, data from these satellites are s t i l l being exploited. Landsat imagery has demonstrated the u t i l i t y of remotely sensed, s a t e l l i t e data in earth Resolution increases with a decrease in i t s numerical value. 7 (b) GROUND TRACK FIGURE 1,1 TYPICAL SATELLITE ORBIT AND GROUND TRACK 8 resources applications such as geology, a g r i c u l t u r e , f o r e s t r y and oceanography. A p p l i c a t i o n i s also found i n cartography. The Landsat s e r i e s , although not designed f o r mapping, was capable of generating imagery for presentation at a scale of 1:250,000. P r e c i s i o n processing of Landsat imagery has been a f i r s t step i n mapping from s a t e l l i t e borne imagery. Also i n the l a t e 1970s, a c i v i l i a n SAR system (Seasat [MDA76]) was deployed and t h i s gave a view of the earth i n the radio frequency portion of the electromagnetic spectrum. Seasat SAR data has been p r e c i s i o n processed to give 25 metres r e s o l u t i o n . Although Seasat was operational f o r only four months, a l l the data c o l l e c t e d has yet to be completely processed. Landsat-D [GSFC82] was launched i n July 1982 becoming Landsat-4. A key feature i n Landsat-4 i s i t s b i d i r e c t i o n a l scanning TM sensor 1 with seven s p e c t r a l bands and a s p a t i a l r e s o l u t i o n of 30 metres i n s i x of these channels. The data should support mapping at a scale of 1:100,000 and w i l l be used for thematic i n t e r p r e t a t i o n . Due to the b i d i r e c t i o n a l scanning of the TM sensor and i t s high r e s o l u t i o n , the r e s u l t i n g imagery s u f f e r s from geometric d i s t o r t i o n more severely than imagery obtained from the e a r l i e r Landsat s a t e l l i t e s . The next generation of s a t e l l i t e s w i l l include SPOT2 [CNES81] which w i l l provide a s p a t i a l r e s o l u t i o n from 10 to 27 metres. SPOT's design also provides stereo c a p a b i l i t y . Large areas of the earth remain to be mapped p r e c i s e l y and SPOT imagery w i l l be used to prepare and update topography maps at a scale of 1:100,000. SAR systems have been carried on 1 See Appendix A. 2 See Appendix B. 9 the space s h u t t l e s . In Canada, a SAR system (RADARSAT) i s planned for i c e study. The trend to use s a t e l l i t e imagery i n mapping i n addition to conventional a e r i a l photography i s gaining support. S a t e l l i t e s are a more cost e f f e c t i v e source of cartographic information compared to conventional a e r i a l photography. The d i g i t a l nature of s a t e l l i t e imagery can also render i t easier to handle and to process. Colvocoresses [COLV79] has advocated a mapping s a t e l l i t e system (MAPSAT) to be deployed i n the 1980s. It w i l l also have the stereo c a p a b i l i t y . Both SPOT and the proposed MAPSAT w i l l use s o l i d state l i n e a r array detectors which can produce imagery of high geometric f i d e l i t y . Another synthetic aperture radar system has been proposed by NASA to map the surface of the planet Venus [JAME82]. The proposed system w i l l be launched i n the la t e 1980s or i n the 1990s. In a short time span, s a t e l l i t e imaging systems have evolved from low s p a t i a l r e s o l u t i o n to high s p a t i a l resolution, have increased the number of s p e c t r a l channels and s p e c t r a l resolution, have extended outside the v i s i b l e portion of the spectrum to include near i n f r a r e d , thermal in f r a r e d and radio frequencies, and have begun to use s o l i d state l i n e a r array detectors. This evolution proceeds i n cycles as i l l u s t r a t e d i n Figure 1.2. Each cycle s t a r t s with the user d e f i n i t i o n of sensor requirements, in c l u d i n g the number of s p e c t r a l bands, the s p e c t r a l response and s p a t i a l r e s o l u t i o n . Then a s a t e l l i t e sensor i s developed and launched. Data acquired by the sensor are transmitted to a ground s t a t i o n where they 10 \ \ \ SATELLITE SENSOR DEFINITION SATELLITE MISSION GROUND STATION \ \ \ \ \ \ \ IMAGE PROCESSING DISTRIBUTION TO USER FIGURE 1.2 AN EVOLUTION CYCLE OF SATELLITE SENSOR DEFINITION 11 are received and processed. F i n a l l y , image products are d i s t r i b u t e d to the user. Experience gained using the imagery allows the user to define requirements f o r the next generation of s a t e l l i t e s . 1.4 Thesis Outline This introductory chapter has defined the problem. The rest of the the s i s i s organized as follows: Chapter 2 gives the methodology for the r e c t i f i c a t i o n of remotely sensed imagery. It presents methods for geometric transformation and i n t e r p o l a t i o n and demonstrates the use of DTMs i n the r e c t i f i c a t i o n process. The methods y i e l d high throughput. The o v e r a l l geometric transformation i s c o e f f i c i e n t driven. Image i n t e r p o l a t i o n i s performed by one-dimensional processing. D i f f e r e n t types of imagery are handled. Chapter 3 describes DTMs, t h e i r a c q u i s i t i o n , types, accuracies and methods of i n t e r p o l a t i o n . Chapter 4 discusses scanner and l i n e a r array imagery. It contains a review of related work i n the area of image r e c t i f i c a t i o n and how GCPs are used. The attitud e modelling method (the at t i t u d e information i s i n the r e c t i f i c a t i o n geometric transformation) i s tested on both Landsat and TIROS-N imagery. F i n a l l y , synthesized Landsat-4 and SPOT imagery i s used to show the s p e c i f i c geometric d i s t o r t i o n s including r e l i e f displacement e f f e c t s . R e c t i f i c a t i o n s using the proposed algorithm are performed on the synthesized imagery as well as on Landsat-2 MSS imagery. 12 Chapter 5 discusses SAR imagery. Unlike scanner and linear array imagery, the raw data of SAR imagery need to be processed to be presented in image form. This chapter starts with an introduction to the theory of SAR and then i t outlines the geometric transformation. It also shows how GCPs are used to update the s a t e l l i t e orbit data. Finally, the problem of relief displacement is examined and results of rectifying a Seasat SAR scene using a DTM of the same area are presented. Chapter 6 presents some applications with rectified imagery such as image mosaicking and multisensor integration. It also extends the method to imagery obtainable from a proposed planetary mission. Chapter 7 concludes the thesis with discussions of the method presented, the experimental results and recommendations for future research work. 13 CHAPTER 2 INVESTIGATION OF AN APPROACH 2.1 Introduction This chapter f i r s t discusses and compares various r e c t i f i c a t i o n approaches, then proposes and investigates one which can handle remotely sensed imagery from d i f f e r e n t s a t e l l i t e s . R e c t i f i c a t i o n i s the warping of an input image to an output g r i d . The input image i s seen by a s a t e l l i t e sensor and hence contains inherent geometric d i s t o r t i o n s caused by earth curvature and r o t a t i o n , s a t e l l i t e a t t i t u d e v a r i a t i o n across the image, panoramic and r e l i e f d i s t o r t i o n s . The output gri d i s a user defined map p r o j e c t i o n . R e c t i f i c a t i o n involves two steps: geometric transformation and i n t erpolation. Since the f i n a l r e s u l t of r e c t i f i c a t i o n i s an image i n the output coordinates, i t i s necessary f o r the geometric transformation to map the output image coordinates (x,y) to the input image coordinates (u,v). The transformation i s represented by the general form: The functions f and g can be represented i n many forms, the simplest being a p a i r of low order polynomials and the more elaborate being a p a i r of complicated algebraic functions e x p l i c i t l y taking parameters of the imaging geometry and earth r o t a t i o n i n t o account. Regardless of the representation, function parameters and c o e f f i c i e n t s need to be determined u = f a ( x , y ) (2.1a) v = g a(x,y) (2.1b) 14 and t h i s normally requires the use of GCPs. In t h i s t h e s i s the determination of these parameters i s c a l l e d n a v i g a t i o n 1 . A t y p i c a l example of parameters i s the set of a t t i t u d e angles of the s a t e l l i t e p l a t f o r m . The output image may contains m i l l i o n s of data p o i n t s . To perform the mapping defined by Equations 2.1a and b f o r a l l the data points may be too time consuming f o r p r a c t i c a l purposes. Therefore i t i s u s e f u l to f i n d a f a s t mapping: u = f b ( x , y ) (2.2a) v = 8 b(x,y) (2.2b) such that f, and g. approximate f and g and t h i s new form renders a b & b a °a f a s t computation. In t h i s t h e s i s the establishment of f ^ and g^ i s c a l l e d remapping. I t must be emphasized that f and g^ are not i n the same f u n c t i o n a l form as f and g , e s p e c i a l l y when f and g are complicated 3 3 3 3 a l g e b r a i c f u n c t i o n s . But f and g must be e s t a b l i s h e d f i r s t and they are 3 3 used to determine f, and g, b b . I n t e r p o l a t i o n determines an i n t e n s i t y value at each output image i n t e g r a l coordinate ( x , y ) . I n t e g r a l values of (x,y) do not n e c e s s a r i l y map i n t o i n t e g r a l values of ( u , v ) . I n t e r p o l a t i o n determines the i n t e n s i t y at each mapped, but n o n - i n t e g r a l ( u , v ) . The convolution method i s used i n t h i s t h e s i s . I f the input imagery data are band l i m i t e d , then the c o n v o l u t i o n k e r n e l ( a l s o c a l l e d the i n t e r p o l a t i o n k e r n e l ) should correspond t o that of an i d e a l low pass f i l t e r [LATR65]. Such a f i l t e r can only be r e a l i z e d w i t h an i n f i n i t e k e r n e l s i z e and hence i s not 1 The term " n a v i g a t i o n " i s used i n t h i s t h e s i s to avoid a lengthy phrase "Using GCPs to determine the parameter values i n the geometric t r a n s f o r m a t i o n given by f and g i n Equation 2.1". 3 3 15 p r a c t i c a l for d i r e c t implementation. A p r a c t i c a l kernel has to be designed so that high radiometric accuracy i s retained i n the i n t e r p o l a t i o n process. Image i n t e r p o l a t i o n i s a two-dimensional convolution process ( i . e . , the convolution kernel i s a two-dimensional s p a t i a l f i l t e r ) . This can be a slow process for two reasons: the number of arithmetic operations i s large and many phys i c a l image read/write operations may be required. However by c a r e f u l l y choosing the remapping functions f ^ and g^ i n t e r p o l a t i o n i s reduced to several one-dimensional processes while r e t a i n i n g the same degree of radiometric accuracy as i n the two-dimensional i n t e r p o l a t i o n case. This chapter presents a navigation algorithm which y i e l d s high geometric p r e c i s i o n and which uses r e l a t i v e l y few GCPs per image scene. T h i s chapter also presents a fast i n t e r p o l a t i o n algorithm, including the use of DTM data. Remapping i s table (or c o e f f i c i e n t ) driven so that i t has the f l e x i b i l i t y to handle a v a r i e t y of s a t e l l i t e imagery. Remapping and i n t e r p o l a t i o n are implemented i n several one-dimensional processing passes. The result i s a high throughput system for r e c t i f y i n g various kinds of s a t e l l i t e borne imagery, including r e l i e f displacement c o r r e c t i o n . Furthermore, radiometric f i d e l i t y i s preserved since one-dimensional processing allows the use of a large i n t e r p o l a t i o n kernel. The image i s stored i n a conventional manner ( i . e . , successive l i n e s of imagery are stored i n consecutive records on disk or tape). Sections 2.2 and 2.3 present and compare various methods. Section 2.4 discusses remapping and i n t e r p o l a t i o n . navigation Section 2.5 16 develops a suitable choice of i n t e r p o l a t i o n kernel. F i n a l l y , Section 2.6 concludes with an i l l u s t r a t i o n of the flow of data (image and DTM) through a t y p i c a l system configured to use the algorithms developed. Though the geometry of d i f f e r e n t sensors i s d i f f e r e n t , the methodology and software performing the r e c t i f i c a t i o n i s the same. 2.2 Geometric Transformation / Navigation Procedures developed for r e c t i f i c a t i o n of s a t e l l i t e imagery are based on the computation of geometric c o e f f i c i e n t s of mathematical models that characterize the d i s t o r t i o n s i n the image. The d i s t o r t i o n s are represented by functions f and g i n Equations 2.1a and b. After obtaining values f o r these c o e f f i c i e n t s with the aid of GCPs i n the navigation process, the mapping function to be used i n the remapping process i s determined. In the l i t e r a t u r e , numerous algorithms have been reported [BAKE75, BERN76, EBNE76, HORN79, KONE76, KRAT72, KRAU78, LARS80, LECK80, LITT'80, 0RTH78a, 0RTH78b, and SIM075a] , In general, these algorithms e i t h e r use a polynomial approach or set up an attitu d e time series model. This section describes some of the approaches, ,and i d e n t i f i e s the p a r t i c u l a r approach adopted i n t h i s t h e s i s . The transformations i n Sections 2.2.1 to 2.2,5 r e f e r to f and e a 6 a in Equation 2.1. 17 2.2.1 Weighted Arithmetic Mean The weighted arithmetic mean method calculates f or any i n t e r p o l a t i o n point (x,y) i n the output grid the displacements Ax = u - x and Ay = v -y, each as a function of displacements of a l l GCPs over the en t i r e image. S p e c i f i c a l l y : 1 n A x = - i W (x,y,x ,y ) Ax (2.3a) K i=l I n Ay = - • Z W (x,y,x ,y ) Ay (2.3b) K 1=1 where i varies over a l l GCPs, hx^ and Ay^ are displacements of the i - t h GCP, n i s the t o t a l number of GCPs, W_^ i s a weighting function and should be a monotonically decreasing function of distance between the i - t h GCP and the i n t e r p o l a t i o n point (x,y), and K i s the normalization constant given by: n K = £ W. • ( (2.A) 1 = 1 1 This i s known as Shepard's method [SCHU76a] i n the i n t e r p o l a t i o n l i t e r a t u r e . While t h i s method i s simple, i t requires that the distance weight be independent of the d i r e c t i o n and p o s i t i o n of the i n t e r p o l a t i o n point [BAKE75 and KONE76]. This requirement i s not n e c e s s a r i l y obeyed i n s a t e l l i t e imagery. 18 2.2.2 Polynomial Transformation The polynomial approach [BERN76, KONE76, LECK80] extrapolates the displacement made at the GCPs to a l l other points in the image by modelling the distortion function as two polynomial functions of position within the image. The coefficients of the polynomials are estimated by least square methods. In effect, this method is a "rubber sheet f i t t i n g " technique between the image and the reference projection. This approach is independent of sa t e l l i t e attitude, scanner geometry and earth geometry. The polynomial can be written as : N N-p u = Z Z a x P y q (2.5a) p=0 q=0 P q N N-p v = Z I b x p y q (2.5b) p=0 q=0 P q where N is the order of the polynomial transformation. The GCPs are used to compute the coefficients a and b by a least squares method. pq pq Sometimes the image is semi-corrected .for system errors such as earth rotation, panoramic distortion, nonlinear sensor sweep and scan skew. Let the semi-corrected image coordinates be (x',y'). In this case, polynomials are applied to displacement errors between (x,y) and (x',y'): N N-p Ax = x' - x = Z Z a x P y q (2.6a) p=0 q=0 P q N N-p AY = y' - y = Z . Z b p q x P y q (2.6b) p=0 q=0 19 The weighted arithmetic mean of Section 2.2.1 i s a s p e c i a l case of the polynomial transformation as can be demonstrated by rewriting Equations 2.3a and b as follows: n u = x + Z f.(x,y) (u - x ) (2.7a) 1=1 n and v = y + Z f.(x,y) (v.- y.) (2.7b) 1=1 1 1 1 n with f.(x,y) = W. / Z W (2.7c) 1 1 i=l The powers of the polynomials i n x and y then depend on the functional forms of f ^ ( x , y ) . 2.2.3 Spline Functions Transformation by polynomials over the entire image often has the disadvantage that higher order polynomials are needed to ensure good f i t . It i s possible to achieve the same degree of f i t with lower order polynomials by f i r s t segmenting the image into separate regions. A polynomial i s determined f o r each region. These polynomials can be made continuous up to order N-l (N i s the polynomial order) at the region boundaries. These polynomials define sp l i n e functions whose c o e f f i c i e n t s are determined using GCPs and a least squares f i t t i n g method. 20 2.2.A General Time Series Model The time series approach models the displacement terms mathematically based on platform and sensor geometry. This t y p i c a l l y involves trigonometry and other complex mathematical functions. If a l l the parameters are known, the s a t e l l i t e data can be geometrically corrected by transformation equations of the form 1 : x = x [p , p„, p , t(u,v)] (2.8a) i I n y = y [Pj, P2> P n» t(u,v)] (2.8b) where t i s the imaging time of an input p i x e l having coordinates (u,v) and the P^'s (i=l to n ) are parameters of scanner geometry, earth ro t a t i o n , s a t e l l i t e o r b i t and a t t i t u d e , mirror sweep rate and the l i k e . Unfortunately, a l l of the parameters are not known a p r i o r i . Unknown parameters can be expressed as polynomials i n time. For example, the i - t h parameter i n the group i s represented by an N^-th degree polynomial: N i 1 p = I k t 1 (2.9) 1=0 The c o e f f i c i e n t s k are estimated using GCPs i n the following manner [K0NE76] using Newton's method: (a) Assume i n i t i a l values of the c o e f f i c i e n t s . (b) Compute the transformation error Ax and Ay for a l l GCP locations.; (c) At the same GCP locations, evaluate the p a r t i a l derivatives of x and y with respect to the c o e f f i c i e n t s . (d) Steps (b) and (c) determine a set of simultaneous l i n e a r equations i n A k ^ which can be solved by a least squares method. (e) Update k.^ by -Ak l For the discussion of determining the p. values, i t i s more convenient to express x and y as functions of u and v than v i c e versa. 21 ( f ) Repeat Steps (b) to (e) u n t i l the values of the c o e f f i c i e n t s i n the time series no longer change enough to affect the o v e r a l l geometric accuracy of the transformation. 2.2.5 Time Series Model with Recursive Estimator Another approach [MDA78, CAR075 and ORTH78a] i s s i m i l a r to that of Section 2.2.4 except that the GCPs drive a recursive estimator to r e f i n e the time series c o e f f i c i e n t s . The recursive estimator i s a s p e c i a l case of the Kalman f i l t e r [Appendix A, MDA78, DUPL67, MILL71, LE0N70 and UCLA79] and i s discussed i n Appendix C. The recursive estimate incorporates the GCPs one at a time into the estimator, updating the error model (time series) a f t e r each. The f i r s t GCP allows removal of some of the image errors, and each following GCP r e f i n e s the geometric p r e c i s i o n . More importantly, the uncertainty i n the predicted l o c a t i o n of each new GCP decreases due to the incrementally improved estimate. Th i s has d i s t i n c t operational advantages, both i n searching for the next GCP l o c a t i o n and i n detecting pointing errors associated with the current GCP. The estimator i s a set of recursive equations for optimally, i n the least squares sense,estimating the state vector containing va ri ab le s which are obtained from noise corrupted measurement data. In the a p p l i c a t i o n here, the state vector consists of the time series c o e f f i c i e n t s . The recursive equations (derived i n Appendix C) are summarized i n what fol l o w s . 22 An earth centred r o t a t i n g (ECR) coordinate system i s chosen to reference a GCP on the earth e l l i p s o i d . In t h i s coordinate system the o r i g i n i s at the earth's centre, the x-axis points to 0° longitude, the z-azis points to the north pole and the y-axis completes a right-handed coordinate system. Let X be the state vector, P be the covariance matrix of the state vector, E be the error covariance matrix projected on the ECR coordinates i n the GCP measurement, H be the measurement matrix to be defined and I be the i d e n t i t y matrix. The recursive equations are then given by: (a) Compute optimal weight b = P ( i ) H T [H P ( i ) H T + E ] " 1 (2.10) (b) Make optimal estimate *•(!) = m - l ) + b [R ( i ) - R (Y-(i-l))] (2.11) " -g -g -(c) Update covariance matrix P ( i ) = ( I - bH) P ( i ) (2.12) Here, R ( i ) are the GCP's ECR coordinates ( i . e . , from a mapsheet), and R g(|) a r e the estimated ECR coordinates f o r the same GCP using Y i n the transformation. The components of E represent GCP designation accuracies i n the three ECR coordinate axes. The matrix E i s taken to be diagonal since the designation accuracies i n the axes are uncorrelated. Hence: E ( i , j ) = 5 2 (2.13) Therefore the value of £ r e f l e c t s the confidence l e v e l i n the GCP 23 measurement, and this i s determined by the operator's pointing accuracy on the displayed image and the GCP accuracy on the mapsheet. The components of H are the partial derivatives of with respect to the state vector parameters. The recursive estimation procedure is shown in Figure 2.1. Landsat-1, -2 and -3 imagery is routinely rectified using this technique [ORTH77, ORTH78a, and ORTH78b] and subpixel accuracy is achievable using 3 to 12 GCPs in a Landsat scene. The number of GCPs required decreases with increased r e l i a b i l i t y in GCP map coordinates and GCP identification. Here the parameters to be estimated are s a t e l l i t e attitude angles (pitch, r o l l and yaw) and s a t e l l i t e altitude. 2.3 Comparision of Geometric Transformation / Navigation Methods In this section comparison is made between the polynomial and the time series methods. In general, one concludes that the latter technique is superior to the former. Comparison along specific dimensions i s presented i n the subsections which follow. 2.3.1 Distribution of GCPs Transformation by the polynomial method has the advantage that no prior information on the geometric distortions is included. Hence, no trigonometric or other complex algebraic transformation equations are involved. However, while polynomials may be made to f i t well at the GCPs from which the coefficients are determined, they can deviate significantly i n areas where no GCPs are present. Therefore, a uniform GCP L i n e / p i x e l c o o r d i n a t e s .Measured ECR c o o r d i n a t e s Updated * J R e c u r s i v e f i l t e r FIGURE 2.1 RECURSIVE ESTIMATION 25 spatial distribution of the GCPs i s v i t a l to the overall accuracy of the polynomial transformation method. In the case of the attitude time series model, the GCP distribution i s constrained only by the requirement that yaw is best estimated by GCPs at the frame edge, since yaw corresponds to rotation about the frame centre lin e . Similarly, the effect of r o l l i s most noticeable at the edges of a scan line or frame. Finally the maximum effect of pitch is to be found at the top and bottom of a frame. An equivalent polynomial order for the time series method is reflected i n the complexity of Equations 2.8a and b. While this is undoubtedly a high order transformation, advantage is taken i n the formulation of the interdependence of coefficients to reduce the total number of degrees of freedom. The fact that this interdependence is taken into account explicitly in the time series method guarantees mapping accuracy i n areas with no GCPs, insofar as no discontinuous changes in the time series parameters can occur. 2.3.2 Transformation Accuracy Image rectification accuracy is often required to subpixel precision, RMS, across the scene. Simon and Caron [SIM075] point out that the polynomial approach cannot be expected to reach this precision, with few (less than 20) GCPs. To achieve subpixel rectification accuracy for Landsat-1, -2 and -3 imagery, each attitude component should be known to the nearest tenth of a milliradian. To this end, Caron and Simon 26 developed a l i n e a r recursive estimator f o r the time series model. The transformation accuracy obtained i s at the subpixel l e v e l with about 10 GCPs, under favourable circumstances ( i . e . , easy i d e n t i f i c a t i o n of GCP and r e l i a b l e GCP map coordinates). Knoecny's [K0NE76] survey of the polynomial approach shows that to achieve the same degree of accuracy, a 12 degree polynomial and 64 GCPs would have to be used. Konecny's survey also shows that i f an a f f i n e transformation i s used, 15 GCPs would give an RMS error of 160 metres i n Landsat imagery. By increasing the number of GCPs to 33, the RMS error drops to 150 metres. 2.3.3 Unique Advantages of Time Series Approach The time series approach allows c e r t a i n parameters to be determined e x p l i c i t l y (rather than being embedded i n polynomial c o e f f i c i e n t s ) . T h i s i s advantageous when those parameters are required f or other considerations. In the case of scanner imagery, r o l l also a f f e c t s r e l i e f displacement (see Section 4.4.4). Hence, for t h i s c o r r e c t i o n , i t i s important' to have the r o l l information i n add i t i o n to t e r r a i n data. The polynomial approach does not give the r o l l angle whereas the time seri e s modelling approach gives each attitude angle as a function of time. Therefore, f o r each p i x e l p o s i t i o n i n an image, the r o l l angle can be determined. 27 2.3.4 Unique Advantages of a Recursive Estimator The recursive estimator approach has operational advantages derived from the sequential measurement of each GCP followed by an estimator update. This sequential structure allows an estimate of the accuracy of a GCP measurement before i t i s used to update the current a t t i t u d e estimates. Erroneous GCP measurements can be rejected. Moreover, as errors diminish at each stage, i t becomes increasingly easy to locate new GCPs i n the image. No such information can be generated by the polynomial or time series model with an i t e r a t i v e approach. Also i n the recursive estimator approach, the GCP deviations can be l i s t e d following the f i t t i n g procedures. Excessive deviations are noted by the operator, who then must edit the input data to the estimator to delete undesirable GCPs. The operator would also have to remeasure GCPs, i f required. 2.3.5 Summary A time series model with a recursive estimator i s preferred for the following reasons: (a) It i s superior to the polynomial approach i n terms of o v e r a l l geometric accuracy. (b) In the navigation process, fewer GCPs are required. This i s of great operational importance since high q u a l i t y GCPs are rare and d i f f i c u l t to obtain. (c) A recursive estimator processes one GCP at a time. A badly measured GCP i s immediately noticeable and can be rejected. (d) In the case of scanner imagery, the r o l l angle across the 28 scene is "reconstructed". This angle is of importance in the r e l i e f displacement correction phase. This thesis extends this technique to TIROS-N AVHRR (Section 4.3.2) which suffers from a much higher panoramic1 distortion, covers a larger area for a scene and has a considerably lower resolution than Landsat imagery. 2.4 Remapping and Interpolation 2.4.1 Background Remapping is the determination of a transformation, given by f, and b g f c in Equation 2.2, from the output grid (x,y) to the input grid (u,v). This transformation, when combined with interpolation, must give a fast computation of the intensity at (x,y) and store i t away without excessive disk t r a f f i c (both read and write). Remapping and interpolation are highly related. Therefore, when designing a remapping algorithm, one must take into account the eventual interpolation method. Existing methods are processing by blocks [ZOBR82] or by strips [RAMA77], In processing by blocks, the output grid is segmented into blocks and one block at a time is processed in core. The input area containing the output block is read into core while processing the block. In processing by strips a l l the input lines required to process an 1 Since pixel sampling in a scan line proceeds in equal time intervals, the ground segment imaged is proportional to the tangent of the scan angle. This effect, combined with the earth's curvature, produces a panoramic distortion in the along scan direction. 29 output l i n e are stored and processed i n core. The core size requirement i s then a function of output image s i z e and of the rotation angle between input and output. This section discusses a processing method i n which remapping and i n t e r p o l a t i o n are performed i n one-dimensional passes over the image. It o f f e r s the following advantages: (a) For each pass, either a l l the column numbers or a l l the row numbers remain the same between input and output. This means that one row/column of input data corresponds to the same row/column of output data. (b) One-dimensional.interpolation i s implemented i n each pass. (c) The algorithm can be performed in-place, hence reducing storage requirements. (d) R e l i e f displacement c o r r e c t i o n can be incorporated i n t o the f i r s t pass which operates i n the along scan d i r e c t i o n . 2.4.2 B i l i n e a r Remapping Fast transformation from the output image coordinates (x,y) to input image coordinates (u,v) can be achieved by segmenting the output image i n t o blocks. The transformation, for each block i s approximated by two piecewise b i l i n e a r functions of the output coordinates. The c o e f f i c i e n t s f o r each block are determined by f i r s t mapping the four corner points from (x,y) to (u,v). T h i s , c a l l e d inverse transformation, cannot be expressed i n closed form. Derivation of t h i s inverse transformation s t a r t s with the forward transformation from (u,v) to (x,y) which can be expressed i n 30 closed form and Is a function of: - sensor geometry, - s a t e l l i t e orbit ephemeris, - sa t e l l i t e attitude (not applicable in SAR imagery), - method of processing, and - earth spheroid and earth rotation. Once the transformation parameters are determined (from navigation), i t is possible to transform from image line and pixel coordinates to ECR coordinates and eventually to any desired output map projection, a l l in closed form. Piecewise bilinear functions are the simplest functions that guarantee continuity at the block boundaries. The block size is set by the required geometric accuracy. The bilinear functions are developed in Appendix D. Alternatively, an affine transformation can be used instead of a bilinear transformation. But an affine mapping is in general discontinuous from one block to the next. This can introduce image artifacts and consequently i s not discussed further here. While remapping by the bilinear transformation between (x,y) and (u,v) is straightforward, a direct usage of the transformation between these coordinates w i l l c a l l for two-dimensional interpolation. The remainder of Section 2.4 shows that a bilinear transformation can be 31 established between the input and output of each one-dimensional processing pass, and that the b i l i n e a r transformation i n each pass can be reduced to an equivalent one-dimensional l i n e a r transformation. 2.4.3 Two-Pass Processing Two-dimensional i n t e r p o l a t i o n i s slow and expensive to implement. A further complication arises when the i n t e r p o l a t i o n kernel increases in s i z e i n both d i r e c t i o n s . However, i t may be possible to modify the process into two one-dimensional i n t e r p o l a t i o n processes 1 which are easier to implement and more e f f i c i e n t with respect to i n t e r p o l a t i o n time requirement. Also, i t has the added advantage that r e l i e f displacement and other high frequency d i s t o r t i o n s can be corrected i n the f i r s t pass. Let p,q be the coordinates of an intermediate image as a r e s u l t of the f i r s t remapping and i n t e r p o l a t i o n pass. It i s po s s i b l e to s p l i t the b i l i n e a r transformation into two parts [TRW77] as follows: F i r s t pass: u = f ^ p . q ) (2.14a) v = q ' (2.14b) Second pass: p = x (2.15a) q = g^x.y) (2.15b) This i s possible i f the i n t e r p o l a t i o n kernel H(x,y) i s separable. That i s , H(x,y) can be decomposed into two functions H (x) and H^(y) where H(x,y) = H 1 ( x ) H_(y). For example, i f H(x,y) = smc(x) sinc(y) , then H(x,y) i s separable with H.(x)=sinc(x) andH-(y) = s i n c ( y ) . 32 The intermediate image and the output image are part i t i o n e d into blocks, and and g^ are b i l i n e a r polynomials. Since i n the value of q i s fixed and i n g^ the value of x i s f i x e d , f j and gj degenerate to l i n e a r functions. Geometrically i t means that a l l pix e l s are moved h o r i z o n t a l l y to t h e i r f i n a l x - p o s i t i o n f i r s t . Also i n t h i s pass, a l l along scan high frequency e r r o r s , such as r e l i e f displacement, sensor o f f s e t , l i n e length v a r i a t i o n , earth r o t a t i o n and p i x e l spacing n o n l i n e a r i t y can be corrected. Furthermore, i t has the s i g n i f i c a n t advantage that a two-dimensional i n t e r p o l a t i o n kernel reduces e f f e c t i v e l y to a one-dimensional kernel since the kernel weight i n the v - d i r e c t i o n i s zero at a l l values of v when v = q. In the second pass which i s performed along columns of data, the kernel again reduces to one-dimension. Lines of an image are stored i n consecutive records i n a t y p i c a l system. Therefore, to perform the second pass remapping and i n t e r p o l a t i o n , the intermediate image (p,q) should be rotated through 90° f i r s t . Fast matrix transpose methods f o r large f i l e s have been described by Eklundh [EKLU72] and Goldbogan [G0LD81], Note that matrix r o t a t i o n through 90° can be achieved by matrix transposition and reading the r e s u l t i n g rows of data backwards. . Figure 2.2 i l l u s t r a t e s the remapping and i n t e r p o l a t i o n d i r e c t i o n s i n each pass. 2.4.4 Three-Pass Processing The output image (x,y) can assume any or i e n t a t i o n which i s user 33 y - DIRECTION (SECOND PASS RESAMPLING DIRECTION) AND q- DIRECTIONS u- DIRECTION (FIRST PASS RESAMPLING DIRECTION ) x - DIRECTION • ORIGINAL SAMPLES O SAMPLES AFTER FIRST PASS A FINAL SAMPLES AFTER SECOND PASS 2.2 TWO ONE-DIMENSIONAL REMAPPING AND INTERPOLATION PROCESSES 34 specified. If there is a rotation with respect to the raw image (u,v), aliasing i s introduced in the two-pass processing method resulting in incorrect radiometric values. Aliasing i s explained as follows. In the two-dimensional method, the intensity at a A-marked point i s obtained by interpolating on the neighbouring D-marked points. However, in the two-pass method, the intensity at the same A-marked point is obtained by interpolating on the 0-marked points. Let l j and 1^ be the across scan spacing between the Q-marked points and 0-marked points respectively as shown in the figure. Since 1^ > l j , aliasing occurs as a result of interpolation. Friedmann has solved the aliasing problem in the case of image rotation by modifying the two-pass one-dimensional processing method into a three-pass one-dimensional processing method [FRIE81]. An oversampling procedure is introduced in the f i r s t pass (along scan) to push apart the spectrum replica of the original image in the frequency domain. The oversampling s is given by: where 6 is the angle of rotation and 0 < s < 1. The number of pixels per line in the f i r s t intermediate image is increased by a factor 1/s. The rotation matrix can be expressed as: s > 1 / (1 + tane) (2.16) (2.17) where (x',y') is the unrotated image coordinate system (the y'-axis aligns 35 with s a t e l l i t e v e l o c i t y ) and 6 i s measured i n the counter clockwise sense. In Friedmann's three-pass method, image r o t a t i o n can be performed properly with an oversampling operation i n the f i r s t pass, and the whole process can be written as: X -s 0-0 1. ' 1 0 If -s tan 0 sec 6 J L cos 6 /s s i n 6 /s"| 0 yj (2.18) v. M, where i s a row d i r e c t i o n (along scan) oversampling operation with s c a l i n g factor s along a row of data, i s a column d i r e c t i o n (across scan) operation and i s another row d i r e c t i o n operation. Operations M , M 2 and are performed i n order, r e s u l t i n g in an image a f t e r each operation. Figure 2.3 shows the three-pass r o t a t i o n method. This thesis extends Friedmann's idea to the remapping and i n t e r p o l a t i o n of s a t e l l i t e imagery i n which, i n addition to r o t a t i o n , along scan and across scan d i s t o r t i o n s have to be corrected. The whole process i s conceptually depicted i n Figure 2.A. At f i r s t , assume that there i s no r o t a t i o n and the two-pass processing given by Equations 2.14 and 2.15 (with x replaced by x' and y by y') w i l l produce an image free of a l i a s i n g . Then the three-pass r o t a t i o n i s performed. The two steps give a five-pass process as shown i n Figure 2.4a. Let the f i r s t pass operation be denoted by H and the second pass operation by H_ The oversampling operation can be coupled with r e s u l t i n g i n modifications to H , M„ and M_ which are now H ' , M 1 and M ' as shown i n SAMPLES AFTER THE FIRST PASS O AND n SECOND PASS * THIRD PASS + NOTE: O ARE THE OVERSAMPLED'PIXELS FIGURE 2.3 THREE ONE-DIMENSIONAL REMAPPING AND INTERPOLATION PROCESSES 37 FIGURE 2.4 COUPLING OF PASSES 38 Figure 2.4b. Since R^ ' and are neighbouring column operations, they can be coupled to form one operation. The f i n a l result is a three-pass process as shown in Figure 2.4c. 2.5 Interpolation Kernel Interpolation is the computation of the intensity of each output pixel with an interpolation kernel. This section examines different kernels and proposes one which is most suited to the three-pass process presented in Section 2.4. 2.5.1 Outline of Interpolation Methods 2.5.1.1 Nearest Neighbour and Linear Interpolations Nearest neighbour interpolation is the simplest form of a l l interpolation methods. The interpolation kernel is given by: 1 -1/2 < x < 1/2 (2.19) o elsewhere Its frequency response is a sine function given by: F n ( f ) = sinc(TTf) (2.20) where sinc(x) is defined as sin(x)/x. Another simple form is linear interpolation given by: 1 - x x < 1 f x ( x ) = < (2.21) 0 elsewhere 39 Its frequency response is given by: Fj(f) = [ F Q 2 ( f ) ] = sinc 2(TTf) (2.22) The interpolation kernels and their frequency responses are plotted in Figure 2.5. These two interpolation kernels perform poorly in terms of radiometric accuracy, as shown by their frequency response deviations from the ideal low pass f i l t e r , when applied to image interpolation. But they are adequate for DTM interpolation as w i l l be discussed in Chapter 3. 2.5.1.2 Cubic Spline A spline i s a smooth curve passing through a set of discrete sample points which are not necessarily uniformly spaced. A common spline is the cubic spline which is defined as follows. Let the sample points be t^ where i is the sample point index. Then i n the interval t^ to t^ +^ » a cubic polynomial i s fitted and the f i t t i n g function, i t s f i r s t and second derivatives are continuous at both end points. The interpolation in this interval is given by: (t-t ) 2 (t-t ) 3 f ( t ) = C Q. + C, .(t-t.) + C 2 1 — ^ + C 3.—- 6 (2.23) From the continuity conditions, the four coefficients C^ to C^ .. can; be found [DEB079] and the result is presented below without proof here. Denote the data points by F^ , the slope of f(t) at t^ by s^, the spacing between t^ and t 1 + 1 by I K ; i.e., h. = t... - t. (2.24) l l+l I (a) I n t e r p o l a t i o n K e r n e l s Magnitude (b) Frequency 0.5 1.0 1.5 2.0 Frequency i n C y c l e s p e r P i x e l N e a r e s t Neighbour I n t e r p o l a t i o n L i n e a r I n t e r p o l a t i o n FIGURE 2. 5 NEAREST NEIGHBOUR AND LINEAR INTERPOLATION - KERNELS AND FREQUENCY RESPONSES 41 and the forward diffe r e n c e by F [ t . , t.,,] where 3 1 l+l nt±i t . + 1 ] = ( F. + 1 - F.) / h. then: C n. = F. Oi l C. . = s . I i l C 2. = 2 {3 F [ V t 1 + 1 ] - 2 s. - B 1 + 1-} / h. C_. = 6 {-2 F [ t . , r . ] + s. + s.^.} / h." 3i I l+l I l+l I (2.25) (2.26a) (2.26b) (2.26c) (2.26d) It remains to determine the slopes at a l l data points. This can be done by solving the set of l i n e a r equations: h. = s. , + 2 (h. + h. .) s. + h. . s... = b. i i-1 l l - l l l - l l+l l where b i - 3 { h i F [ t i - 1 « C i ] + V l F t t i ' ^ + 1 ^ (2.27) (2.28) and 2 < i < n-1 where n i s the t o t a l number of points. In matrix form, t h i s i s : " s l " r b 2 -1 • • • • • • s _ n _ _ V 2 . (n-2) x n nxl (n-2)xl (2.29) It requires 0(n) operations to solve for the s_.'s , Note that there are n-2 equations with n unknowns. It i s necessary to choose two of the n 42 values of s i . This can be done by taking s^ as the slope at t^ of the quadratic polynomial which agrees with F at » t^ and t ^ + j • Then: h i F [ V l > t i ] + h i F [ V V l 1 s = r—T-T (2.30) h i + h i + l The error in interpolation is given by: I e | < (5/384) h ± A | f ( 4 ) | (2.31) (4) where f is the fourth derivative of f ( t ) . 2.5.1.3 B-Spline Another class of splines i s the set of B-splines. A B-spline i s a polynomial f i t i n each i n t e r v a l t^ to t ^, and the 0-th to ( n - l ) t h d e r i v a t i v e s are continuous where n i s the polynomial order. Interpolation i s performed with the aid of a set of basis function B: n f ( t ) = E a. B. (t) (2.32) i=l 1 1 , k -where n i s the t o t a l number of points, k-1 i s the polynomial order, a i s a set of c o e f f i c i e n t s to be determined, and B. (t) belongs to the set B. The set B i s defined r e c u r s i v e l y : ( 1 , t < x < t B (t) =( (2.33a) ' \ 0 , elsewhere B i k ( t ) • t ~ ~ — = V B i k - i ( t ) + [ 1 + k - \ B 1 + i k - i ( t ) ( 2 - 3 3 b ) 1 , l c c i + k - i t i 1 » k 1 h+k h-i ' 1 with k > 1 and 0/0 = 1. The evaluation of the B-spline by the 43 recurrence r e l a t i o n i s very stable since only additions of non*-negative qu a n t i t i e s are involved. For a given t say t j > t > t > then only the following B-splines are non-zero: B - i V l , 2 B j , 2 Bj-2,3 Bj-1,3 B j , 3 It remains to determine a . This i s done by solving a set of simultaneous l i n e a r equations: J , " i B i , k ( t ) = F i (2.34) In matrix form t h i s i s : rF,"| 1 l • = • • a F n n (2.35) The matrix 0 i s a function of B. and i s diagonally dominant. 1 , K 2.5.1.4 The Sampling Theorem Let F^ be an i n f i n i t e s t r i n g of pulses sampled on a function f ( t ) at equal time i n t e r v a l s iT where T i s the time spacing between samples. The o r i g i n a l function f ( t ) and i t s spectrum can be recovered completely i f the following conditions are s a t i s f i e d : 44 (a) The s i g n a l f ( t ) i s band limited with bandwidth W; i . e . , i t s highest s i g n i f i c a n t frequency i s less than W cycles/second. (b) The s i g n a l i s sampled at more than the Nyquist rate of 2 W samples per second. Th i s i s known as the sampling theorem [LATH65, JERR77 and SAKR68]. For a given s i g n a l , the highest frequency 1 that can be represented i s 0.5 c y c l e / p i x e l and the i n t e r p o l a t i o n i s given by: f ( t ) = E F. sine [ T r ( t - i ) ] (2.36) i = -co 2.5.1.5 Cubic Convolution The i n t e r p o l a t i o n given by Equation 2.36 requires an i n f i n i t e number of data points. Cubic convolution uses a truncated sine function which extends over four points. The kernel i s modified so that the f i r s t d e r i v a t i v e , as a result of i n t e r p o l a t i o n , i s continuous. The kernel i s [BERN76] : (2.37) 1 - 2 |x| + |x| 3 , 0 < |x| < 1 f(x) = { A - 8 |x| + 5 |x| 2 - |x| 3 , 1 < |x| < 2 0 , elsewhere Cubic convolution i s advocated by Bernstein [BERN76] and Simon [SIM075b] because the kernel s i z e i n two-dimensional i n t e r p o l a t i o n i s only 4X4. 1 For ease i n discussion, frequency here i s expressed as c y c l e s / p i x e l . This can e a s i l y be converted to cycles/second i f the sampling rate i s known. 45 Shlien [SHLI79] has also investigated other splines with different kernel sizes, which are approximations to the ideal sine function. 2.5.1.5 Comparison of Interpolation Kernels The spline methods are not pursued further here. The more traditional engineering approach of using the convolution of sample points with the sine function (Equation 2.36) is used in the remainder of the thesis. Figure 2.6 shows the frequency response of the cubic convolution kernel. The ideal low pass f i l t e r , which the sine function yields, is inserted in the figure for comparison purposes. Interpolation experiments have been performed with sine waves of different frequencies and the RMS interpolation error i s also shown in the figure. Since the frequency response crosses 1 at 0.28 cycle per pixel, a dip occurs in the RMS error around this frequency. Beyond this frequency, the response deviates drastically from the ideal low pass f i l t e r , and this results in a rapid growth in the RMS error. With frequency below 0.28 cycle per pixel, the RMS error ranges from 0 to 7 levels out of 256 levels (assuming 8-bit image data). Note that RMS error and frequency are related as follows. Let R(f) be the frequency response for frequency f. Then the response of a sine wave with magnitude of 1 unit w i l l have a magnitude of R(f) units after interpolation. The RMS error for frequency f is measured between the RMS ERROR 12.0-f FREQUENCY RESPONSE 10.0. 8.0-6.0-4.0. 2.0 0.0 0.1 0.2 FREQUENCY IN CYCLE PER PIXEL 0.3 -0.98 \ Note: The RMS error includes the effect of tabulating the interpolation kernel at every 1/32 data point interval, See Section 2.5.2.1. FIGURE 2.6 RMS ERROR AND FREQUENCY RESPONSE FOR CUBIC CONVOLUTION 47 output and input sine waves, which have the same frequency and phase angle but different amplitudes. The RMS error can be reduced by expanding the f i l t e r size. Since the three-pass process only calls for one-dimensional interpolation, the number of arithmetic operations increases only linearly with kernel size. Figure 2.7 shows the frequency response of a 16-point sine function (see Section 2.5.2.2) modified by a Kaiser window (see Section 2.5.2.3) and the RMS error as a result of interpolation on sine waves of different frequencies. A comparison with Figure 2.6 shows that this 16-point sine function gives a much lower RMS interpolation error. In particular, the MSS and TM modulation transfer functions have significant frequency responses up to 0.16 cycle per pixel. Below this frequency, the overall RMS interpolation error is about 1 level. Based on these results, this kernel is used in the interpolation phase. 2.5.2 Preprocessing Preprocessing of the sine function consists of the following three steps: 2.5.2.1 Tabulation of Sine Function For an interpolation of f ( t ) , i t i s necessary to compute the value of the sine function at the data points. Such computation can be avoided i f the sine function is replaced by a table and the f i l t e r weight can be obtained by a simple table look up technique. This requires the sine 48 Note: The RMS error includes the effect of tabulating the interpolation kernel at every 1/32 data point interval. See Section 2.5.2.1. FIGURE 2.7 RMS ERROR AND FREQUENCY RESPONSE FOR SINC FUNCTION INTERPOLATION 49 function to be divided into s t r i p s and a step function i s used i n each s t r i p to approximate the sine function as shown in Figure 2.8a. Therefore, preprocessing i s required to compute the table. The i n t e r p o l a t i o n accuracy i s then a function of the s t r i p s i z e . The worst case i s that each F l i e s at the edge of a s t r i p . The result i s that the entire function i s s h i f t e d by 1/2 s t r i p as shown i n Figure 2.8b. The s t r i p width can be determined as follows. The modulation t r a n s f e r function (MTF) of a t y p i c a l mechanical scanner dictates that i t takes about three p i x e l s to respond from one i n t e n s i t y extreme to the other ( t h i s corresponds to about 0.16 cycle per p i x e l ) . For th i s worst case on 8-bit imagery, the error i s : E r r o r = 0.5 x s t r i p width x 256 / 3 l e v e l s , where 256/3 i s the slope of the response i n lev e l s per p i x e l . Therefore: E r r o r = s t r i p width x 43 l e v e l s . Assuming uniform d i s t r i b u t i o n i n i n t e n s i t y slope, the RMS error i s then: RMS error = s t r i p width x 43 / 3 l e v e l s . I f the s t r i p width i s 1/32 p i x e l , the RMS error i s then below one l e v e l which i s acceptable i n many a p p l i c a t i o n s . The MTFs of SPOT's l i n e a r array sensors are not known at present, but the s t r i p width can be determined i n the same manner once the MTFs become known. s \ i n c f u n c t i o n / S t e p f u n c t i o n (a) Approximation FIGURE 2. 8 APPROXIMATION OF SINC FUNCTION BY STEP FUNCTION 51 2.5.2.2 Kernel S i z e Determination It takes an i n f i n i t e number of points to synthesize the i d e a l low pass f i l t e r . T h i s i s re f l e c t e d i n Equation 2.36 i n which the summation i s over a l l data points. However, the maximum magnitude of each s i d e l o b e 1 of h(t) decreases as a function of distance away from the centre. For example, i n the 8-th sidelobe the maximum magnitude i s about 1/27 of the magnitude at the centre. This suggests that the i d e a l but i n f i n i t e sine function can be approximated by a truncated sine f u n c t i o n . Any approximation would introduce error i n h(t) i n the frequency domain where the spectrum of the truncated sine function i s no longer an i d e a l low pass f i l t e r . The most important tradeoff i n performing image i n t e r p o l a t i o n i s speed versus accuracy. The larger the f i l t e r , the greater the i n t e r p o l a t i o n accuracy and the slower the computation speed. The accuracy obtainable with any given kernel si z e can be computed i n a s t a t i s t i c a l sense [SHLI79]. Let the function F(t) have autocorrelation R: E [F(t) F(t+/r)] = R(T) (2.38) where E denotes expected value. Let f ( t ) be an estimate of F(t) given by: f ( t ) = Z F. h ( t - i ) (2.39) i where h(t) i s the sine function multiplied by a window. Then the RMS Sidelobes i n the time domain must not be confused with sidelobes i n the frequency domain. It should be clear which domain i s referred to i n the t ext. 52 error i n i n t e r p o l a t i o n i s given by: e 2 ( t ) = E{[F(t) - f ( t ) ] 2 } (2.40) which a f t e r expanding gives: ( t - i ) + e 2 ( t ) = R(0) - 2 Z R ( t - i ) h + ZZ R( i - j ) h ( t - i ) h ( t - j ) (2.41) Thus the RMS error can be calculated f o r any kernel s i z e provided the autocorrelation function R of the data i s known. The autocorrelation function i s scene and t e r r a i n type dependent. A Seasat scene (Orbit 230) and Landsat scene (I.D.2921-18025) are used to study t h i s function. Both scenes were imaged over the Vancouver area and cover d i f f e r e n t t e r r a i n types such as mountain, sea and c i t y . Figure 2.9 i l l u s t r a t e s the autocorrelation function for a few t e r r a i n types. Using the function f o r various t e r r a i n types, the RMS i n t e r p o l a t i o n error versus kernel s i z e i s derived and i s shown i n Figure 2.10. From t h i s f i g u r e i t i s seen that f o r a kernel s i z e of 16, the maximum RMS error i s less than 0.5 l e v e l out of 256 f o r both scenes. The kernel s i z e used i n the remaining part of the thesis i s 16 unless otherwise stated. 2.5.2.3 Windowing The sine function i n t e r p o l a t i o n kernel has alternate algebraic signs for successive weights i n the main lobe. T h i s , together with truncation i n the sine function, produces a " r i n g i n g " e f f e c t , i . e . , Gibb's phenomenon, for sharp bends i n the data values. In the frequency domain, the truncated sine function no longer corresponds to an i d e a l low FIGURE 2.9 EXAMPLES OF AUTOCORRELATION FUNCTION FIGURE 2.10 RMS INTERPOLATION ERROR VERSUS KERNEL SIZE .55 pass f i l t e r , but a f i l t e r with sidelobes. A windowing function can,be used to minimize the r i p p l e s i n these sidelobes, at the expense of narrowing the main lobe. This w i l l reduce ringing i n the time domain. One example i s the Kaiser window given by: where N i s the t o t a l number of points i n the sine function, J Q ( B ) i s the zeroth order Bessel function [K0RN61] and g controls the tradeoff between sidelobe peaks and width of the main lobe. The i n t e r p o l a t i o n kernel i s multiplied by the Kaiser window p r i o r to the i n t e r p o l a t i o n process. Figures 2.11 1 and 2.12 are plo t s of the 4-point truncated sine function without and with Kaiser windowing. The ef f e c t on the sidelobes i s immediately obvious by comparing these two f i g u r e s . Figure 2.13 shows the frequency response of the 16-point sine function with Kaiser windowing. Another example i s a Dolph-Chebyshev window which makes the sidelobe peaks a l l equal [RABI75]. When the sine function i s divided into s t r i p s , the m u l t i p l i c a t i o n by a window imposes very l i t t l e overhead i n the preprocessing. 2.5.3 Nonuniformly Spaced Samples The above i n t e r p o l a t i o n kernels are only applicable to samples 1 Figures 2.11 to 2.13 were provided by Mr. D. Friedmann at MDA. J, 0 w(t) = (2.42) 56 0.25 0.50 0.75 1.00 1.25 1.50 CYCLE/PIXEL FIGURE 2.11 FOUR-POINT SINC FUNCTION AND ITS FREQUENCY RESPONSE CYCLE/PIXEL FIGURE 2.12 FOUR-POINT SINC FUNCTION WITH KAISER WINDOW AND ITS FREQUENCY RESPONSE , • • u l_f 1 -f, I — I 0.25 0.50 0.75 1.00 1.25 1.50 CYCLE/PIXEL FIGURE 2.13 SIXTEEN-POINT SINC FUNCTION WITH KAISER WINDOW AND ITS FREQUENCY RESPONSE 59 which are uniformly spaced. However, for Landsat-4 TM imagery, the data samples are no longer equally spaced i n the across scan d i r e c t i o n due to the presence of j i t t e r , v a r i a t i o n s i n s a t e l l i t e a l t i t u d e and ground speed, earth's obliqueness and bowtie e f f e c t [GSFC82]. A l l these d i s t o r t i o n s w i l l create "scan gaps" which r e s u l t i n deviations from nominal sample spacing. The sine kernel i s no longer s u i t a b l e f or the across scan i n t e r p o l a t i o n . YEN [YEN56] has shown that the s i g n a l can s t i l l be reconstructed from the d i s c r e t e samples. In the case of no scan gap, the kernel he proposes degenerates to the sine function. A modification of Yen's method i s applied to the present problem. Since t h i s i s s p e c i f i c only to Landsat-4 TM imagery i n the across scan pass, detailed discussion on t h i s method i s deferred to Chapter 3. 2.6 System Data Flow Figure 2.14 depicts o v e r a l l system data flow. It i s divided into three parts: (a) Processing of a u x i l i a r y data so that the useful data are stored i n disk f i l e s . (b) Generation of the three-pass remapping and i n t e r p o l a t i o n c o e f f i c i e n t s . (c) Use of remapping c o e f f i c i e n t s and the DTM. REMAPPING AND INTERPOLATION AUXILIARY DATA PROCESSING REMAPPING AND INTERPOLATION COEFFICIENTS AND INTERPOLATED. DTM EARTH PARAMETERS GENERATE 3-PASS REMAPPING AND| INTERPOLATION COEFFICIENTS 3-PASS REMAPPING & INTERPOLATION COEFFICIENTS DTM A.D. = AUXILIARY DATA FIRST PASS REMAPPING & INTERPOLATION ROTATE 90° SECOND PASS REMAPPING & INTERPOLATION ROTATE 90< THIRD PASS REMAPPING & INTERPOLATION FIGURE 2.14 SYSTEM DATA FLOW o 61 Generation and usage of remapping and interpolation coefficients have been discussed above. This section discusses auxiliary data processing and the relief displacement correction in the f i r s t pass of remapping and interpolation. 2.6.1 Auxiliary Data Processing Auxiliary data contains various kinds of correction information that the user can apply to produce a geometrically and radiometrically corrected image. The contents of auxiliary data can be divided into three areas: - s a t e l l i t e data for geometric correction, - sensor data for geometric and radiometric correction, and - earth and projection dependent data for placing the image on a map grid. Sa t e l l i t e data contains orbit information and can be in the form of orbital elements (which are parameters to describe the orbit) and/or sampled sa t e l l i t e position and velocity. Orbit information is used to determine the s a t e l l i t e position as a function of time. Satellite data also contains attitude (pitch, r o l l and yaw) information. Sensor geometric data contains information which, when combined with s a t e l l i t e data, can be used to determine the geolocation of any pixel in the image. Both sets of data are obtained from telemetry design specifications or telex from the appropriate agencies. These data have to be processed so that the orbit and attitude information are stored in 62 data f i l e s which w i l l be used i n the generation of remapping c o e f f i c i e n t s . In t h i s way the system has the f l e x i b i l i t y to handle imagery from various s a t e l l i t e s since the transformation i s table driven and has become s a t e l l i t e independent. Remapping accuracy depends upon the accuracies of the orbit and a t t i t u d e information. These data can usually be refined by ground control points (GCPs). The transformation accuracy w i l l have an impact on the eventual i n t e r p o l a t i o n of the DTM. How t h e i r accuracies a f f e c t the use of the DTM i s discussed i n Chapter 3. Earth parameters are used to determine where an instantaneous view vector i n t e r s e c t s the earth spheroid. Projection parameters contain information f or mapping a point on the earth spheroid to a two-dimensional surface. The map pro j e c t i o n parameters depend on the chosen map p r o j e c t i o n . Earth parameters also include t e r r a i n f l u c t u a t i o n s . A DTM must be available i n order to correct f o r r e l i e f displacement. 2.6.2 R e l i e f Displacement Correction The usual geometric errors i n s a t e l l i t e borne imagery are slow varying and thus the output image can be segmented into blocks as has already been argued. R e l i e f displacements are fast varying and, therefore, i t i s e s s e n t i a l to separate r e l i e f displacement correction from the b i l i n e a r model. To correct f o r r e l i e f displacement, r e l i e f information must be av a i l a b l e i n the form of a DTM. 63 A simple, f a s t and operationally f e a s i b l e . c o r r e c t i o n f or r e l i e f displacement can be performed on the raw image data. Here, the d i s t o r t i o n i s only in the along scan d i r e c t i o n , because i t i s in t h i s d i r e c t i o n only that the displacement due to c e n t r a l perspective occurs. The b i l i n e a r model for geometric c o r r e c t i o n can be implemented i n three passes over the same image as shown i n Section 2.4. The f i r s t pass i s i n the along scan remapping and i n t e r p o l a t i o n process which i s i l l u s t r a t e d i n Figure 2.15. R e l i e f displacement c o r r e c t i o n i s incorporated into t h i s pass. Other fast varying errors such as l i n e length v a r i a t i o n , sensor o f f s e t and earth r o t a t i o n are also incorporated i n t o the f i r s t pass. As shown i n the f i g u r e , the raw image (Box A) i s read one l i n e (Box B) at a time into core. It i s desired to correct f or a l l the along scan d i s t o r t i o n s and hence map t h i s l i n e into an output buffer (Box C). The mapping consists of both low frequency and high frequency components. The low frequency component i s modelled by the along scan c o e f f i c i e n t s (Box D). The high frequency components, r e l i e f displacement (Box E), l i n e length v a r i a t i o n and sensor o f f s e t (Box F) depend on p i x e l l o c a t i o n . The DTM i s used to correct f or r e l i e f displacement i n the following manner. To prevent a two-dimensional i n t e r p o l a t i o n on the DTM, the i d e a l case would be to i n t e r p o l a t e i t as an o f f - l i n e process to , the image coordinates defined a f t e r the f i r s t pass. In t h i s way the t e r r a i n height of the point being processed can be read e a s i l y and the r e l i e f displacement can then be computed. However, t h i s i s operationally unwise since the interpolated DTM i s only s p e c i f i c to one p a r t i c u l a r scene. Furthermore, i n the case of TM imagery where scan gaps occur, i t 64 RAW IMAGE LINE OF IMAGE D PIXEL SHIFT I ALONG SCAN RESAMPLING DTM PROJECTION COEFFICIENTS IMAGE LINE AFTER THE FIRST PASS E CONVERSION TO DISPLACEMENT H HEIGHT F DTM SCROLLING SENSOR OFFSET LINE LENGTH ORRECTION ETC BUFFER RESAMPLED DTM FIGURE 2.15 FIRST REMAPPING AND INTERPOLATION PASS 65 i s d i f f i c u l t to get a DTM to such a p r o j e c t i o n . A more f e a s i b l e way i s to i n t e r p o l a t e the DTM to a map p r o j e c t i o n close to the input image coordinates (Box G); e.g., assume nominal orbit parameters and zero a t t i t u d e angles. Then a s c r o l l i n g buffer (Box H) i s used i n main memory during the along scan remapping and i n t e r p o l a t i o n process. The buffer should contain enough l i n e s so that the heights of a l l the p i x e l s of the current image l i n e can be found i n the buffer. Indexing to the buffer i s done by predetermined c o e f f i c i e n t s (Box I) which vary from scene to scene. T h i s i s operationally f e a s i b l e since the orientations of a l l scenes over the same area do not d i f f e r s i g n i f i c a n t l y . Hence, the interpolated DTM can be used i n a l l these scenes and only the indexing c o e f f i c i e n t s have to be changed. Inte r p o l a t i o n of the DTM can also be done by the three-pass processing method since the DTM can be treated as an image with the height values representing the p i x e l i n t e n s i t i e s . The remapping c o e f f i c i e n t s can r e a d i l y be determined since each point i n the DTM represents a c e r t a i n geographic location and each l i n e , p i x e l coordinate i n the "nominal" image can be mapped into a geographic l o c a t i o n . Therefore, i t i s possible to develop a transformation between the DTM and the raw "nominal" image coordinates. The s i z e of the DTM s c r o l l i n g buffer i s then a function of the yaw angle. 2.7 Chapter Summary A u n i f i e d approach to the r e c t i f i c a t i o n of remotely sensed imagery has been developed. The key features are: Geometric transformation parameters, e.g., s a t e l l i t e attitude angles, are refined i n the navigation process from a p r i o r i estimates using a recursive estimator and GCPs. Only a reasonable number of GCPs per scene are required to achieve a high geometric accuracy. Remapping and i n t e r p o l a t i o n i s performed with three one-dimensional passes over the image. The- one-dimensional processing i s the key to y i e l d a high throughput. Remapping i s performed v i a a pair of piecewise b i l i n e a r polynomials and fast computation can be achieved. A reasonable i n t e r p o l a t i o n kernel size (say 16) i s used to r e t a i n radiometric accuracy. The kernel i s tabulated, again f o r speed i n computation. R e l i e f displacement and other along scan d i s t o r t i o n s can be corrected i n the f i r s t pass over the image. 67 CHAPTER 3 DIGITAL TERRAIN MODELS 3.1 Introduction In order to correct remote sensing data f o r r e l i e f displacement, ele v a t i o n information must be av a i l a b l e i n d i g i t a l form. A DTM i s an ordered array of numbers that represents the s p a t i a l d i s t r i b u t i o n of t e r r a i n c h a r a c t e r i s t i c s . The usual c h a r a c t e r i s t i c i s the t e r r a i n height at the given x,y p o s i t i o n (when t h i s i s the case, the DTM i s also referred to as DEM or D i g i t a l E l e v a t i o n Model). The x,y coordinate system i s usually i n some map projection or i n lat i t u d e / l o n g i t u d e . T h i s chapter describes the types of DTMs, t h e i r a c q u i s i t i o n methods, s p a t i a l accuracies and i n t e r p o l a t i o n s . F i n a l l y t h e i r accuracy requirements and i n t e r p o l a t i o n methods with respect to the r e c t i f i c a t i o n process are inve s t i g a t e d . 3.2 A c q u i s i t i o n of DTMs DTM data may be acquired by several methods [DOYL78]. The usual and least expensive method i s to manually d i g i t i z e e x i s t i n g mapsheets. This i s done on a commercially available d i g i t i z i n g t a b l e . However, t h i s method i s very laborious and demands a l o t of patience by the operator and i s thus susceptible to human error. An improved method i s to use automatic contour l i n e following instruments. In t h i s case, a contour map must be prepared by removing a l l contour labels and f i l l i n g i n the gaps. Th i s i s d e f i n i t e l y an improvement over the manual method except 68 f that the instruments may have d i f f i c u l t y with l i n e s of d i f f e r e n t weights such as indexed contours. Furthermore, there i s also d i f f i c u l t y i n automatically assigning elevations to each contour. Another automatic DTM a c q u i s i t i o n method i s to use a scanning device on e x i s t i n g maps. Depending on the type of scanner used, one l i n e or several l i n e s can be scanned at a time. The contour sheet i s placed on a drive which either rotates under a fixed l i g h t source or has the l i g h t source r o t a t i n g with respect to the d r i v e . Each time the scanner crosses a contour, the x and y coordinates are recorded. E l e v a t i o n assignment i s s t i l l a problem i n this method. DTMs acquired by the above approaches may leave large empty spaces between contours. To render these DTMs use f u l , time consuming i n t e r p o l a t i o n has to be applied to f i l l i n the spaces. The t y p i c a l accuracy of data from topographic maps i s : Scale C l a s s 1 S p a t i a l Accuracy E l e v a t i o n Accuracy a R90 a 1:50000 Al 25 m 16 m 10 m 6 m 1 :50000 B2 50 m 31 m 20 m 13 m 1:250000 Al 125 m 78 m 50 m 31 m 1:250000 B2 250 m 156 m 100 m 63 m where a i s the standard deviation and i s the value which 90% of the map f a l l s within. To obtain better DTM accuracies, a photogr ammetric stereo model can be used. 1 NATO standards, also adopted i n Canada. 2 Hi 1.6 o assuming Gaussian d i s t r i b u t i o n . 69 An example of a photogrammetric stereo model Is the Gestalt Photomapper described by K e l l y et a l [KELL77]. It i s a highly automated photogrammetric system designed to produce DTMs by using e l e c t r o n i c image c o r r e l a t i o n to measure parallax between the l e f t and the right images. The measured parallax can be transformed into height. The spacing between the rows and columns i n the DTM i s 182 urn on the o r i g i n a l photographs. On t y p i c a l 1:50,000 a e r i a l photographs, the p i x e l spacing corresponds to 9 metres. From these photographs, the r e s u l t i n g RMS errors i n the X, Y and Z axes i n the DTM are 0.6, 0.85 and 1 metre r e s p e c t i v e l y . The errors are even smaller i n the newer version of the Gestalt Photomapper (GPM I I - l ) . Another possible source of DTMs i s data obtained from altimeters on board spacecraft. Examples include the altimeter on board Seasat. i At the present time t e r r a i n data i s very sparse. The U.S. Defence Mapping Agency has d i g i t i z e d contour data from 1:250,000 maps for the whole U.S., and USGS plans to do the same from 1:24,000 maps. Other mapping organizations of Canada, U.K. and A u s t r a l i a are producing some DTMs. In the future i t i s expected that a world data base with good accuracy w i l l be b u i l t up. However, i t has been estimated that t h i s would require 2 x 10*~* b i t s of data to cover every square metre on land where each elevation i s d i g i t i z e d to 16 b i t s . 3.3 Types and In t e r p o l a t i o n of DTMs ( D i f f e r e n t types of DTMs may be distinguished by t h e i r coordinate 70 spacings. Since DTM i n t e r p o l a t i o n i s required (as an o f f - l i n e process) i n the r e c t i f i c a t i o n algorithms, t h i s section presents d i f f e r e n t DTM 1 types and t h e i r i n t e r p o l a t i o n methods. For a detai l e d d e s c r i p t i o n , r e f e r to [SCHU76b, ALLA78, PEUC78 andHEIL78]. 3.3.1 Irregular DTMs In an i r r e g u l a r DTM, reference points are not equally and regularly spaced. In some cases they are selected at random p o s i t i o n s , while i n other cases they are selected at positions where there i s a change i n slope or some d i s t i n c t feature. The l a t t e r method gives an adequate d e s c r i p t i o n of the topography with the smallest number of reference points. The contour DTM, acquired by d i g i t i z i n g contour maps, gives equally spaced elevation points. However, the s p a t i a l d i s t r i b u t i o n of reference points i s not regular, but i s selected at equal i n t e r v a l s along the contour or spaced more c l o s e l y along contours i n areas where the t e r r a i n i s h i l l y than when i t i s f l a t . P r i o r to using an i r r e g u l a r DTM, the general p r a c t i c e i s to transform i t to a regular DTM; i.e.,- to compute t e r r a i n heights at the nodes of a square g r i d . One method [SCHU76] i s to define the t e r r a i n height by: 2 2 h = + a^x + a^y + b^x + b 2xy + b^y (3.1) For the i n t e r p o l a t i o n of a point, the coordinates and the height of each of the surrounding reference points are substituted into t h i s 71 equation. Also each reference height i s given a weight that i s a monotonically decreasing function of the distance to the reference point. Then the parameters i n the equation are solved by a method of least squares. Taking the o r i g i n of the coordinate system to be the interpolated point, where i t s height i s required, then only a^ needs to be computed. To keep the computation time within reasonable bounds, only reference points within a s p e c i f i e d maximum distance from an interpolated point w i l l be used. Generally, going from one interpolated point to an adjacent one, the surface defined by Equation 3.1 w i l l change i t s o r i e n t a t i o n and, possibly, i t s shape. For t h i s reason, t h i s method i s also referred to as the "moving surface" method. Another method [HARD71] i s to construct around each reference point a fix e d surface which has a v e r t i c a l axis of symmetry at that point. I n t e r p o l a t i o n i s performed by summing the heights of a l l surfaces computed at that point. 3.3.2 Regular DTMs In a regular DTM, spacings of the s p a t i a l coordinates are regular. Hence, reference points are at the nodes of a regular g r i d . Since the regular DTM consists of an e a s i l y addressable array of elevation values, i t i s simpler than the i r r e g u l a r DTM to store and access i n computer storage. According to the sampling theorem, t e r r a i n v a r i a t i o n s which correspond to a wave length less than twice the grid spacing are not represented. 72 If the elevation of a point other than the nodes i s required, i t can be obtained by the i n t e r p o l a t i o n methods discussed i n Section 2.5, such as nearest neighbour, b i l i n e a r , cubic convolution or convolution with a l a r g e r kernel. If the DTM i s dense enough to represent a l l s i g n i f i c a n t changes i n slope, then b i l i n e a r or even nearest neighbour i n t e r p o l a t i o n i s adequate. 3.3.3 Triangulated Irregular Network Usually, t e r r a i n i s not regular but changes from one land area to another. If a regular DTM i s used to represent the t e r r a i n , i t has to be adjusted to be adequate f o r the roughest t e r r a i n . But then the data i s highly redundant i n smooth t e r r a i n . A triangulated i r r e g u l a r network [PEUC78 and HEIL78] t r i e s to represent t e r r a i n adequately with the smallest number of reference points. T h i s i s done by covering the t e r r a i n with a network of t r i a n g l e s . Each t r i a n g l e i s a " f a c e t t e d " representation of the t e r r a i n surface. The c o n t i n u i t y of the surface representation i s maintained by r e s t r i c t i n g the construction such that each t r i a n g l e has only one other t r i a n g l e sharing a s i d e . Vertices are chosen at l o c a l maximum or minimum points. The sides of the constructed t r i a n g l e s follow ridges and channels. Interpolation of a point i n s i d e a t r i a n g l e follows from the condition that the point l i e s on the three-dimensional t r i a n g l e formed by i t s surrounding v e r t i c e s : 73 Det X y z 1 x l y i Z l 1 X2 y2 Z2 1 x 3 y3 Z3 1 = 0 where (x^, y.^, z,), i = 1 to 3, are the v e r t i c e s of the t r i a n g l e . 3.4 DTM Accuracy 3.4.1 Accuracy Requirement i n Interpolated Coordinates The DTM i s warped as an o f f - l i n e process close to the input image coordinates (u,v). A DTM has three coordinates: x, y and z. Inaccuracies i n x, y and z w i l l cause an error i n the i n t e r p o l a t i o n p o s i t i o n i n the input image. The acceptable error i s less than 1/2 p i x e l . This subsection deals with the accuracy requirements i n x, y and z for the interpolated DTM, as a function of sensor and s a t e l l i t e parameters. The displacement d on ground i s related to t e r r a i n height h by: d = fh (3.2) where the f a c t o r f i s s a t e l l i t e and sensor s p e c i f i c and i s also a function of p i x e l p o s i t i o n . Then: Ad = f Ah (3.3) Let the p i x e l grid spacing be r . The usual requirement i s that d < 1/2 r . Hence: f Ah < r/2 and 74 Ah < r/(2f) (3.4) I 2 2 But the contributions to Ah are: spatial accuracy Ap = /Ax + Ay , and elevation accuracy Az. Let 6 be the slope. Then Ah can be expressed in the RMS sense by: Ah = ^Ap 2 t a n 2 e + A z ^ (3.5) Substituting into Equation 3.4, we have / a p 2 tan 26 + Az 2 < r/(2f) (3.6) To obtain the most stringent requirement, we have to allow for the worst possible case of rel i e f displacement; i.e., maximum f. The maximum value of f and the, value of r for various sensors are given in Table 3.1. Plots of the accuracy requirement for various sensors Figures 3.1 to 3.4 for terrain slopes of 20°, 30° and decrease in slope, the accuracy Ax (and hence Ay) can anti cipated. In Figures 3.1 and 3.2 the points A and B, for class Al and B2 maps, l i e to the left of the 6 = 45° line. This shows that DTMs obtained by digitizing 1:50,000 class Al or B2 maps are adequate for Ad < 1/2 pixel position accuracy for these two types of imagery. are shown in 45°. With a be relaxed as For SPOT side looking (26° resolution DTMs made by devices off nadir) and Seasat SAR imagery, such as the Gestalt photomapper have high to be TABLE 3.1 RESOLUTION AND WORST RELIEF DISPLACEMENT FACTOR Sen s o r D i s p l a c e m e n t d = f h R e s o l u t i o n r L a n d s a t - 4 TM d = 0.14 h 30 m SPOT MLA ( n a d i r l o o k i n g ) d = 0.042 h 20 m SPOT MLA ( s i d e x l o o k i n g ) d = 0.60 h 26 m SPOT PLA ( n a d i r l o o k i n g ) d = 0.042 h 10 m SPOT PLA ( s i d e l o o k i n g x) d = 0.60 h 13 m Sea s a t SAR d = 2.3 h 25 m 1 In the SPOT side looking mode, the off nadir angle is 26°. 76 Ap METRES FIGURE 3.1 DTM ACCURACY REQUIREMENT FOR LANDSAT-4 TM SENSOR (A AND B ARE TYPICAL Ax AND Ap FOR A ACCURACY AND B ACCURACY 1:50,000 MAPS RESPECTIVELY) FIGURE 3.2 DTM ACCURACY REQUIREMENT FOR SPOT SENSORS, NADIR LOOKING (A AND B ARE TYPICAL Ax AND Ap FOR A AND B ACCURACY MAPS RESPECTIVELY) 77 FIGURE 3.3 DTM ACCURACY REQUIREMENT FOR SPOT SENSORS, OFF-NADIR LOOKING ANGLE = 26° Ap METRES Ax METRES FIGURE 3.4 DTM ACCURACY REQUIREMENT FOR SEASAT SAR SENSOR 78 used i f half a p i x e l accuracy i s desired. This i s shown i n Figures 3.3 and 3.4 where the corresponding points A and B shown i n the other two f i g u r e s are now to the right of the 9 = 20° l i n e . Using a good approximation (say a 16x16 point kernel) of the sine function to sample the raw DTM to the output coordinates, the inaccuracies i n Ax, Ay, and Az should be about the same i n both coordinate systems. Hence, the accuracy requirements just discussed should also apply to the raw DTM. 3.4.2 DTM Interpolation i n S c r o l l i n g Buffer The DTM s c r o l l i n g buffer shown i n Figure 2.15 contains enough lin e s to ensure that a l l the heights of an output l i n e can be found. However, the DTM pointer f o r each output point i s usually not an integer and therefore, DTM i n t e r p o l a t i o n i s required to compute the heights. If nearest neighbour i n t e r p o l a t i o n i s used, the height extraction process w i l l be very f a s t . A method i s presented here, with examples, to determine when nearest neighbour i n t e r p o l a t i o n can be used without s i g n i f i c a n t l y a f f e c t i n g the geometric mapping accuracy. (Caution: If the DTM i s also used for radiometric cor r e c t i o n , nearest neighbour i n t e r p o l a t i o n may be inadequate.) The method consists of two steps: (a) Determine the error due to nearest neighbour i n t e r p o l a t i o n . If t h i s error i s less than a p i x e l , then nearest neighbour 79 interpolation can be used. Let this error be e^ . (b) If error e^ is greater than a pixel, compare i t with map induced error &^ • If masks out e^ , then use nearest neighbour interpolation, otherwise use bilinear interpolation. The error e^ , can be computed as follows. The maximum error made by using nearest neighbour interpolation in x and y i s : Ax = Ay = r'/2 (3.7) where r' now is the resolution of the interpolated DTM. Therefore maximum e^ i s , from Equation 3.2, given by: e1 = f / ( A x ) 2 + ( A y ) 2 tan 6 In terms of pixels: /Z r' ej = f V | | tan 6 (3.8) In the f i r s t example, assume that: (a) r' = r , i.e., the interpolated DTM has the same pixel spacing as the raw image. (b) Ax and Ay have uniform distributions. (c) 6 has a maximum value of 45° and tan 9 has uniform distribution. Then: f2 r' tan 45' 2 r/3 /5 = f 12 I 6 pixels RMS 80 The computation of i s given by Equations 3.3 and 3.5 in which Ap = J 2 2 V A X + Ay , and Az as a function of map class is given in Section 3.1. Two DTM sources are used, namely, photogrammetric stereo and 1:50000 Class Al maps. The results of e^ and e.^ a r e shown in Table 3.2. Here nearest neighbour interpolation can be used for a l l the imagery mentioned since e^ is less than 1 pixel RMS. Note that in the last row, subpixel rectification accuracy can no longer by achieved for Seasat imagery due to the magnitude of e^-A second example is presented here with the following changes: (a) r' = 2r (b) 1:50,000 Class B2 maps are used instead of Class Al maps. The results are presented i n Table 3.3. Of particular interest are the Seasat SAR cases. If the DTM source i s photogrammetric stereo, nearest neighbour interpolation should not be used since the error e^ is now more than 1 pixel RMS. However i f the DTM source is from the Class B2 map, then nearest neighbour can be used since error masks out e^. Similar analysis can be performed for various slopes and DTM resolutions. TABLE 3.2 DTM RESAMPLING IN SCROLLING BUFFER ( I ) r' = r DTM Source S a t e l l i t e Sensor e. P i x e l s RMS e 2 P i x e l s RMS To t a l E r r o r P i x e l s RMS Recommended DTM I n t e r p o l a t i o n Method Landsat-4 TM 0.03 0.03 Nearest neighbour Photo-gramme-t r i c s t e r e o SPOT ( n a d i r l o o k i n g ) SPOT ( s i d e l o o k i n g ) 0.01 0.14 n e g l i g i b l e 0.01 0.14 Nearest neighbour Nearest neighbour Seasat SAR 0.54 0.54 Nearest neighbour 1:50,000 Cla s s A l maps Landsat-4 TM SPOT MLA ( n a d i r l o o k i n g ) 0.03 0.01 0.08 0.036 0.09 0.037 Nearest neighbour Nearest neighbour SPOT PLA ( n a d i r l o o k i n g ) 0.01 0.072 0.073 Nearest neighbour SPOT MLA ( s i d e l o o k i n g ) 0.14 0.39 0.41 Nearest neighbour SPOT PLA ( s i d e l o o k i n g ) 0.14 0.78 0.79 Nearest neighbour Seasat SAR 0.54 1.56 1.65 Nearest neighbour TABLE 3.3 DTM RESAMPLING IN SCROLLING BUFFER ( I I ) r ' =2r DTM Source S a t e l 1 i t e Sensor e a P i x e l s RMS ea P i x e l s RMS T o t a l E r r o r P i x e l s RMS Recommended DTM I n t e r p o l a t i o n Method Photo-gramme-t r i c s t e r e o Landsat-4 TM SPOT ( n a d i r l o o k i n g 0.06 0.02 N e g l i g i b l e 0.06 0.02 N e a r e s t n e i g h b o u r N e a r e s t n e i g h b o u r SPOT ( s i d e l o o k i n g ) 0.28 0.28 N e a r e s t neighbour Seasat SAR 1.08 1.08 B i l i n e a r 1:50,000 C l a s s B2 Maps Landsat-4 TM SPOT M L A ( n a d i r l o o k i n g ) 0.06 0.02 0.16 0.072 0.10 0.04 N e a r e s t neighbour N e a r e s t n e i g h b o u r SPOT PLA ( n a d i r l o o k i n g ) 0.02 0.14 0.14 N e a r e s t n e i g h b o u r SPOT MLA ( s i d e l o o k i n g ) 0.28 0.78 0.83 N e a r e s t neighbour SPOT PLA ( s i d e l o o k i n g ) 0.28 1.56 1.58 N e a r e s t n e i g h b o u r Seasat"SAR 1.08 3.10 3.28 Ne a r e s t n e i g h b o u r 83 CHAPTER 4 SCANNER IMAGERY 4.1 Introduction Scanner imagery i s of two types: mechanical and e l e c t r o n i c . Examples of the former are Landsat multispectral scanner (MSS) and TIROS-N imagery while an example of the l a t t e r i s SPOT imagery. A multispectral scanner (MSS) system consists e s s e n t i a l l y of three parts: a mechanical, an o p t i c a l and an e l e c t r o n i c system as i l l u s t r a t e d i n Figure 4.1. The l i n e scan i s performed i n the mechanical system by a rot a t i n g or o s c i l l a t i n g mirror. The o p t i c a l system consists of a beam s p l i t t i n g device which divides the incident l i g h t i n t o several s p e c t r a l bands. The e l e c t r o n i c system i s an array of detectors. Each s p e c t r a l band i s focused onto a set of detectors. The s a t e l l i t e moves forward during the scanning process, hence producing a two-dimensional image. Figure 4.2 shows a t y p i c a l way i n which imagery i s obtained. An e l e c t r o n i c scanning sensor employs a large number of detectors ( i n the thousands), as shown i n Figure 4.3,. i n the f o c a l plane so that no o s c i l l a t i n g mirror i s required. The detectors are sampled e l e c t r o n i c a l l y i n sequence. For t h i s reason, the e l e c t r o n i c scanner i s also c a l l e d a "push-broom" scanner. The main difference between the e l e c t r o n i c scanner and the mechanical scanner l i e s i n the physical arrangement of the detectors. In the former, each detector "sees" only one Mechanical O p t i c a l «. E l e c t r o n i c R o t a t i n g or O s c i l l a t i n g M i r r o r Beam S p l i t t e r D e t e c t o r s * /7\ Recorder J i e l d of View Ground FIGURE 4.1 SCANNER SENSOR SYSTEM 85 FIGURE 4.2 TYPICAL MECHANICAL SCANNING SENSOR FOCAL POINT SAMPLING SEQUENCE EARTH SURFACE FIGURE 4.3 TYPICAL ELECTRONIC SCANNING SENSOR 87 spot on the ground along a line of the image, while in the latter, a detector "sees" spots on the ground along a column of the image. Raw scanner data do not need to be processed to be rendered as an image. However without processing, the raw imagery w i l l be geometrically distorted due to earth rotation, nonlinear sensor sweep, panoramic distortion, variation in sa t e l l i t e attitude angles (pitch, r o l l and yaw), scan skew and the lik e . Furthermore, relief displacement i s significant in Landsat-4 TM imagery and w i l l be significant in SPOT imagery. Digital rectification corrects these distortions using information determined from a priori data and measurement data such as GCPs. Panoramic distortion and distortions due to nonlinear scanner sweep and scan skew are systematic; i.e., their effects can be predicted in advance. Distortions due to earth rotation are a function of s a t e l l i t e latitude and orbit and thus can be predicted from tracking data. While distortions due to attitude and altitude changes.are also systematic, the parameters required to account for them are not known. If these distortions are to be corrected, their effects must be measured for each image. The measurement technique involves identifying in the raw image recognizable geographic features (GCPs) whose map coordinates are known. Using the image locations of these points, the image distortion caused by the sensor platform can be determined. The purposes of this chapter are as follows: (a) Review and compare various rectification algorithms pertinent to 88 MSS imagery. The study shows that a t t i t u d e modelling i s superior to other methods i n terms of accuracy with a reasonable number of GCPs. (b) E s t a b l i s h how GCPs can be used to r e f i n e the at t i t u d e angles and or b i t parameters. (c) Identify low and high frequency d i s t o r t i o n s common to a l l scanner imagery. (d) Id e n t i f y geometric d i s t o r t i o n s s p e c i f i c to Landsat-4 TM and SPOT imagery. (e) Identify a t t i t u d e and o r b i t parameters to be used i n the recursive estimator i n the navigation phase. The recursive estimator i s tested with Landsat and TIROS-N imagery. These two types of imagery have d i f f e r e n t resolutions and panoramic d i s t o r t i o n s . ( f ) Investigate the problem of r e l i e f displacement, e s p e c i a l l y i n Landsat-4 TM and SPOT imagery. Synthesize r e l i e f d i s t o r t e d imagery and r e c t i f y the r e s u l t i n g d i s t o r t e d imagery. (g) Develop an i n t e r p o l a t i o n algorithm to handle imagery f o r which p i x e l spacings are not uniform. This i s applicable to Landsat-4 TM imagery where scan gaps are present. 4.2 Overview of Geometric E f f e c t s 4.2.1 Common Di s t o r t i o n s Raw scanner imagery i s geometrically d i s t o r t e d due to earth shape, 89 scanner geometry and s a t e l l i t e motion. These c h a r a c t e r i s t i c d i s t o r t i o n s can be c l a s s i f i e d i n t o two types: low frequency and high frequency. The low frequency d i s t o r t i o n s are: (a) Panoramic. The p i x e l s i z e increases with the scan angle which the imaging sensor makes with the v e r t i c a l . (b) A t t i t u d e . P i t c h causes a l i n e s h i f t , r o l l causes a p i x e l s h i f t while yaw causes an image r o t a t i o n . (c) Earth r o t a t i o n . Earth r o t a t i o n causes an image skew. (d) Mirror scan. The scanning rate of a mechanical scanning sensor i s nonuniform, hence the sensor pointing angle cannot be determined as a simple l i n e a r function of the p i x e l number. For e l e c t r o n i c scanners, a s i m i l a r problem may occur because scanning i s not done at a rate corresponding to equal increments i n scan angle. (e) Band o f f s e t . There exists a time delay i n sampling between bands and the sensors of d i f f e r e n t bands are positioned d i f f e r e n t l y on the f o c a l plane. These factors contribute to an o f f s e t between bands. The high frequency components are: (a) Sensor o f f s e t . In MSS and TM the detectors are deployed i n a matrix pattern so that one swath contains more than one l i n e of imagery data. There e x i s t s a time delay i n sampling the corresponding p i x e l between the detectors. In addition, the 90 corresponding detectors are p o s i t i o n e d d i f f e r e n t l y on the f o c a l plane. These f a c t o r s c o n t r i b u t e to sensor o f f s e t which i s d i f f e r e n t from l i n e to l i n e w i t h i n a swath. (b) L i n e length v a r i a t i o n . The l i n e length i n a mechanical scanner v a r i e s due t o v a r i a t i o n s i n the scanning mirror speed. T h i s problem does not e x i s t i n l i n e a r array imagery which i s obtained by e l e c t r o n i c scanners. But, i n t h i s case, a s i m i l a r problem a r i s e s from the f a c t that gaps occur between groups of d e t e c t o r s . Thus the output p i x e l spacing i s regu l a r only w i t h i n each group of d e t e c t o r s . (c) E a r t h r o t a t i o n . E a r t h r o t a t i o n between de t e c t o r s w i t h i n a swath has to be taken i n t o c o n s i d e r a t i o n . E a r t h r o t a t i o n i s not l i n e a r from l i n e to l i n e as i l l u s t r a t e d i n Figure 4.4. This problem only occurs i n imagery with more than one l i n e per swath. (d) S a t e l l i t e a l t i t u d e and v e l o c i t y . S a t e l l i t e a l t i t u d e v a r i a t i o n produces a change i n the scan l i n e spacing w i t h i n a swath, w h i l e s a t e l l i t e v e l o c i t y v a r i a t i o n produces a change i n swath centre spacing. The two f a c t o r s when coupled together cause i r r e g u l a r l i n e spacings i n the along track d i r e c t i o n . (e) R e l i e f displacement. T h i s causes a p i x e l s h i f t i n the along scan d i r e c t i o n and i s discussed i n more d e t a i l i n S e c t i o n 4.4.4. The b i l i n e a r c o e f f i c i e n t s used i n remapping can only account f o r low frequency d i s t o r t i o n s . High frequency components are added to the r e s u l t of the b i l i n e a r transformation i n the f i r s t and second image Ground c o v e r a g e due t o b i l i n e a r c o e f f i c i e n t s assumption' A c t u a l swath c o v e r a g e on ground FIGURE 4.4 EARTH ROTATION EFFECT 92 processing passes. 4.2.2 Landsat-4 TM Imagery S p e c i f i c D i s t o r t i o n s The major TM imagery d i s t o r t i o n , i n addition to the ones mentioned i n Section 4.2.1, i s the underlap and overlap of swath data. This i s caused by v a r i a t i o n s i n s a t e l l i t e a l t i t u d e (from 685 to 740 km) and v e l o c i t y , j i t t e r , bowtie e f f e c t and by earth r o t a t i o n as described i n Appendix A. The overlap and underlap i s known as scan gap and i s also described and analyzed i n d e t a i l i n Appendix A. The scan gap size varies between ±2 data points and t h i s has a s i g n i f i c a n t impact on i n t e r p o l a t i o n ( i n the second pass) since the data samples are no longer uniformly spaced. A s p e c i a l i n t e r p o l a t i o n kernel to minimize the radiometric error i n i n t e r p o l a t i o n w i l l be discussed i n Section 4.5. 4.2.3 SPOT Imagery S p e c i f i c D i s t o r t i o n s The SPOT sensor i s described i n Appendix B. The geometric d i s t o r t i o n s associated with SPOT are as mentioned i n Section 4.2.1, but the following points should be noted: (a) Since each swath contains only one l i n e of imagery data, the , high frequency earth r o t a t i o n effect i n the along scan d i r e c t i o n does not e x i s t . Also a l t i t u d e and v e l o c i t y v a r i a t i o n s no longer produce nonuniform l i n e spacings i n the across scan d i r e c t i o n . (b) P i x e l r e s o l u t i o n i n the along scan d i r e c t i o n i s a function of j o f f nadir angle as demonstrated i n Appendix B. 93 SPOT imagery suf f e r s no high frequency d i s t o r t i o n s other than r e l i e f displacement. 4.2.A R e l i e f Displacement Any imaging system introduces r e l i e f displacement, the magnitude of which depends on the imaging parameters. The exact formulation f o r r e l i e f displacement i n units of raw input p i x e l s can be derived with the aid of Figure A.5. The displacement i n scan angle 8 can be expressed e x p l i c i t l y by: Hz s i n B tan A3 = — • ,y (A.l) R H cos 6 - (R + z ) \ 1 - (H s i n 3 / R) where z i s the t e r r a i n height, R i s the l o c a l earth radius and H i s the l o c a l o r b i t radius. Note that R >> z and Ag i s small enough so that tan A3 - A3. Then: A3 H s i n 8 (A.2) R H cos 3 - R Jl - (H s i n 3 / R)' A3 where — i s i n units of angular measure per unit height. This can e a s i l y be converted i n t o p i x e l s per unit height by: Ap AP A3 z A3 z where Ap i s r e s o l u t i o n . the displacement i n pix e l s and Ap/A3 i s the angular S a t e l l i t e E a r t h s u r f a c e FIGURE 4.5 RELIEF DISPLACEMENT IN SCANNER IMAGERY 95 The angle 3 i s a function of p i x e l number, mirror scan nonlinearity and r o l l angle. The r o l l angle cannot be ignored e s p e c i a l l y i n the case of SPOT PLA imagery when the off nadir angle i s at i t s maximum, i . e . 26°. In t h i s case a r o l l angle of 0.5° would introduce a 10 metre s h i f t per km height. 4.3 R e c t i f i c a t i o n Transformation 4.3.1 Closed Form So l u t i o n T h i s section develops a closed form transformation between the instantaneous sensor viewing angle and the geographic coordinates of any point on the earth's surface. This i s referred to as the forward transformation. Derivation of the transformation, det a i l e d i n Appendix E, s t a r t s by postulating a unit scan vector viewing a GCP i n the sensor platform coordinate system. This vector, shown i n Figure 4.6, i s then transformed into the s a t e l l i t e coordinate system and f i n a l l y to the ECR system. Let u be the unit vector pointing from the s a t e l l i t e to a point on the earth, R the ground point vector and R the s a t e l l i t e p o s i t i o n 8 -s vector as shown i n Figure 4.6. The transformation from the scan vector Du to the ground point at R i s given by: where D i s the distance between the s a t e l l i t e and the image point. Appendix E shows that D can be e x p l i c i t l y expressed as: R = R + Du -g -s (4.3) (A.4) where 2 2 e' = (a - b ) / b 2 EARTH CENTRE FIGURE 4.6 SCAN VECTOR GEOMETRY 97 a = semi-major axis of the earth b = semi-minor axis of the earth = z-component of Rg u = z-component of u z -and T = (uT.R + e'R u ) 2 - (1 + e' u 2)(|R | 2 + e'R 2 - a 2) (4.5) -s z z z -s z The functional form of u is developed in Appendix E where i t is expressed in terms of sa t e l l i t e orbit position, s a t e l l i t e attitude angles, sensor deployment angle and scan angle. 4.3.2 Block Size Once the attitude and orbit parameters are known i t is quite straight-forward to transform from line, pixel coordinates to ECR coordinates. But the inverse transformation is required to perform the remapping which is performed into three passes. The transformation in each pass is performed by a pair of piecewise bilinear polynomials. The f i r s t and second intermediate images are segmented into blocks and a set of bilinear remapping coefficients is associated with each block. The remapping accuracy decreases with increased block size. A trade<-off study has been performed between remapping accuracy and block size and this section gives the results of the study. In the investigation, nominal sensor and orbit parameters were used. The block size and RMS error quoted in the figures to be presented are in ground distance in metres, which when combined with the 98 required output pixel size would give terms of output pixel and output determined experimentally. Figure 4.7 shows the results for Landsat-1, -2 and -3 MSS imagery with resolution of about 79 metres X 79 metres. Hence a reasonable output pixel size would be 50 metres X 50 metres. From the figure, i f a 0.1 pixel (5 metres) RMS error is required then a block size of 160 X 160 (8 km X 8 km) has to be used. Figure 4.7 also shows the corresponding results for Landsat-4 MSS and TM imagery. The two graphs ill u s t r a t e that the block size required for Landsat-4 MSS can be larger than that for Landsat-1, -2 and -3 for the same RMS error. This i s because in -2 Landsat-1, -2 and -3 the d r i f t in each attitude angle is 10 degree per second whereas in Landsat-4 i t is 10 ^ degree per second. The attitude angles are thus a significant factor in the block size. Figure 4.8 shows the results for SPOT PLA and MLA imagery (Section 4.2.3). In this figure the block size is also a function of the off nadir pointing angle. From this the block size for any off nadir pointing angle can be determined for a fixed RMS error. 4.4 Recursive Estimator in Navigation 4.4.1 State Vector Parameters Rectification accuracy of s a t e l l i t e borne line scan imagery is a function of s a t e l l i t e orbit and attitude accuracies. Satellite orbit the block size and RMS error in line spacings. The graphs were RMS ERROR IN METRES 12 + BLOCK SIZE IN METRES FIGURE 4.7 ACCURACY VERSUS BLOCK SIZE FOR LANDSAT SENSORS RMS ERROR IN METRES 25 2 0 -15.. 10--X. 0 = 0.45 radian a 0 = 0.34 radian o 9 = 0.23 radian A 9 = 0.11 radian 9 = 0.00 radian 4-I 1 2000 4000 6000 BLOCK SIZE IN METRES 8000 FIGURE 4.8 ACCURACY VERSUS BLOCK SIZE FOR SPOT SENSORS 101 information can be determined from ground based tracking measurements and a t t i t u d e data are derived from onboard measurements and are telemetered to the ground. When such a t t i t u d e measurements are made, they are not precise enough f o r a p r i o r i subpixel l e v e l r e c t i f i c a t i o n . The s t a t e vector J£ (see Section 2.2.5) then contains the time s e r i e s c o e f f i c i e n t s of the a t t i t u d e angles. A further v a r i a b l e i n the state vector i s the s a t e l l i t e a l t i t u d e above the earth's surface. Landsat-1 and -2 imagery has been r e c t i f i e d to subpixel accuracy, with a cubic polynomial i n each a t t i t u d e angle and a constant i n s a t e l l i t e a l t i t u d e i n the state vector. 4.4.2 Extension to TIROS-N Imagery 4.4.2.1 Objective Guided by Konecny's survey [KONE76] and the findings i n the r e c t i f i c a t i o n of Landsat-1 and -2 imagery, the recursive estimator technique has been extended i n t h i s thesis to the r e c t i f i c a t i o n of TIROS-N imagery. Such imagery su f f e r s from a much higher panoramic d i s t o r t i o n , covers a much larger area for a t y p i c a l scene and has a considerably lower r e s o l u t i o n than the Landsat imagery. The objective i s to v e r i f y t h i s method independent of r e s o l u t i o n and panoramic d i s t o r t i o n . As an example, a TIROS-N Advanced Very High Resolution Radiometric (AVHRR) image was chosen f o r study. Simulation r e s u l t s have shown that the technique could y i e l d subpixel r e c t i f i c a t i o n geometric accuracy 102 with about 10 to 20 GCPs under favourable conditions. A.A.2.2 TIROS-N C h a r a c t e r i s t i c s In order to analyze the transformation accuracy, i t i s important to understand the s a t e l l i t e o r b i t and i t s AVHRR sensor c h a r a c t e r i s t i c s . The nominal values of the or b i t and sensor parameters f o r the AVHRR are shown i n Table 4.1. Transformation accuracy with a time series model was shown by Caron and Simon [CAR075] to be within one p i x e l (79 metres) f o r Landsat MSS imagery. However, o r b i t and sensor s p e c i f i c a t i o n s f o r the TIROS-N AVHRR are somewhat d i f f e r e n t from those for Landsat MSS. The main di f f e r e n c e r e l a t i n g to transformation accuracy i s the viewing angle of 110.8" (compared to 11.5 i n Landsat MSS). Therefore, a scan l i n e i n a TIROS-N image covers about 3000 km compared to 185 km i n a Landsat MSS image. The wider area of coverage introduces more panoramic d i s t o r t i o n i n an AVHRR image than i n a Landsat MSS image. The ef f e c t of earth curvature i s also considerably larger than i n Landsat. T y p i c a l l y f o r AVHRR imagery a p i x e l at the scene centre covers 0.78 km compared with 4.5 km at the edge, giving a r a t i o of about 5.8 to 1; the corresponding r a t i o f o r Landsat imagery i s 1.01 to 1. A graph of resol u t i o n versus p i x e l coordinates i s shown i n Figure 4.9. The l i n e spacing i s about 1.1 km as compared to 79 metres i n Landsat MSS. 1 0 3 TABLE 4.1 TIROS-N CHARACTERISTICS Semi-major a x i s E c c e n t r i c i t y I n c l i n a t i o n Nodal p e r i o d F i e l d o f view Scan p e r i o d Number p i x e l s p e r scan P i x e l s a m p l i n g r a t e 7211 km 0.001 98.7° 101.58 minutes 110.8° 1/6 second 2048 39936 samples p e r second 104 250" 50"b" : 75*0 1000 DISTANCE FROM CENTRE IN PIXELS FIGURE 4.9 RESOLUTION VERSUS DISTANCE FROM CENTRE FOR TIROS-N IMAGERY (ASSUMING ZERO ATTITUDE ANGLES) 105 4.4.2.3 Orbit Accuracy Analysis Geometric transformation is a function of the following parameters: (a) S a t e l l i t e inclination — Inaccuracy w i l l give an image rotation. (b) Nodal period and altitude — These give a scaling change in the image. These two parameters are related by Kepler's laws of motion. (c) Eccentricity — — Since this value is about 0.001, the s a t e l l i t e altitude can vary by about 14 km, taking the earth centre as a focus of the orbit. (d) Equatorial crossing — — This is the longitude at which the s a t e l l i t e crosses the equator. Any inaccuracy w i l l give an along scan error. (e) Elapsed time between equatorial crossing and image acquisition •' Inaccuracy w i l l give an along track error. The equatorial crossing has to be known to within 0.005° resulting in no more than 0.5 pixel error in the along scan direction. Also, the elapsed time has to be known to within 0.08 seconds resulting in no more than 0.5 line error in the along track direction. The altitude is assumed to be the major contribution to geometric transformation error in the orbital parameters and requires an update in the recursive f i l t e r . 106 4.A.2.4 Attitude Accuracy Analysis The recursive f i l t e r i n g technique estimates the attitude angles only to a certain accuracy. It can be shown that in the absence of yaw ( K) and pitch ( <f>), r o l l (to) w i l l have to be determined to a precision of I Aio | < 0.48 mrad (4.6a) i n order for the along scan pixel coordinate to be accurate to +0.5 pixel. Similarly for pitch (<{>), |A<fi| < 0.66 mrad (4.6b) w i l l result in ±0.5 line accuracy i n the across scan direction. Finally for yaw (K) , |AKI < 0.45 mrad (4.6c) w i l l result in ±0.5 line accuracy. If these accuracies for the attitude angles are attainable, then the geometric accuracy as a result of the combined effect w i l l be 73 x 0.5 = 0.86 pixels, thus subpixel accuracy w i l l be achieved. It is possible to relate the RMS deviations Au, Ao> and AK to the accuracy of designating GCPs, and thus arrive at an estimated mapping accuracy and the number of GCPs required. Assume that any GCP can be located to within 1 pixel (i.e.,1 instantaneous f i e l d of view,or IFOV) for designation purposes. Then for equally probable r o l l , pitch and yaw, a GCP measurement results in AOJ ~ ±0.96 mrad, Aa> 2- ±1.32 mrad and AK 2> ±0.90 mrad. For n GCP measurements accurate to ±1 pixel, the individual error of each attitude angle must be combined by the Gaussian rule so 107 that: Aio = +0.96 / J~n mrad n A<b = +1.32 / J~h~ mrad and AK = +0.90 / /n mrad n v To a t t a i n the a t t i t u d e accuracies shown i n Equations 4.6a, b and c f o r subpixel geometric accuracy, i t follows that n > 4. That i s , i n t h i s example at least 4 GCPs are required to reduce the attitude angles to acceptable l e v e l s . In p r a c t i c e , GCP designation of 1 p i x e l accuracy may be d i f f i c u l t to achieve. If each GCP measurement of m IFOV i s attained, 2 4m GCPs w i l l be necessary f o r the acceptable a t t i t u d e angles. Thus, i f i n p r a c t i c e ±2 IFOVs p r e c i s i o n per GCP measurement i s made, a t o t a l of 16 GCPs w i l l be needed. 4.4.2.5 Experimental Results and Discussions Experiments using the recursive estimator technique were performed with a TIROS-N image covering the western part of the United States. In the experiments, error i s defined: n 2 2 Residual RMS error = / I (Ax. + Ay. ) / n (4.7) 1=1 1 1 where n i s the number of GCPs processed, Ax_^ and Ay., are the differences i n p i x e l and l i n e coordinates f o r a p a r t i c u l a r GCP i , between the actual and estimated (using the updated state vector a f t e r processing n GCPs) values. 108 Verification Test 1 In this experiment each attitude angle was allowed to vary between ±0.5°as a cubic function of time. The time duration was selected to be about 200 seconds which corresponds to 1,200 lines of imagery data and to about 1,300 km in the along track direction. Nineteen GCPs were simulated. The object of the experiment was to study the performance as a function of the number of GCPs. The selection of the i n i t i a l covariance matrix P(0) is based upon the following assumptions: (a) The components of V. are uncorrelated. (b) The uncertainty in the zeroth order component is the same for each attitude angle and is bounded by the attitude measurement uncertainty given in the TIROS-N specification. This is 0.5°. (c) The uncertainty in the j-th order component is the same for each attitude angle and i t s contribution at any time to the attitude angle does not exceed the contribution made by the ( j - l ) t h component. (d) The uncertainty in the semi-major axis follows from the TIROS-N specification. This is 5 km. In the experiment the confidence value of GCP measurement K is f i r s t set equal to 0.5 pixel. Figure 4.10 shows the result of the experiment. . The f i n a l residual error reduces to less than 0.02 pixel. The value of £ was then set equal to 4 pixels and the result is also shown in Figure 4.10. After processing the GCPs, the residual error reduces to about 0.7 pixel. This error reduces further by processing more 109 P i t ch Uncerta inty 10"2 rad RMS Ro l l Uncerta inty 10"2 rad RMS Yaw Uncertainty 10"2 rad RMS A l t i t u d e Uncerta inty km RMS Residual fc»«.0 C-0.5 e r r o r 1.0 0.050 1n pixels WS 0.5 0.025 4-i 0.0 0.0 S 10 15 Number of GCPs processed 20 40 60 FIGURE 4.10 UNCERTAINTIES AND RESIDUAL ERROR VERSUS NUMBER OF GCPs PROCESSED, VERIFICATION TEST 1 110 simulated GCPs as shown i n the f i g u r e . Also i n both cases, i t was found that the at t i t u d e time seri e s and the semi-major axis a f t e r processing these GCPs agreed with the simulated values. A further test was performed i n the following manner. The geographic coordinates of the GCPs were corrupted by random noise equivalent to 3 p i x e l s RMS. This simulates combined GCP designation error and measurement e r r o r . The r e s u l t s show that the at t i t u d e time series and semi-major axis obtained a f t e r processing a l l the GCPs agreed well with the simulated values and the re s i d u a l error was about 3 p i x e l s RMS. This further shows that although the uncertainty of the attitud e angles can be reduced to any desired l i m i t by processing a large enough number of GCPs (shown i n Equations 4.6a, b and c ) , the r e s i d u a l error as defined i n Equation 4.7 may not necessarily reduce to zero due to measurement error. V e r i f i c a t i o n Test 2 V e r i f i c a t i o n test 1 was repeated but a bias of 0.1° was introduced i n t o the equatorial crossing and a s h i f t of 2 seconds was introduced i n t o the elapsed time from equatorial crossing. These perturbations correspond to 10 p i x e l s and 12 l i n e s s h i f t , or to an addit i o n a l p i t c h of 0.015 rad and an a d d i t i o n a l r o l l of 0.0096 rad. Figure 4.11 shows the RMS r e s i d u a l error as a function of the number of GCPs processed. Af t e r processing the 19 GCPs the RMS r e s i d u a l error of about 1.6 pixels was obtained and the ad d i t i o n a l r o l l and p i t c h were c o r r e c t l y r e f l e c t e d i n the output attitu d e time s e r i e s . I l l RESIDUAL ERROR IN PIXELS RMS 5 10 15 NUMBER OF GCPs PROCESSED FIGURE 4.11 VERIFICATION TEST 2 RESIDUAL ERROR IN PIXELS RMS 0 5 10 15 20 NUMBER OF GCPs PROCESSED FIGURE 4.12 OPERATIONAL RUN 112 The equatorial crossing was then s h i f t e d by 0.2° and the elapsed time from equatorial crossing s h i f t e d by 4 seconds. The RMS r e s i d u a l error increased to 4 p i x e l s . These experiments show that i f the measurement geometry, i . e . , matrix H i n the recursive f i l t e r , i s i n error i t i s not possible to absorb a l l the errors i n t o the state vector elements. One so l u t i o n i s to place a l l the uncertain terms which w i l l a f f e c t the transformation accuracy i n t o the state vector. Experiments with Landsat MSS imagery show that such perturbations can be absorbed into the p i t c h , r o l l and yaw components with very l i t t l e e f f e c t on the transformation accuracy. This i s due to the fact that a Landsat image covers a much smaller area than a TIROS-N image. An experiment was performed with a Baltimore / Washington scene (I.D. 20598-14590) i n which the scene centre was shifted by 30 p i x e l s and 22 l i n e s and the s a t e l l i t e heading was assumed to be constant over the e n t i r e scene. The r e s i d u a l error thus obtained was less than 50 metres RMS with 15 GCPs. This i s quite acceptable i f the user i s not interested i n how the a t t i t u d e angles vary. However, the yaw angle v a r i a t i o n over the scene was found to be 0.25°. This value exceeds the s p e c i f i c a t i o n and i s due to the fact that the s a t e l l i t e heading was erroneously assumed to be a constant. The s a t e l l i t e heading was then allowed to vary as i t should, and as a result the yaw v a r i a t i o n dropped to an acceptable value of 0.017°. The r e s i d u a l error was also less than 50 metres RMS. This suggests that the attitu d e angles cited by Caron and Simon are the a t t i t u d e angles plus perturbation e f f e c t s due to other parameters, such as 113 s a t e l l i t e i n c l i n a t i o n , equatorial crossing, elapsed time and o r b i t semi-major a x i s . Operational Run In t h i s experiment 19 GCPs were marked and d i s t r i b u t e d randomly on the image from 38°N to 48°N i n l a t i t u d e . This area corresponds to about the same 1,200 l i n e s as i n the f i r s t two experiments. Mapsheets of scale 1:2,500,000 were used to mark the geographic coordinates of the GCPs (better r e s o l u t i o n maps were not a v a i l a b l e ) . The GCPs were marked at r i v e r confluences, points on lakes and r e s e r v o i r s , that i s a l l on water/land boundaries. It should be noted that water l e v e l s , e s p e c i a l l y i n the U.S. midwest, are susceptible to seasonal changes. More det e r m i n i s t i c and time invariant GCPs such as highway in t e r s e c t i o n s could not be i d e n t i f i e d due to the r e s o l u t i o n of TIROS-N AVHRR imagery. The value of P(0) was i n i t i a l i z e d i n the same way as i n V e r i f i c a t i o n Test 1. Results i n Figure 4.12 show that the residual error i s 2.8 a f t e r 19 GCPs. The value of K l i e s between 0.5 and 4 p i x e l s and each value depended on the confidence i n designating the GCP. Also i n designating each GCP the image was zoomed by cubic convolution by 7 times. Furthermore i n an a d d i t i o n a l test the image was extended to 2,000 lin e s and the RMS r e s i d u a l error remained at 2.8 p i x e l s . 114 Error Analysis The estimated error sources noted i n Table 4.2 can be combined to f u l l y account for the measurement error of 2.8 pixels RMS in the operational run. The table shows a combined error of 2.6 pixels RMS. In addition there could be two more sources of error. It is doubtful i f the equatorial crossing of the s a t e l l i t e and the elapsed time were accurate to the limit as discussed in Section 4.4.2.3. It has been shown that errors in these two parameters cannot be absorbed into the current state vector. Also each attitude angle may be of a higher degree than a cubic polynomial as a function of time 1. In future work, the state vector should contain more elements such as equatorial crossing, elapsed time and more coefficients for the attitude angles and altitude. Although the residual error in the operational test is about 2.8 pixels RMS, i t is fel t that the technique with high quality ground control points and maps has the potential for measurement errors of about 1 to 2 pixels RMS. 4.5 Experiments — — Interpolation i n the Presence of Scan Gap 4.5.1 Objective Appendix A shows that data samples for Landsat-4 TM imagery are not uniformly spaced in the across scan direction. If the sine """At the time of writing we are s t i l l waiting for the results of the attitude measurements on board the s a t e l l i t e . TABLE 4.2 ERROR ANALYSIS OF TIROS-N NAVIGATION E r r o r Source Map a c c u r a c y GCP mar k i n g on map s h e e t GCP marking on image A l t i t u d e e r r o r RMS E r r o r Comments 2.0 km 1:2,500,000 s c a l e 1.3 km 0.5 mm measurement e r r o r 0.5 km 0.5 p i x e l e r r o r i n image 1.0 km Due t o c o n s t a n t h e i g h t a s s u m p t i o n o v e r 1300 km 116 i n t e r p o l a t i o n kernel i s used f o r image i n t e r p o l a t i o n , then the corresponding low pass f i l t e r w i l l operate on a spectrum which does not accurately represent the image i r r a d i a n c e . Another i n t e r p o l a t i o n kernel w i l l have to be used. This kernel must be able to y i e l d a high throughput i n image i n t e r p o l a t i o n . One technique, based on cubic convolution and on creating uniformly spaced samples has been proposed by Prakash and Beyer [PRAK81] and simultaneously by Avery and Hsieh [AVER81]. In t h i s method, as in the case of uniformly spaced samples, cubic convolution i s used. Therefore the kernel covers four samples and i s a weighted and truncated sine f u n c t i o n . When these four p i x e l s span two swaths, at most two samples are not uniformly spaced. Two, at most, uniform samples are then created by using two p i x e l s on each swath. A 4-point cubic s p l i n e i s then used to performed the i n t e r p o l a t i o n to give the values of the created samples. Then cubic convolution i s applied to the uniform samples which include the created samples. Th e i r method suffers the following drawbacks: (a) A 4-point cubic convolution does not give an accurate ' i n t e r p o l a t i o n compared to a 16-point truncated sine function. . with Kaiser windowing (as demonstrated i n Section 2.5.1.6). (b) The 4-point cubic s p l i n e makes no attempt to preserve the s i g n a l spectrum. 117 This section presents a new method of i n t e r p o l a t i n g Landsat-4 TM imagery i n the presence of scan gaps. Experiments are performed on sine waves with various frequencies and scan gap s i z e s , on synthesized sine wave imagery and on synthesized TM imagery. Experimental r e s u l t s have v e r i f i e d that t h i s method y i e l d s an RMS i n t e r p o l a t i o n error of one l e v e l out of 256 l e v e l s i n 8-bit imagery data f o r the case of overlap. F a i l u r e to use a method such as the one presented here can lead to s i g n i f i c a n t radiometric e r r o r s , thus greatly reducing the u t i l i t y of geometrically corrected imagery data. 4.5.2 Theory of Proposed Method Let a band l i m i t e d s i g n a l f ( t ) be sampled at the Nyquist rate and divided into two uniformly sampled groups. S p e c i f i c a l l y , l e t the sample po s i t i o n s be denoted by: where T i s the sample spacing and p i s a non-negative integer. For convenience, separation of the two groups by an amount At i s assumed to occur at t = 0. Following Linfoot and Shepherd's work [LINF39] . Yen [YEN56] has shown that the s i g n a l f ( t ) can be accurately constructed by: t < 0 t > 0 (4.8) 00 f ( t ) = • Z f ( x ) * (t) (4.9) m=-°o where 118 (-1)"' r(At/T + p) P!r(t/T)rc(At-t)/Ti P + (At-t)/T P + t/T 1 1 » T T ID m 2 T (4.10) Therefore the kernel ^ (t) can be expressed as Gamma functions [K0RN61] of At, scan gap where At occurs, sample number m and sampling point. The Gamma function has the recursive property: For At = 0, i.e. , the sample points are uniformly spaced, degenerates to the sine function as expected. This method is applicable here for the following reason. Conforming to the size of the 16-point sine function in the case of uniformly spaced samples, the kernel size In this case also spans 16 points. Since each swath contains 16 samples, the kernel covers at most one scan gap. The sample further away from the interpolation position usually carries less weight. Also, the contribution by samples not under the kernel is insignificant. Therefore the sequence of sample points under the kernel can be viewed to consist of two groups separated by a gap. Each group can be assumed to extend i n f i n i t e l y since points after the 16th sample from the interpolation location are not under the kernel. The two in f i n i t e groups assumption then satisfies the condition cited in Yen's method. m T(n+1) = n T(n) (4.11) Figure 4.13 shows plots of the kernel for various values of sample (a) Gap size « -0.9 ( O Gap sl2e « -0.5 (e) 1 I 1 J Gap slze-^ «= 0.5 A « ' . 1 , 1 , i ' i ' i * (9) (1) T T 1 Gap slze^ .• 1.5 1 J 1 | 1 1 " Note: 1: * • resample position 2. T •1.0 (b) / Gap size • -0.9 (d) Gap size - -0.5 (f) 6ap size —v « 0.5 t » (h) 1 1 1 /"Gap size 1 ./:••.• r i , I • . . . (j) V Gap size » 0.9 .Gap size « 1.5 FIGURE 4.13 INTERPOLATION KERNEL OF YEN'S METHOD number and scan gap size. The following observations are made: (a) When the scan gap i s negative, the kernel weight decreases with distance of sample point away from the interpolation position except at the scan gap (e.g., Figure 4.13a). (b) When the scan gap is negative and interpolation is performed in the gap, the two nearest sample points weigh more and a l l the other sample points weigh less than the sine function kernel, as expected. (c) When the scan gap is larger than 1 sample and interpolation is performed in the gap, the kernel weights do not decrease as a function of sample distance within the 16 samples. The figure shows that two modifications have to be made: (a) When the gap is positive, extra samples have to be created so that the adjusted scan gap is reduced to a size when the kernel weights decrease appreciably as a function of sample distance. The adjusted scan gap is then always less than a predetermined size (see Section 4.5.4 for experimental results). (b) When the scan gap is negative and the gap is close to either end of the 16 sample points, the sine interpolation is used. This i s because , in Yen's method, the kernel weights are heavy at both ends of the gap and the kernel weights are s t i l l significant for a few more samples beyond the gap. Hence, i f the kernel is truncated to 16 points, significant kernel weights w i l l be eliminated. This has been proved in the experiment 121 (Section 4.5.4) to y i e l d a larger i n t e r p o l a t i o n error than the sine function i n t e r p o l a t i o n . Section 4.5.4 shows the experimental r e s u l t s which determine the scan gap p o s i t i o n where switching to the sine function i n t e r p o l a t i o n should occur. Extra samples can be created by a l o c a l cubic f i t using the four nearest points, i . e . , two points on each side of the scan gap. Let these data points be F^, F 2 , F^ and F^ at positions t j , t ^ and t ^ r e s p e c t i v e l y . Then i n t e r p o l a t i o n can be done by the divided differences method: 4 i-1 f ( t ) = E f [ t t ] n (t - t ) (4.12) i-1 j=l J where F[ ] i s the divided d i f f e r e n c e which can be defined r e c u r s i v e l y by: F[t.] = F. F [ t i + 1 , ... , t i + k ] - F [ t . , ... , t 1 + k _ j ] C i + k " h F [ t . , ... , t . + k ] = The cubic polynomial f i t can equivalently be represented by a Lagrange polynomial: 4 f ( t ) = E F * (t) (4.13) 1=1 4 t - t. where 4>.(t) = II _ 1 (4.14) 1 j=l,j=i i J 122 The additional interpolation error due to creating new samples is determined as follows. Let X '»i = 1 to n, be the values of the n new samples and X^ be the actual values of the function to be interpolated at the same locations. Also let , i=l to n, be the kernel weights at the n positions. Then additional interpolation error i s : " w (X - X ') 1=1 This error i s larger when interpolating in the scan gap than away from the scan gap because the kernel weights are heavier when the interpolation position i s closer to the gap. Experimental results (Section 4.5.4) also show that a Kaiser window is necessary for a better performance of the method. 4.5.3 Implementation Considerations The two groups of uniformly spaced samples are symmetrical about the midpoint of the scan. The symmetrical property reduces the computation of the kernel weights by a factor of two. It is not necessary to compute the kernel weight for every sample each time the interpolation is called. The kernel weights can be stored in a table the size of which is determined as follows. The kernel is a function of sample number, t and scan gap size. If: (a) sample spacing is divided into 32 subintervals and the kernel 123 size i s 16 samples, (b) a new sample is created when the gap size i s larger than 0.3 line (discussed in Section 4.5.4), (c) a sample is deleted when the gap size is less than or equal to -1 line, (d) uniform sample spacing is assumed when the gap size is between -0.1 and 0.1 line (discussed in Section 4.5.4), (e) no interpolation i s necessary i f the interpolation position i s within 1/64 of any sample point, and (f) scan gap size i s varied in steps of 0.05 line (see discussion below), then the number of coefficients, taking advantage of the symmetrical property, i s (1/0,5) x 10 x 31 x 16 x 8 = 79360 (I) (ID (III) (IV) where the term: (I) = number of quantized gap sizes, (II) = number of subpixel intervals minus 1, (III) = number of data points, and (IV) = number of scan gap size positions (taking advantage of the symmetrical property). The scan gap quantization step of 0.05 line i s obtained by the following experiment. Scan gaps of ±0.5 line were chosen for a sine wave with frequency of 6 pixels per cycle and interpolation .was done 124 with erroneous scan gap sizes. The results, presented in Table 4.3, show that the RMS error in interpolation due to a 0.05 line error in scan gap size i s s t i l l acceptable (about 1 level for underlap and 0.7 for overlap). Note that in Canada overlap occurs more often than underlap and the scan gap size is expected to be less than 0.3 line. The above calculation assumes that the new samples are actually created when required. By examining Equation 4.13, i t is obvious that the new samples need not be physically created, but the kernel weights in the two samples on each side of the scan gap need to be adjusted by <(> as shown in the equation. Therefore, the scan gap size for 12 of the 16 kernel weights varies from -0.9 to 0.3 line, and i t varies from -2 to 2 lines for the remaining four kernel weights. Therefore, the number of coefficients, i f the new samples are not physically created, is about: (1/0.5) x 31 x 8 x (10 x 16 + 40 x 4) - 158720 The required number of coefficients cited above w i l l be reduced by 20% i f switching to sine kernel is adopted under certain conditions as discussed i n Section 4.5.2. 4.5.4 Experiments with Sine Waves Sine waves are f i r s t tested for the following reasons: (a) The frequency and scan gap can be controlled easily and hence the interpolation can be evaluated as a function of frequency 125 TABLE 4.3 RMS ERROR VERSUS SCAN GAP SIZE QUANTIZATION ACTUAL SCAN ERRONEOUS SCAN GAP RMS ERROR GAP SIZE SCAN GAP SIZE SIZE ERROR -0.50 0.00 0.67 -0.51 0.01 0.67 -0.5 -0.52 0.02 0.68 ( o v e r l a p ) -0.53 0.03 0.69 -0.54 0.04 0.69 -0.55 0.05 0.70 -0.50 0.00 0.60 -0.51 0.01 0.63 0.5 -0.52 0.02 0.71 (underlap) -0.53 0.03 0.83 -0.54 0.04 0.99 -0.55 0.05 1.17 Note: Test was performed on s i n e wave w i t h a frequency of 6 p i x e l s per c y c l e . 126 and scan gap s i z e . (b) Using the sine waves, i t can be shown how other methods (e.g., i n t e r p o l a t i o n by sine function kernel 3") compare to the chosen method. (c) In t h i s part of the experiment, the following parameters i n the chosen method can be determined: the scan gap p o s i t i o n f o r switching to the sine kernel i n t e r p o l a t i o n , and the scan gap size above which new samples are created. Each sine wave i n t h i s part of the experiment consists of 32 data points. The f i r s t 16 and the l a s t 16 sample points are equally spaced (both groups have the same spacing A ^ ) . The spacing between the 16th and 17th data point ( A ^ ) r e f l e c t s the magnitude of the scan gap and scan gap size i s defined as A^ - A^ . In the experiment, the scan gap s i z e v a r i e s from -0.9 to 2.0 samples. Scan gap s i z e less than or equal to -1.0 sample represents redundant data and w i l l not a f f e c t the experimental r e s u l t . The TM modulation transfer function (MTF) has a s i g n i f i c a n t frequency response up to 0.16 cycle per p i x e l (the maximum i s 0.5 cycle •', per p i x e l ) . The experiment hence places a strong emphasis around t h i s p a r t i c u l a r frequency. Performance of the sine i n t e r p o l a t i o n and Yen's method are shown i n Figure 4.14. Frequency of the sine wave i s 0.16 cycle per p i x e l . The 1 Here the sine function i s referred to as the 16-point truncated sine function with Kaiser windowing. 127 (a) SCAN GAP SIZE * -0.9 (b) SCAN GAP SIZE = -0.6 10 8 6 4 2 RMS ERROR 9 10 11 12 13 14 15 16 (c) SCAN GAP SIZE = -0.3 i 0 BESAMPLE POSITION , .9 10 11 12 13 14 15 16 (d) SCAN GAP SIZE * -0.1 10 T R M S IERROR 8 . 6 •• 4 •• 2 .. V _ 7 RESAMPLE P O S I T I O N — i 1 1 10 8 .. 6 .. 4 •• 2 .. -0 RMS ERROR 9 10 11 12 13 14 15 16 (e) SCAN GAP SIZE = 0.1 BFsarPi F iPns i jr inni , h 9 10 11 12 13 14 15 16 ( f ) SCAN GAP SIZE = 0.4 10 4 RMS ERROR 8 e • 4 2 10 8 6 + RESAMPLE POSITION —3 10 l ' l li 13 l'4 IS* ll 0 RMS ERROR RESAMPLE POSITION —§ id ii '—fz—ft—fr-tt—rr YEN'S METHOD WITH YEN'S METHOD WITHOUT KAISER WINDOWING SINC INTERPOLATION KAISER WINDOWING FIGURE 4.14 RMS ERROR VERSUS INTERPOLATION POSITION 128 abscissa of each plot i s the i n t e r p o l a t i o n p o s i t i o n from 8 to 16 + half of scan gap. The kernel weights of other i n t e r p o l a t i o n positions are r e p l i c a s of the weights i n t h i s segment. Each plot i s the envelop of the RMS e r r o r . The following conclusions can be drawn from the f i g u r e : (a) Yen's method without Kaiser windowing gives a higher RMS error than Yen's method with Kaiser windowing. Therefore, the former method i s discarded. This i s an expected result since t h i s i s also true f o r the sine function i n t e r p o l a t i o n and the o v e r a l l shape of the Yen's method kernel i s s i m i l a r to that of the sine function. (b) The sine function i n t e r p o l a t i o n and Yen's method give approximately the same RMS error when the scan gap size i s between ±0.1 l i n e . Hence, f or t h i s range of gap s i z e , the sample spacings can be assumed to be uniform and the sine i n t e r p o l a t i o n i s used. (c) When the i n t e r p o l a t i o n p o s i t i o n i s such that the scan gap i s close to the end of the kernel and i t s s i z e i s negative, Yen's method gives a very large error because the kernel weights a f t e r the gap are s t i l l s i g n i f i c a n t . In t h i s case the sine i n t e r p o l a t i o n performs better. The point of switching i s when the scan gap i s at the following (marked with x) positions i n the kernel: x x x x —I 1 1 1— • • • —I 1 1 h 1 2 3 4 13 14 15 16 The RMS error i s zero when the i n t e r p o l a t i o n p o s i t i o n f a l l s on a data point. The envelop i s the locus following the maximum error between pairs of data points. 129 It remains to determine when new samples have to be created. Experimental results show that one sample should be created when the scan gap is between 0.3 and 1.3 samples, and two samples should be created i f the scan gap is larger than 1.3 samples. It is not necessary to create more than two samples because the scan gap is always less than 2.0 samples. Figure 4.15 gives comparisons between the sine interpolation, Yen's method and the proposed method (i.e., Yen + Kaiser windowing + creating of new samples) for four different frequencies. The figure shows that the chosen method is superior to the other two methods. Sharp turns exist at the gap size equal to 0.3 and 1.3 samples for the chosen method before new samples are added. Without adding the new samples the RMS error increases rapidly due to the fact that the kernel requires more than 16 points, i.e., the kernel weights have not subsided yet at the end of the 16-point kernel. 4.5.5 Scan Gap Simulation in Imagery The scan gaps were assumed to be due to, the following factors: s a t e l l i t e altitude variation, j i t t e r and scan skew. This section describes how these scan gap constituents were simulated. The simulated scan gap is used in the synthesis of distorted imagery (Section 4.5.6). 130 3.25 1.0 D.75 0.5 0.25 0.0 (a) FREQUENCY • 0.02 CYCLE PER PIXEL RMS ERROR (b) FREQUENCY • 0.07 CYCLE PER PIXEL 4.0 3.0 2.0 1.0 0.0 RMS ERROR / -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 SCAN GAF SIZE SCAN GAP SIZE 12.5 10.0 7.5 5.0 2.5 0.0 (c) FREQUENCY « 0.12 CYCLE PER PIXEL j 2 . 5 RMS ERROR 10.0 7.5 5.0 2.5 0.5 (d) FREQUENCY • 0.15 CYCLE PER P/IXEL RMS ERROR -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 " •-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 SCAN GAF SJZ£ SCAN GAP SIZE PROPOSED METHOD YEN'S METHOD SINC INTERPOLATION FIGURE 4.15 COMPARISON OF METHODS - RMS ERROR VERSUS SCAN GAP SIZE 131 1. A l t i t u d e V a r i a t i o n This causes a constant scan gap s i z e , i n the absence of j i t t e r and scan skew, given by: where denotes the i n t e r l i n e spacing within a swath and Aj denotes the spacing between the two swaths. The value of G i n the experiment varied between ±2.0 l i n e s . A slanted road due to a p o s i t i v e and negative value of G i s shown i n Figure 4.16. A negative G means the s a t e l l i t e i s at a higher orbi t than nominally, and hence the sensor images a wider portion of the road with i t s 16 detectors. As a r e s u l t , the road i s segmented as shown i n Figure 4.16a. This corresponds to overlap. On the contrary, a p o s i t i v e value of G shows underlap and the slant road becomes segmented i n the way shown i n Figure 4.16b. 2. J i t t e r J i t t e r has frequencies i n the range from 2 to 125 Hertz, i t can be modelled as a ser i e s of sine waves at the output of a t y p i c a l band pass f i l t e r centred at 60 Hertz. The 3 db points at 30 and 90 Hertz are shown i n Figure 4.17. Hence, the scan gap i s : N 2 d(t) = c z {exp[(-60 + 12001/N) /1300]} sin[-2.ir( 120it/N) ] i=l where t = time, N = t o t a l number of sine waves, and c = m u l t i p l i c a t i v e constant to adjust the gap s i z e . ORIGINAL ROAD 4.16 EFFECT OF OVERLAP AND UNDERLAP ON A ROADWAY 1 Hertz 0 36 60 90 120 FIGURE 4.17 JITTER MODEL 134 In the experiment, the value of c was chosen such that the j i t t e r contribution to the scan gap varied from -0.5 to 0.5 line. 3. Scan Skew The scan skew was chosen so that i t s contribution to the scan gap was between +0.1 line. 4.5.6 Experiments with Sine Wave Imagery Sine wave patterns were used i n the simulation of imagery because they can synthesize linear features such as roadways and hence the effect of scan gap on these features can be visualized. Figure 4.18 depicts linear features on the ground running at an angle of 45° to the ve r t i c a l . The sine wave has a frequency of 6 pixels per cycle. This frequency was chosen because the MTF of the TM sensor dictates that i t takes about 3 pixels to respond from one intensity extreme to the other. The simulated TM image is shown in Figure 4.19, where the scan gap ranges from -1 to 1 line. Hence, this corresponds to overlap and the sine patterns are segmented in a manner similar to Figure 4.16a. Two methods of reconstruction were used - — sine interpolation and the proposed method. The results are shown in Figures 4.20 and 4.21. The difference images, each multiplied by a factor of 4, are shown in Figures 4.22 and 4.23. FIGURE 4.18 SINE WAVE IMAGE FIGURE 4.19 DISTORTED SINE WAVE IMAGE WITH OVERLAP FIGURE 4.20 RECTIFIED SINE WAVE IMAGE (OVERLAP) WITH SINC INTERPOLATION FIGURE 4.21 RECTIFIED SINE WAVE IMAGE (OVERLAP) WITH PROPOSED METHOD 139 FIGURE 4.22 RADIOMETRIC ERROR IN RECTIFIED SINE WAVE IMAGE (OVERLAP) WITH SINC INTERPOLATION FIGURE 4.23 RADIOMETRIC ERROR IN RECTIFIED SINE WAVE IMAGE (OVERLAP) WITH PROPOSED METHOD 141 The same image, but with underlap, was also simulated with scan gap s i z e s from 0.3 to 1.4 l i n e s . The result i s shown i n Figure 4.24 i n which the sine patterns are segmented i n a manner s i m i l a r to Figure 4.25 to 4.28. Figure 4.25 (sine i n t e r p o l a t i o n ) shows more obvious seams than Figure 4.26 (proposed method). The d i f f e r e n c e images c l e a r l y show that the proposed method i s better then the sine i n t e r p o l a t i o n f o r both the overlap and underlap cases. The quantitative r e s u l t s are given by the RMS value i n the d i f f e r e n c e images: Overlap Underlap Sine i n t e r p o l a t i o n 2.8 6.9 Proposed method 1.1 2.6 Thus, the proposed method i s about 2.5 times better than the 16-point sine i n t e r p o l a t i o n i n terms of RMS e r r o r . A plot of absolute error f o r a sequence of data points i n the d i f f e r e n c e images i s given i n Figure 4.29. T h i s figure c l e a r l y demonstrates the s u p e r i o r i t y of the proposed method. The method performs about two to three times better than the sine i n t e r p o l a t i o n . However, by just considering the neighbourhood of the scan gaps, the proposed method i s even better. The sine i n t e r p o l a t i o n w i l l introduce considerable a r t i f a c t s even i n the case of overlap. FIGURE 4.24 DISTORTED SINE WAVE IMAGE WITH UNDERLAP FIGURE 4.25 RECTIFIED SINE WAVE IMAGE (UNDERLAP) WITH SINC INTERPOLATION FIGURE 4.26 RECTIFIED SINE WAVE IMAGE (UNDERLAP) WITH PROPOSED METHOD FIGURE 4.27 RADIOMETRIC ERROR IN RECTIFIED SINE WAVE IMAGE (UNDERLAP) WITH SINC INTERPOLATION FIGURE 4.28 RADIOMETRIC ERROR IN RECTIFIED SINE WAVE IMAGE (UNDERLAP) WITH PROPOSED METHOD 10 1., ERROR MAGNITUDE G * -0 . 4 G = -0.7 G = -0.3 G = -0.6 • ^ v ^ w * w £ y . P i X E i I POSITION G » -0 . 4 (a) OVERLAP ERROR '; MAGNITUDE , ' l (b) UNDERLAP — Si n e I n t e r p o l a t i o n — Proposed Method FIGURE 4.29 ERROR MAGNITUDE 148 4.5.7 Synthesized TM Imagery TM Imagery with scan gaps can be synthesized using existing Landsat MSS imagery. A 512 x 512 Landsat-2 MSS image over Vancouver Island, B.C., Canada was selected for experiment. Figure 4.30 shows the scene which i s imaged over rugged terrain. A scan gap of 0.0 to -1.0 line corresponding to overlap was used to simulate the jittered image. In the simulation, a 16-point sine interpolation was used. Reconstruction was done by both the sine interpolation method and the proposed method. Figures 4.31 and 4.32 show the error images, each multiplied by a factor of 4. The error images show that the proposed method yields less interpolation error than the sine interpolation. The RMS error for the proposed method i s 0.8 level and that for the sine interpolation is 1.7 levels. RMS error also includes the error introduced i n simulating the jittered image by a 16-point sine interpolation. The error in the simulation is about 0.5 level RMS. When this simulation error is subtracted (in the Gaussian sense), the proposed method performs two to three times better than the sine interpolation and this is consistent with the sine wave experimental results. For the sine interpolation, the worst errors occur in the scan gap neighbourhoods; while in the proposed method (overlap case), the error is independent of positions relative to scan gaps. This again agrees with the sine . wave experimental results. The larger errors in the scan gap neighbours are due to the following: (a) Kernel weights in larger scan gaps do not taper off as fast as those in smaller scan gaps. (b) Error is attributed to creating new samples. When interpolating FIGURE 4.30 MSS IMAGE USED TO SYNTHESIZE TM IMAGE FIGURE 4.31 RADIOMETRIC ERROR IN RECTIFIED SYNTHESIZED TM IMAGE WITH SINC INTERPOLATION FIGURE 4.32 RADIOMETRIC ERROR IN RECTIFIED SYNTHESIZED TM IMAGE WITH PROPOSED METHOD 152 i n a gap neighbourhood, the kernel weights on these new samples are heavy. D i f f e r e n t scan gap sizes are tested with the scene and the RMS errors are shown i n Figure 4.33. The following should be noted: (a) The d i s t o r t e d images with scan gaps were simulated by using 16-point sine i n t e r p o l a t i o n which would introduce an RMS error of about 0.5 l e v e l . (b) In the overlap area, the proposed method y i e l d s an RMS error two to three times lower than that of the sine i n t e r p o l a t i o n method ( a f t e r subtracting the error i n the RMS sense introduced i n simulating the d i s t o r t e d image). This agrees with the sine wave experimental r e s u l t s . (c) In the underlap area, the RMS errors f o r both methods are about equal. This i s due to errors introduced i n the simulation of the d i s t o r t e d imagery. In addition to the 0.5 l e v e l RMS error, another error source could be that the data spacings i n the o r i g i n a l MSS image are not uniform. A scan gap occurs once every 6 l i n e s of data and can be as big as 0.1 l i n e . In the i n t e r p o l a t i o n process, creating new samples i s another error source. Nevertheless, the imagery experiments i n d i c a t e that the proposed method would y i e l d an i n t e r p o l a t i o n error of less than 1 l e v e l RMS i n the case of overlap. FIGURE 4.33 RMS ERROR VERSUS SCAN GAP SIZE IN SIMULATED TM IMAGE 154 4.6 Experiments — — Relief Displacement 4.6.1 Objective Rectification experiments of the f i r s t pass have been performed with synthesized imagery. The objective of this section is to present the experimental results which show that subpixel geometric accuracy can be achieved. Terrain models consisted of a simulated Gaussian terrain and a real DTM. Both Landsat-4 and SPOT sensor parameters are considered. In synthesizing the imagery, a Lambertian surface reflectance model [WOOD80 and H0RN81b] is assumed where the intensity of any particular point (x,y) is given by: , P P + q q + 1 I (x,y) = S , (4.15) 1 + p + q 1 + p s + q g where p = 9f(x,y) / 8x, q = 9f(x,y) / 9y, f(x,y) = altitude of point (x,y), and (-pg, -q ,1) is a vector pointing in the direction of the sun. Other surface reflectance models would yield different radiometric values for the synthesized imagery but would not significantly affect the study of geometric rectification accuracy. The results obtained with simulated Gaussian terrain are presented f i r s t . The advantage of using this terrain model is that slopes p, q required by the Lambertian model, and the intersection of the sensor 155 viewing vector with the t e r r a i n , can be computed a n a l y t i c a l l y . In the case with the r e a l DTM, the slopes and the i n t e r s e c t i o n points have to be estimated numerically. 4.6.2 Gaussian Terrain with Landsat-4 TM Sensor Parameters The Gaussian t e r r a i n model consists of several Gaussian h i l l s as i l l u s t r a t e d i n Figure 4.34. Here elevation has been shown d i r e c t l y proportional to image brightness. The t e r r a i n varies from 0 to 2,500 metres. A synthesized orthographic image using Landsat-4 TM sensor parameters of the same area i s shown i n Figure 4.35. Here, the sun i s at elevation 46° and azimuth 122°, corresponding to sun p o s i t i o n i n a TM scene. The image area i s about 15 km x 15 km. The synthesized image has a p i x e l s i z e of 30 metres x 30 metres. The image as seen by the sensor i s shown i n Figure 4.36. This distorted image i s synthesized by computing, a n a l y t i c a l l y , the i n t e r s e c t i o n of the sensor vector with the t e r r a i n and then the slopes at the point of i n t e r s e c t i o n . Since the worst case t e r r a i n d i s t o r t i o n i s to be studied, the area i s assumed to appear on the leftmost part of the imaged scene. Here, an elevation of 1 km i s displaced by about 130 metres. The di s t o r t e d image i s r e c t i f i e d by the system depicted i n Figure'2.6 (Section 2.6.2) and the r e c t i f i e d image i s shown i n Figure 4.37. The geometric accuracy of the r e c t i f i e d image i s i l l u s t r a t e d i n Figure 4.38. The t e r r a i n of one of the h i l l s i s shown i n Figure 4.38a and the i n t e n s i t y of the same area i n the orthographic image i s plotted i n Figure 4.38b. Figures 4.38c and d are i d e n t i c a l except for/minor radiometric errors. FIGURE 4.34 GAUSSIAN TERRAIN FIGURE 4 . 3 5 ORTHOGRAPHIC IMAGE (GAUSSIAN TERRAIN) 158 FIGURE 4.36 DISTORTED IMAGE (GAUSSIAN TERRAIN) FIGURE 4.37 RECTIFIED IMAGE (GAUSSIAN TERRAIN) 160 Height Intensity FIGURE 4.38 CROSS-SECTIONAL PLOT OF A HILL FIGURE 4 . 3 9 ERROR IMAGE (GAUSSIAN TERRAIN) 162 The r e l i e f displacement can e a s i l y be seen by comparing Figures 4.38b and c. The d i f f e r e n c e between the orthographic image and the r e c t i f i e d image shows the i n t e n s i t y of each Gaussian h i l l i s accurate to within 0.5 l e v e l RMS out of 256. The d i f f e r e n c e , accentuated by a f a c t o r of 80, i s shown i n Figure 4.39. This d i f f e r e n c e i s mostly due to an i n t e r p o l a t i o n such as ringing (Gibb's phenomenon) at the shadow boundaries. This simple experiment demonstrates the successful result obtained with the proposed r e c t i f i c a t i o n algorithm with an a n a l y t i c a l t e r r a i n model. 4.6.3 Real Terrain with Landsat-4 TM Sensor Parameters A regular gridded DTM of the S t . Mary Lake area i n B.C, Canada i s used with a s p a t i a l r e s o l u t i o n of 30 metres x 30 metres. The DTM covers an area of 15 km x 15 km. Figure 4.40 shows the t e r r a i n and Figure 4.41 shows the orthographic image. The t e r r a i n ranges from 950 metres i n the v a l l e y to 2,860 metres at the peaks. Again, to obtain the worst case t e r r a i n d i s t o r t i o n i n the d i s t o r t e d image, the area i s assumed to be at the leftmost part of the image. The d i s t o r t e d image i s shown i n Figure 4.42. The sun p o s i t i o n i s the same as i n the Gaussian t e r r a i n case, i . e . , from the southwest. (A more accurate sense of the t e r r a i n can be viewed by turning the picture, upside down since human v i s u a l i n t e r p r e t a t i o n i s more accustomed to having the i l l u m i n a t i o n from the top l e f t corner.) In synthesizing the d i s t o r t e d image, the slope at a point i n each 163 FIGURE 4.40 DTM OF ST. MARY LAKE AREA 164 FIGURE 4.41 LANDSAT-4 TM ORTHOGRAPHIC IMAGE 165 FIGURE 4.42 LANDSAT-4 TM RELIEF DISTORTED IMAGE 166 d i r e c t i o n i s determined by f i t t i n g a quadratic curve through the three closest points i n that d i r e c t i o n . Newton's form and divided differences can r e a d i l y be applied here. Let the three closest points to t, where the slope i s required, be ( t ^ , F j ) > ( C2' F2^ a n d ( t3» w n e r e the ^'s a r e the DTM values. Then the slope at t i s : F'(t) = F 2 - F 1 + 0.5 ( F 3 - 2 F 2 + Fj) (2t - t } - t 2> which can be computed very e a s i l y . ' The r e c t i f i e d image i s shown i n Figure 4.43. The geometric accuracy can best be i l l u s t r a t e d by Figure 4.44 which i s a plot of the displacement vectors as a r e s u l t of c o r r e l a t i o n between the o r i g i n a l orthographic image and the r e c t i f i e d image for 100 points. The mean error i s only 0.023 p i x e l and the standard deviation i s 0.077 p i x e l . The maximum error i s only 0.22 p i x e l . Hence subpixel geometric r e c t i f i c a t i o n accuracy has been achieved with the DTM. The RMS d i f f e r e n c e between the orthographic and the r e c t i f i e d images i s about 6 levels RMS out of 256. This radiometric d i f f e r e n c e i s caused by the fact that the i n t e n s i t y function of the synthesized d i s t o r t e d image i s not band l i m i t e d . The d i f f e r e n c e , accentuated by a f a c t o r of 20, i s shown i n Figure 4.45. 4.6.4 Real T e r r a i n with SPOT PLA Sensor Parameters The off nadir pointing angle f o r the SPOT sensor can be as high as 26*. A t e r r a i n height of 1,000 metres i s displaced by about 700 metres at 167 FIGURE 4.43 LANSAT-4 TM RECTIFIED IMAGE FIGURE 4.44 GEOMETRIC RECTIFICATION ACCURACY (LANDSAT-4 TM) FIGURE 4.45 ERROR IMAGE (LANDSAT-4 TM) Note: This results from both geometric and radiometric errors. 170 the far edge of the scene. The same regularly grldded DTM as In the Landsat-4 TM case Is used, but the spatial resolution i s now assumed to be 20 metres x 20 metres. Also, the altitude is multiplied by 2/3 to retain the same slopes. A r e a l i s t i c sun angle is taken to be at elevation 68° and azimuth 104°. The orthographic image is shown i n Figure 4.46 and the distorted image is shown i n Figure 4.47. It can be seen from Figure 4*47 how severely the mountain ranges are distorted. Moreover, since the sensor and the sun are at opposite sides of the image, the brighter sides of the slopes which face the sun are much narrower than the darker sides of the slopes which face away from the sun. This w i l l undoubtedly introduce radiometric errors in the rectified image as shown i n Figure 4.48. However, the displacement vectors for 100 points in Figure 4.49 show that the mean error is only 0.028 pixel and the standard deviation i s 0.07 pixel. Also, the maximum error is 0.42 pixel. Hence, subpixel geometric rectification accuracy can also be achieved in SPOT imagery for the PLA sensor when operating in the side looking mode. The RMS difference between the orthographic and the rectified images i s again about 6 levels out of 256 due to reasons already mentioned for the Landsat-4 TM case. The difference, accentuated by a factor of 20, is shown in Figure 4.50. 4.7 Experiments with Landsat-2 Imagery The Landsat-2 scene (I.D.2921-18025, imaged on 31st July 1977) used i s an image of an area over Vancouver, B.C., Canada. The image is shown in Figure 4.51. 171 ~Z FIGURE 4.46 SPOT ORTHOGRAPHIC IMAGE 172 FIGURE A.47 SPOT RELIEF DISTORTED IMAGE 173 FIGURE 4.48 SPOT RECTIFIED IMAGE FIGURE 4.49 GEOMETRIC RECTIFICATION ACCURACY (SPOT) 175 FIGURE 4.50 ERROR IMAGE (SPOT) Note: This r e s u l t s from both geometric and radiometric errors. FIGURE 4.51 LANDSAT-2 MSS RAW IMAGE 177 Results of the navigation procedure are shown i n Figure 4.52 and Table 4.4. In the fi g u r e , the RMS residual error stays r e l a t i v e l y constant at about 70 metres. The error i s mainly caused by: (a) Marking a target l i n e , p i x e l l o c a t i o n i n the image from which a p i x e l corresponds to 58 metres x 79 metres. In the process the maximum error i s estimated to be 1 p i x e l i n the along and across scan d i r e c t i o n . Assuming uniform d i s t r i b u t i o n , the RMS error i s : ,\2 /_\2 c, =11 — ] + | — | metres = 57 metres 1 11 bl 1/3 (b) Marking of target l o c a t i o n on a mapsheet. The mapsheets used are 1:50,000 sc a l e . A marking error of 1 mm corresponds to 50 metres. Assuming t h i s i s the maximum error, then the RMS error i s : °2 =/[ j^j + ^J^j Inet:res = 21 metres (c) Map accuracy. It i s not known what class the mapsheets belong to . Assume that the best maps (Class Al) are used, then the RMS error i n each d i r e c t i o n i s 16 metres. Denote t h i s by a^. Therefore the t o t a l error budget i s : a =./o\2 + + a 2 = 72 metres which i s close to the experimental r e s i d u a l error. RMS RESIDUAL ERROR METRES 80 ' 60 « 40 •• 20 « 0 8 10 12 14 16 NO. OF GCPs FIGURE 4.52 RESIDUAL ERROR VERSUS NUMBER OF GCPs TABLE 4.4 ATTITUDE AT 9 POINTS OVER THE IMAGE (NOTE: ALL VALUES IN DEGREES) ROLL -0.309 -0.315 -0.319 -0.323 -0.325 -0.326 -0.325 -0.322 -0.319 PITCH 0.124 0.120 0.116 0.120 0.109 0.107 0.107 0.111 0.119 YAW -0.516 -0.517 -0.518 -0.516 -0.514 -0.511 -0.505 -0.500 -0.494 179 The rectified scene is on a UTM projection with the north axis aligned with the vert i c a l . The angle of rotation (from raw image to O rectified image) is about 14 . The rectified image is .shown in Figure 4.53. The pixel spacing i s 50 metres and the scale i s 1:1,000,000. This choice of scaling proves to be very convenient i f a subarea of the image i s required and on a scale of l:l,000,000n or 1:1,000,000/n where n i s an integer. An example of using this rectified image on another scale w i l l be given i n the next chapter. 180 FIGURE 4.53 RECTIFIED LANDSAT-2 MSS IMAGE 181 CHAPTER 5 SYNTHETIC APERTURE RADAR 5.1 Introduction R e l i e f displacement i s a major factor i n synthetic aperture radar (SAR) imagery. It i s necessary to understand how SAR works and how the SAR signals are processed before the r e c t i f i c a t i o n remapping can be derived and r e l i e f displacement correction can be formulated. The theory of SAR can be found i n references [KOVA, MDA76, MDA77a, MDA77b and HOVA79] and i s summarized i n Appendix F. References [BENN78, BENN79, CUMM79, MDA79 and BENS80] show how the SAR signals are processed. This again i s described i n Appendix F. The basic properties of SAR and the range and azimuth resolutions attainable are summarized as follows. A SAR system sends out f i n i t e pulses and receives t h e i r echoes at regular i n t e r v a l s ; meanwhile the s a t e l l i t e carrying the SAR advances i n p o s i t i o n between pulses. The main c h a r a c t e r i s t i c s of the processed SAR imagery i s that i t has high azimuth (along track d i r e c t i o n ) r e s o l u t i o n which i s s a t e l l i t e a l t i t u d e independent, and that t h i s r e s o l u t i o n i s numerically proportional to the antenna aperture ( i . e . , the smaller the antenna aperture, the higher the r e s o l u t i o n ) . By v i r t u e of the s a t e l l i t e motion i n the azimuth d i r e c t i o n , the return s i g n a l of a target i s a chirp ( l i n e a r frequency modulated) signal i n t h i s d i r e c t i o n . By matched f i l t e r i n g t h i s return with i t s r e p l i c a , the azimuth resolution p obtained i s 182 P = D/2 (5.1) a where D i s the antenna aperture. High r e s o l u t i o n i n range (across track d i r e c t i o n ) can be obtained by transmitting a chirp sign a l and matched f i l t e r i n g the return by i t s r e p l i c a . The range r e s o l u t i o n thus obtained i s P r = 0.886 c / (2 KT) (5.2) where c i s the v e l o c i t y of propagation, K i s the l i n e a r frequency modulation (FM) rate and T i s the radar pulse duration. With 4 looks (Appendix F) and sidelobe suppression the azimuth and range r e s o l u t i o n f o r Seasat imagery i s p = 25 metres a p = 25 metres r Without matched f i l t e r i n g , the azimuth res o l u t i o n i s q = X R / D (5.3) a and the range r e s o l u t i o n i s q r = c T/ (2R) (5.4) where X i s the wavelength of operation and R i s the target range. For Seasat they correspond to q = 18.7 km a and q = 5.1 km r 183 These re s u l t s r e a d i l y show that the matched f i l t e r i n g i s important to obtain high resolutions i n both azimuth and range. The purposes of t h i s chapter are to: (a) Develop the r e c t i f i c a t i o n transformation which depends on how the SAR data i s processed. (b) Investigate the geometric d i s t o r t i o n s i n which r e l i e f displacement i s a major problem. (c) Develop o r b i t refinement methods. In SAR imagery, a t t i t u d e angles are not parameters i n the r e c t i f i c a t i o n transformation and hence these angles do not need to be r e f i n e d . (d) Perform r e c t i f i c a t i o n experiments with Seasat data. 5.2 R e c t i f i c a t i o n Transformation The MDA Seasat Processor described i n Appendix F i s capable of producing a high r e s o l u t i o n SAR image whose coordinates are slant range and azimuth. When converted to ground range, the resolution i s 25 X 25 2 metres^. To r e c t i f y t h i s to a map pro j e c t i o n (e.g., UTM), the geometric corrections must take i n t o account s a t e l l i t e ephemeris data, Doppler frequency to which the input data has been compressed (as defined by the Seasat processor), t e r r a i n , earth curvature and earth r o t a t i o n . These geometric corrections together with the output map projection parameters, define the transformation l i n k i n g the input image to the output image. The transformation i s deta i l e d i n Appendix G. 184 It should be noted that the slant range s is not linearly related to ground range g. They are related, from Figure 5.1, by: g = R 6 where cos 6 = (H 2 + R2 - s 2^ / (2HR) A plot of Ag/As as a function of range for a typical Seasat image is shown in Figure 5.2. An important result of the transformation i s that i t is independent of.the s a t e l l i t e attitude angles. These attitude angles only affect the positioning of the SAR beam which can be determined by an independent method. However, the attitude angles and orbit must be determined to within a tolerance limit for focusing purposes. 5.3 Terrain Induced Distortions 5.3.1 Relief Displacement Relief displacement is the major geometric distortion in radar imagery. It has the following attributes: (a) A target migrates towards the radar by an amount approximately equal to the product of i t s height and the tangent of the depression angle (see Figure 5.3). Hence, the displacement effect decreases with depression angle. This can also be explained by the fact that for a larger depression angle the slant range difference between the top and foot of an elevated target is 135 FIGURE 5.1 GEOMETRIC RELATION BETWEEN SLANT RANGE AND GROUND RANGE FIGURE 5.2 PLOT OF Ag/As VERSUS RANGE FIGURE 5.3 RELIEF DISPLACEMENT IN RADAR IMAGERY 188 g r e a t e r t h a n t h e d i f f e r e n c e f o r a s m a l l e r d e p r e s s i o n a n g l e . (b) S l o p e s f a c i n g t h e r a d a r appear t o be s h o r t e r i n t h e image and s l o p e s f a c i n g away f r o m t h e r a d a r appear t o be l o n g e r . T h i s i s r e f e r r e d t o as s l o p e f o r e s h o r t e n i n g . ( c ) Radar l a y o v e r o c c u r s when t a r g e t s f u r t h e r away on t h e ground appear t o be c l o s e r t o t h e r a d a r t h a n n e a r e r t a r g e t s . T h i s o c c u r s when t h e s l o p e exceeds a l i m i t d e f i n e d by t h e d e p r e s s i o n a n g l e . H ence, i n v e r y mountainous a r e a s , mountain t o p s and v a l l e y s can have t h e same p i x e l number ( i . e . , same s l a n t range) i n t h e image. I n r a d a r i m a g e r y , s a m p l i n g i s done at e q u a l i n c r e m e n t s i n s l a n t r a n g e . Hence, i n t h e remapping p r o c e s s , i t i s p r e f e r a b l e t o d e f i n e r e l i e f d i s p l a c e m e n t i n s l a n t range u n i t s . From F i g u r e 5.3, t h e d i s p l a c e m e n t i s g i v e n by Ad = d' - d. The v a l u e of d' i s known t o be a l i n e a r f u n c t i o n of p i x e l number, and d can be d e t e r m i n e d f r o m t h e geometry o f t h e f i g u r e : d =JR 2 + H 2 - 2RH c o s 6 (5.5) where cos0 = [ H 2 .+ 2(R+z)H - ( d ' ) 2 ] / [2(R+z)H] (5.6) Hence Ad can be e x p r e s s e d i n c l o s e d f o r m by f i n d i n g t h e d i f f e r e n c e of d' - d. S i n c e z << R, an a p p r o x i m a t i o n o f d can be f o r m u l a t e d as f o l l o w s . From E q u a t i o n s 5.5 and 5.6 i t can be shown t h a t Ad = d' - d A ; z (R - H cos 0 ) / d 189 From the geometry of the figure the above equation can be written as: where 3 is the radar look angle of the vertical and the negative sign shows that displacement is towards the radar. Using nominal values of Seasat orbit and look angle 8 » the displacement and the c r i t i c a l slope angle when layover occurs are shown in Figure 5.4. Relief displacement i s a high frequency geometric distortion and can be removed in the f i r s t processing pass as w i l l be demonstrated later. 5.3.2 Importance of Compression to Zero Doppler In the azimuth compression process, the matched filtered output i s placed at location of zero Doppler. This has significant impact on the remapping process for the case that the antenna illumination i s not broadside. The advantage of compression to zero Doppler i s that relief displacement i s always only along scan. If compression i s to nonzero Doppler, then relief displacement also has a nonzero along scan component. For example, the Seasat antenna i s t i l t e d at about 2°. If compression i s performed to the Doppler frequency along the beam centre, then a 1,000 metre high target w i l l be displaced by 2,300 sin 2° = 80 metres in the azimuth direction, assuming that each metre in elevation i s displaced by 2.3 metres on the ground. (5.7) 190 RELIEF DISPLACEMENT COEFFICIENT GROUND SLANT RANGE | RANGE 22.0 - .905 i t " 20 30 40 GROUND RANGE IN KM FROM FIRST PIXEL FIGURE 5.4 RELIEF DISPLACEMENT AND CRITICAL ANGLE NOTE: RELIEF DISPLACEMENT COEFFICIENT IS THE RATIO OF DISPLACEMENT TO TARGET ELEVATION 191 5.3.3 Target Defocuslng Another problem with t e r r a i n i s that the t r a j e c t o r y of an elevated point target i s incorrect i n the azimuth matched f i l t e r i n g process, hence, defocusing the target. The following examines how the tr a j e c t o r y deviates from zero elevation and hence assesses the amount of defocusing. For Seasat imagery, the worst d i s t o r t i o n i n slant range r i s shown in Figure 5.5 as Ar. This occurs at the t i p of the i l l u m i n a t i o n beam. Using the Seasat parameters as shown and nominal o r b i t parameters, i t can be shown that Ar < 0.1 metre per km i n elevation. This amount of Ar has no s i g n i f i c a n t e ffect on range c e l l migration and hence, produces no s i g n i f i c a n t target b l u r r i n g . 5.3.4 Intensity Migration Intensity migration occurs when the same area i s imaged i n more than one pass. It i s caused by a change i n incidence angle to the same l o c a t i o n due to the r e l a t i v e p o s i t i o n of the or b i t i n the passes. The change i n incident angle i s i l l u s t r a t e d i n Figure 5.6. has Intensity an adverse migration, although not s t r i c t l y e f f e c t on the r e g i s t r a t i o n between a geometric d i s t o r t i o n , two radar scenes. 192 TIME n FIGURE 5 5 PREDICTED AND ACTUAL TARGET TRAJECTORIES OF AN ELEVATED TARGET 193 SATELLITE POSITION I NOTE: POINT A IS OF MAXIMUM REFLECTIVITY TO . SATELLITE POSITION I POINT B IS OF MAXIMUM REFLECTIVITY TO SATELLITE POSITION II FIGURE 5.6 INTENSITY MIGRATION IN RADAR IMAGERY 194 5.4 Experiments with Seasat Imagery 5.4.1 Image Location Accuracy Study The purpose of the image l o c a t i o n accuracy study i s to determine the geometric accuracy i n SAR imagery produced by a t y p i c a l SAR processor (which i n t h i s case i s the MDA SAR processor) and hence drive a navigation method. The experiment i n the study consisted of the following procedures: (a) Select a number of c l e a r l y i d e n t i f i a b l e targets i n several Seasat scenes (each scene i s 40 km X 40 km). (b) From topographic maps, determine the l o c a t i o n and elevation of these targets. (c) Using only the location of the target within an image, as well as Goddard Space F l i g h t Centre supplied D e f i n i t i v e Orbit Record data [JFL77], predict the l o c a t i o n of each of these targets, correcting for each target's elevation. (d) Determine the means and standard deviations of the differences derived from steps (b) and ( c ) . The scenes selected were as follows: . - Ottawa, Ontario, Orbit 472, descending, - Pembroke, Ontario, Orbit 1218, descending, and - Vancouver Island, B.C., Orbit 193, ascending. The f i r s t two scenes were received at Shoe Cove, Newfoundland while the l a s t one was received at Goldstone, C a l i f o r n i a . The three chosen scenes vary from gently r o l l i n g t e r r a i n (e.g., Pembroke scene where the t e r r a i n 195 varies between 100 and 200 metres) to very rugged terrain (e.g., Vancouver Island scene where the terrain varies from 0 to 1,000 metres). The maps used were 1:50,000 Class Al maps from the National Topographic Series produced in Canada. Each map i s on a UTM grid. The accuracy of Al maps has already been discussed in Chapter 3. Experimental results are summarized in Table 5.1 which gives the along and across track errors. The table shows two types of errors: - absolute or mean error with bias in the along and across track directions; and - relative error after bias removal and expressed in terms of standard deviation. The across track bias lies below 54 metres for a l l three scenes. The along track bias varies drastically in the Pembroke scene from the other two scenes. In the Pembroke scene i t exceeds 6,000 metres. This substantial bias is believed to result from a timing error in the receiving station clock of approximately one second. The relative error results from a variety of sources : (a) Marking a target line and pixel coordinate in the image. A pixel and a line each corresponds to 12.5 metres on the ground. (b) Marking of target on a mapsheet. In a 1:50,000 scale mapsheet, a measurement error of 1/2 mm would correspond to 25 TABLE 5.1 ALONG TRACK AND ACROSS TRACK ERRORS Scene No. o f Targets Along t r a c k b i a s i n Metres Across t r a c k b i a s i n Metres Along t r a c k a i n metres E x p e r i m e n t a l / Budgeted Across t r a c k a i n metres E x p e r i m e n t a l / Budgeted Ottawa 17 47 35 23/30 38/27 Pembroke 19 6186 54 24/30 38/27 Vancouver Island 12 117 44 23/30 13/29 metres. Map accuracy. For Class A l 1:50,000 mapsheets, the s p a t i a l RMS error i s 15 metres i n easting and northing d i r e c t i o n s combined (equivalent to 15 metres i n along track and across track combined). Orbit ephemeris data e r r o r . The D e f i n i t i v e Orbit Record data supplied by the Goddard Space F l i g h t Centre contains s a t e l l i t e p o s i t i o n information once every minute. The error f o r each p o s i t i o n i s [JPL77]: a i n along track = 17 metres a i n across track = 10 metres a i n r a d i a l d i r e c t i o n = 10 metres For each scene, a quartic polynomial f i t to the s a t e l l i t e p o s i t i o n was used to define the s a t e l l i t e p o s i t i o n . This can be shown to give an error of less than 1.2 metres i n each d i r e c t i o n . Hence the ephemeris i n t e r p o l a t i o n error i s n e g l i g i b l e . R e l i e f displacement. This has e f f e c t only i n the across track d i r e c t i o n . As has been mentioned, o f o r elevation error i s 0.3 contour i n t e r v a l f or Class Al maps. On gently r o l l i n g areas, the contour i n t e r v a l i s 25 f e e t . On h i l l y areas, i t can be 50 feet or 100 f e e t . Then f o r gently r o l l i n g areas, a (Class Al map) = 2.29 metres. This corresponds to 5 metres i n ground displacement, since f o r every metre i n height, r e l i e f displacement i s approximately 2.3 metres on, the ground. Contour I n t e r v a l Interpolation e r r o r . This again has e f f e c t only i n the across track d i r e c t i o n . Usually a target does not 1 9 8 f a l l onto a contour line in the mapsheet and the target height has to be interpolated. An interpolation accuracy of 1/4 contour interval is estimated. Thus, for a contour interval of 25 feet and assuming uniform distribution i n interpolation error, o (Class Al map) = 1.1 metres. This corresponds to 3 metres in relief displacement. The budgeted overall a was calculated for each scene and the result i s also shown in Table 5.1. The contour intervals of various mapsheets have to be taken into account in the calculation. For the Ottawa and Pembroke scenes, the maps used were with both 25 and 50 foot contours; and for the Vancouver Island scene, the mapsheets were with 20 metre contour interval. The experimental errors agree f a i r l y well with the predicted values for the f i r s t two scenes and l i e within the predicted value for the third scene. 5.4.2 Orbit Refinement The transformation described in Appendix G is a function of earth parameters and orbit ephemeris, and is attitude independent. Results in Section 5.4.1 show that the orbit ephemeris has to be refined for precision rectification of SAR imagery. In Seasat SAR imagery, the orbit latitude $ and longitude X are expressed as quartic polynomials: * i (j) = E <j> t (5.8a) 1=0 199 X = E X. t 1 (5.8b) 1=0 1 For convenience t = 0 is chosen as the image scene centre. Results i n the image location accuracy study clearly show that by adding a bias to both <()Q and XQ , the resulting rectified image has a residual error in the along and across track directions of about 40 metres in which measurement error is included. Therefore, the true rectification error, free from map induced measurement error, is less than 40 metres in the two directions. An even simpler way is to translate the image with respect to the map projection grid lines. The worst translations occur in the Pembroke scene (about 6 km in the along track direction). The resulting image w i l l s t i l l have the accuracy cited as above. Although the two methods (adding biases <j) ^ and X Q and image translation) are not mathematically equivalent, i t can be shown that the 2 differences are negligible over a 40 X 40 km scene. Such a method of rectification of Seasat SAR imagery is possible due to the high degree of accuracy in tracking data supplied by the Goddard Space Flight Centre. The along track error and across track error are most probably due to timing errors and hence can be removed by an image translation i n the manner described above as shown in Figure 5.7. These plots do not exhibit any residual height error or 200 (A) OTTAWA SCENE (B) PEMBROKE SCENE \ (C) VANCOUVER ISLAND SCENE 50 m FIGURE 5.7 GCP ERROR VECTORS AFTER NAVIGATION 201 r e s i d u a l image r o t a t i o n , hence j u s t i f y i n g the t r a n s l a t i o n process. For SAR imagery i n general, o r b i t refinement can be done by a least squares f i l t e r i n g technique s i m i l a r to that presented i n Section 4.3 f o r scanner imagery. However, at t i t u d e angles should not be entered i n the state vector. For t y p i c a l SAR imagery the parameters are s a t e l l i t e l a t i t u d e , longitude and l o c a l o r b i t radius. 5.4.3 R e c t i f i c a t i o n with DTM The Vancouver Island scene used i n the image lo c a t i o n accuracy study i s also used as an i l l u s t r a t i v e example of how a DTM can be used i n the r e c t i f i c a t i o n . This scene i s chosen because i t i s imaged over a very rugged area and a DTM i s a v a i l a b l e . The DTM covers an area of 16 km X 20 km and has s p a t i a l r e s o l u t i o n of 60 metres X 60 metres. It was obtained from a triangulated network [PENC78 and HEIL78] on a 1:50,000 Class Al topographic map with a 20 metre contour i n t e r v a l . More data points were d i g i t i z e d i n steep areas than r e l a t i v e l y f l a t areas. About 5,000 points were manually d i g i t i z e d . The navigation process has been performed i n the image lo c a t i o n accuracy study. The remapping and i n t e r p o l a t i o n was performed i n three one-dimensional passes and the northing d i r e c t i o n i n the r e s u l t i n g image was to be aligned with the v e r t i c a l . The block size versus remapping accuracy i s shown i n Figure 5.8. It was chosen to be 100 p i x e l s X 100 202 0 100 200 300 400 SCO BLOCK SIZE IN PIXELS NOTE: PIXEL SIZE IS 12.5 METRES FIGURE 5.8 BLOCK SIZE VERSUS REMAPPING ACCURACY 203 l i n e s i n the f i r s t and second passes, r e s u l t i n g i n a 0.1 p i x e l RMS remapping e r r o r i n the f i n a l product and a l l p i x e l s l i e w i t h i n 0.5 p i x e l remapping e r r o r . The o r i g i n a l image i s shown i n Figure 5.9. The image appears to be "turned around" s i n c e i t was acquired I n an ascending pass. Figure 5.10 shows the r e c t i f i e d image w i t h the use of the DTM while Figure 5.11 shows the r e c t i f i e d image without using the DTM. The sharp boundaries i n Figure 5.10 show where the DTM ends. Note that i n the l a s t two f i g u r e s the north i s not a l i g n e d w i t h the northing a x i s due to a nonzero convergence angle [MALI73] i n the UTM p r o j e c t i o n . The r e c t i f i c a t i o n e r r o r , by measuring s e v e r a l c o n t r o l points on the lakes i n the mountainous areas, was found to be 35 metres RMS. T h i s i s w i t h i n the e r r o r budget c o n s i d e r i n g GCP marking e r r o r and map accuracy.. A comparison of the l a s t two f i g u r e s c l e a r l y demonstrates the advantage of using a DTM. Many areas i n the r e c t i f i e d image i n Figure 5.10 appear to be orthographic w h i l e i n Figure 5.11 the same areas appear to be s h i f t e d towards the radar. The "noisy" nature i n areas w i t h steep slope i s due to r e l i e f layover and slope f o r e s h o r t e n i n g e f f e c t s . A very narrow band of bright p i x e l s ( f a c i n g the radar) i n the raw image has to be mapped to a wide area. Furthermore, speckles are present w i t h i n t h i s narrow band. The s i n e f u n c t i o n i n t e r p o l a n t k e r n e l does not take slope i n t o account. The removal of r a d i o m e t r i c e f f e c t s due to slope i s l e f t as a f u r t h e r research t o p i c . 204 mnw i zo FIGURE 5 . 9 RAW SEASAT IMAGE FIGURE 5.10 RECTIFIED SEASAT IMAGE WITH DTM 206 UTM ZONE 4 3 0 0 6 0 1 0 4 4 0 0 0 0 s s s 00 m in IS <s ( 5 0 0 0 0 1 0 4 3 0 0 0 0 4 4 0 0 0 0 4 5 0 0 0 0 FIGURE 5.11 RECTIFIED SEASAT IMAGE WITHOUT DTM 207 An undesirable geometric error r e s u l t i n g from using the DTM i n a steep area i s i l l u s t r a t e d by Grant Lake (coordinates: 443,000 Easting and 5,385,000 Northing). Its potato shape appearance i n the raw image has changed to an L shape i n the r e c t i f i e d image. The protruding part corresponds to an area of extreme high slope of 38° as measured from the given topographic map. T h i s has f a r exceeded the c r i t i c a l slope angle shown i n Figure 5.4, thus r e l i e f layover has occurred. The r e c t i f i c a t i o n error i n t h i s part i s due to DTM inaccuracy. Experimental r e s u l t s show that the remapping error i s about 60 metres and t h i s corresponds to an elevation error of 26 metres (since the r e l i e f displacement i s about 2.3 metres per metre height). Indeed, some DTM values are found to d i f f e r from the topographic map by as much as 20 metres i n the h i l l s around Grant Lake. The DTM error i s due to the following f a c t o r s : (a) D i g i t i z a t i o n e r r o r . In a 1:50,000 map with 20 metre contour i n t e r v a l , 4 mm crosses 2 contour i n t e r v a l s at 38° slope. Hence, a 0.5 mm measurement error would correspond to 20 metres. (b) I n t e r p o l a t i o n e r r o r . The i n i t i a l DTM was set up i n a triangulated network and then interpolated to a regular DTM. In' a steep area, the i n i t i a l DTM must therefore be dense enough to r e t a i n accuracy i n the regular DTM. (c) Water l e v e l . It is.found that the width of the lake i n the image d i f f e r s from that i n the map. This may be caused by seasonal changes i n the water l e v e l . 208 5.5 Navigation by Automatic Focusing A key feature i n the r e c t i f i c a t i o n of SAR imagery i s the p o s s i b i l i t y of r e f i n i n g s a t e l l i t e t r a j e c t o r y data without using GCPs. Th i s method i s automatic focusing. In both range and azimuth, the accuracy of focusing i s dependent on the p r e c i s i o n to which the matched f i l t e r required f o r compression can be determined. In azimuth compression the FM rate i s a function of the earth spheroid and s a t e l l i t e o r b i t parameters. The uncertainties i n these parameters w i l l introduce an error i n FM rate and hence w i l l not bring a target into focus. An example of azimuth misfocus i s shown i n Figure 5.12 where the SAR data i s processed f i r s t with the correct FM rate and then with an FM rate o f f by 10%. Di s t o r t i o n s i n the FM rate with respect to various parameters have been analyzed for Seasat [MDA77] and are presented i n Table 5.2. For two geographic l o c a t i o n s . Focusing targets i n the imagery can thus provide a way to r e f i n e the parameters. One technique i s to process the imagery with the nominal value of the FM rate, so that d i s c r e t e strong targets can be i d e n t i f i e d . These targets are processed again by varying the FM rate aiming at minimizing the h a l f power width of the responses and maximizing the peak magnitudes of the selected targets. Another technique takes advantage of the fa c t that r e g i s t r a t i o n between azimuth looks, i n multi-look processing, i s very s e n s i t i v e to 209 GROUND RfiNGE (A) PROCESSING WITH CORRECT FM RATE GROUND RANGE (B) PROCESSING WITH FM RATE DISTORTED BY 10% FIGURE 5.12 EFFECT OF AZIMUTH MISFOCUSING ( C o u r t e s y o f Mr. P. Hasan and Mr. P. McConnell o f MDA) TABLE 5.2 SENSITIVITY OF FM RATE TO PARAMETER VARIATIONS IN SEASAT IMAGERY % Change i n FM Rate U n i t s T a r g e t l o c a t i o n 60°N 58°W 0°N 9 6 ° W T a r g e t l a t i t u d e 4.62 3.46 /deg O r b i t semi-major a x i s 20.52 x 1 0 - 6 15.62 x 1 0 - 6 /m O r b i t e c c e n t r i c i t y 121.73 0.0 S a t e ! 1 i t e a n g u l a r f r e q u e n c y 3250 3269 /(deg) sec O r b i t i n c l i n a t i o n 0.32 0.22 /deg E a r t h semi-major a x i s 8.85 x 1 0 - 6 13.85 x 1 0 - 6 /m 211 variation in the azimuth matched f i l t e r FM rate [BENN78]. Hence the FM rate can be tuned until the looks are registered. From the refined FM rate at different positions, the orbit trajectory parameters can be refined by, again, a recursive estimator as shown in Figure 5.13. Detailed study of this method is left for future research. 212 AZIMUTH/SLANT RANGE COORDINATES PARAMETERS FM RATE COMPUTED — FUNCTION FM RATE ~V MEASURED FM RATE USING AUTOMATIC FOCUSING UPDATED RECURSIVE FILTER PARAMETERS FIGURE 5.13 NAVIGATION LOOP USING AUTOMATIC FOCUSING 213 CHAPTER 6 APPLICATIONS 6.1 Introduction This chapter presents some applications of p r e c i s i o n r e c t i f i e d s a t e l l i t e imagery and an a p p l i c a t i o n of the r e c t i f i c a t i o n method to a future planetary mission. Applications of p r e c i s i o n r e c t i f i e d imagery presented here include multisensor i n t e g r a t i o n and image mosaicking. These two examples merely i l l u s t r a t e the p o t e n t i a l of r e c t i f y i n g imagery to a common data base, and are by no means exhaustive instances i n the applications area. The r e c t i f i c a t i o n method can be extended to a planned future NASA planetary exploration mission - Venus Or b i t i n g Imaging Radar. The lack of GCPs w i l l make the automatic focusing method promising i n the r e c t i f i c a t i o n process. 6.2 Multisensor Integration A Synthetic Aperture Radar (SAR) can produce very high r e s o l u t i o n imagery which contains f i n e d e t a i l s i n surface roughness. On. the contrary, Landsat imagery i s of lower r e s o l u t i o n but contains information about surface r e f l e c t i v i t y i n d i f f e r e n t s p e c t r a l bands. Since Landsat and SAR imagery complement one another i n t h i s way, t h e i r composite would contain very r i c h and useful information f o r image i n t e r p r e t a t i o n . 214 With the availability of spaceborne SAR imagery (Seasat), the production of Landsat/Seasat colour composites becomes feasible. Such products may be produced routinely by a future sa t e l l i t e system incorporating both an MSS and SAR on the same platform. In this thesis an example is shown by a Landsat/Seasat colour composite image of the Ottawa area in Canada, precision rectified to UTM coordinates. Basically, both the Landsat and Seasat images are registered to UTM coordinates, the procedures of which have been described in previous sections. Then the intensity of the black and white Seasat image is used to modulate the false colour Landsat image (made from bands 4, 5 and 7). The resulting image can be viewed on a video display or written to film by a film recorder. Figures 6.1 to 6.3 are pictures taken directly from a 512 x 512 d i g i t a l display. The dynamic range of the display is 8 bits. Each 2 pixel is 25 x 25 metre . The Seasat image is shown in Figure 6.1 and the Landsat image in Figure 6.2. The Seasat image is used to modulate the Landsat image. The modulation is justified by the fact that image intensity IQ is related to the surface albedo p by a multiplicative rule: I 0 = K p where K is a function of imaging geometry, atmospheric parameters and wave length etc. The picture in Figure 6.3 is produced in the following modulation way: FIGURE 6.1 SEASAT COMPONENT IN LANDSAT/SEASAT COMPOSITE 216 FIGURE 6.2 LANDSAT COMPONENT IN LANDSAT/SEASAT COMPOSITE I FIGURE 6.3 LANDSAT/SEASAT COMPOSITE (12 KM X 12 KM) OF THE OTTAWA AREA 218 Band 1 n = /Band n of Landsat • Seasat where n = 4 , 5 , 7 . Therefore, In each case 3 intermediate products were made. The colour composites were made with the following colour assignment i n each case: Band' 4 blue Band' 5 green Band' 7 red Also for each Band' n, i t s histogram i s used as a guide to set up an i n t e n s i t y transfer function so that the maximum dynamic range can be u t i l i z e d . The geometric mean i s used f o r the following reason. Let two random variables be of the same mean Their product then has a mean < MQ while the geometric mean has a mean of ^ - MQ as i l l u s t r a t e d i n Figure 6.4. The geometric mean composite then has the same "colour balance" as the Landsat. Note that the geometric mean i s equivalent to summation i n the logarithm domain: log(composite) = 0 . 5 [log(Seasat) + log(Landsat)] Figures 6.1 and 6.2 give a comparison between the Seasat and Landsat images. The Landsat image i s of low re s o l u t i o n but contains good surface r e f l e c t i v i t y information since the Landsat sensors operate i n a region of sp e c t r a l bands. On the other hand, the Seasat image i s of HISTOGRAM OF VARIABLE X FIGURE 6.4 HISTOGRAM OF PRODUCT OF TWO RANDOM VARIABLES 220 high r e s o l u t i o n and contains good surface texture information 1. By combining the two, as shown i n Figure 6.3, the advantages of both sets of sensors are integrated into one image. The colour composite shows bridges, roads, highways and t e r r a i n prominently whereas these structures are not as clear i n the Landsat image. In addition, the colour composite has s p e c t r a l reflectance content not present i n the Seasat image alone. To s u i t various applications such as f o r e s t r y , geology or urban planning, the colour assignment and the dynamic range of i n d i v i d u a l components may have to be varied. In p a r t i c u l a r several places of i n t e r e s t are (see Figure 6.5 for region number): Region 1: The inner structure of the square can be seen here and the square as a whole i s c l e a r l y noticeable i n the compositej but t h i s i s not very noticeable i n the Seasat image. While i t i s also noticeable i n the Landsat image, i t does not show high r e s o l u t i o n . Region 2 : The Ottawa River here i s not very clear i n the Seasat component. Region 3: The target here i s not v i s i b l e i n the composite. This i s expected since t h i s target only appears i n the Seasat image and not i n the "Landsat image. Region A: The Seasat image saturates i n t h i s area (due to corner r e f l e c t o r s formed by the downtown Ottawa area) but the Landsat image does not. The colour composite also saturates i n t h i s area because a l i n e a r s t r e t c h function 1 One example i s that Seasat imagery shows geological f a u l t l i n e s prominently. REGION 1 REGION 4 REGION 2 FIGURE 6.5 REGIONS OF INTEREST IN THE COMPOSITE 222 with c l i p p i n g at high i n t e n s i t y l e v e l s i s used i n producing the composite. Saturation can be prevented i n the composite by using the cumulative d i s t r i b u t i o n function as the transfer function. General: As has been mentioned, the composite shows bridges, roads, highways and downtown Ottawa very prominently. "Disappearance" of targets i n the geometric mean composite (Region 3) can be avoided i f the composite i s done by: where a^ ( i = 1 to 4) are su i t a b l y chosen c o e f f i c i e n t s . In other words, a l i n e a r s t r e t c h i s applied to each image f i r s t before taking the product. By examining the terms of the expression, the composite contains m u l t i p l i c a t i v e as well as additive terms. A composite of a larger area around the Ottawa Valley i s shown i n 2 Figure 6.6. Each p i x e l covers an area of 12.5 x 12.5 metres i n t h i s composite. Another example i s the Vancouver Landsat/Seasat composite as shown i n Figure 6.7 (Seasat) and Figure 6.8 (composite). The r e c t i f i e d Landsat image has been shown i n Figure 4.53. The area corresponds to NTS mapsheet number 92G/3 which i s shown i n Figure 6.9. This example demonstrates the benefits of r e c t i f y i n g s a t e l l i t e imagery to a common data base, e.g., to a conventional mapsheet numbering system i n the following 223 FIGURE 6.6 LANDSAT/SEASAT COMPOSITE (40 KM X 40 KM) OF THE OTTAWA AREA 226 FIGURE 6.9 NTS MAPSHEET 92G/3 (Surveys and Mapping Branch, Department o f Energy, Mines and R e s o u r c e s , Ottawa, Canada.) 227 sense. The r e c t i f i e d image can e a s i l y be overlaid to an e x i s t i n g mapsheet. In addition, an "image map" i n which topographic information from a mapsheet i s overlaid onto the image can be made. 6.3 Image Mosaicking An a p p l i c a t i o n of r e c t i f i e d imagery i n neighbouring passes to the same data base i s the a b i l i t y to form mosaics. This subsection discusses d i f f e r e n t methods of mosaicking and presents an example of image mosaicking. 6.3.1 Linear Interpolation V i s u a l d i s c o n t i n u i t i e s w i l l occur i f two images are. joined together without some kind of smoothing operation at the j o i n t . This i s because the l o c a l histograms for the two images at the j o i n t are not balanced even though the histograms f o r both images may be balanced. Smoothing at the j o i n t can be done by l i n e a r i n t e r p o l a t i o n of the i n t e n s i t i e s i n the overlap region. Let the image i n t e n s i t i e s be denoted by 1^ and and the width, of the overlap area denoted by L as shown i n Figure 6.10. Then the i n t e n s i t y i n the overlap area i s given by: I (x) = [(L - x) I.(x) + x I 0 ( x ) ] / L (6.1) c i i where x i s defined i n the f i g u r e . Equation 6.1 produces a gradual t r a n s i t i o n or ramp e f f e c t . —*• x _ i 1 S , 2_ INTENSITY IN THE OVERLAP AREA GIVEN BY; I c = Ix (x) + J L I 2 (x) FIGURE 6.10 RADIOMETRIC MOSAICKING BY LINEAR INTERPOLATION 229 The disadvantage of t h i s method i s image degradation. The ramp slope determines how many points are affected by the smoothing process. A more gradual slope conceals the d i s c o n t i n u i t y better, but i t also has the adverse e f f e c t of degrading information content i n the v i c i n i t y of the j o i n t . Another disadvantage i s that features that appear i n one image but not the other (e.g., clouds) seem to fade away miraculously. However, t h i s method i s fast and easy to implement and gives good v i s u a l r e s u l t s at the boundary. 6.3.2 Feathering In t h i s method, an "optimal" seam point i s located for each image l i n e , hence the seam edge zigzags across the mosaicked image. This method i s proposed by Milgram [MILG75] and consists of three steps. The f i r s t step i s to equalize the gray l e v e l s t a t i s t i c s (mean and standard deviation) between the two images. The second step i s to locate the seam point within the overlap region. In the. mosaic, one image f i t s to the right of the seam point and the other image to the l e f t of the seam point. The t h i r d step i s to smooth across the seam using a ramp function. To equalize the gray l e v e l s t a t i s t i c s i n the overlap area, those regions containing features that occur i n the one image but not the other should be excluded i n computing the gray l e v e l transform. For example, clouds, snow or ground water present i n one image but not the other should be excluded. The next step i s to choose a seam point for each mosaic image 230 l i n e . This choice i s based on the minimum l o c a l d i f f e r e n c e i n the two images i n the overlap area. The seam point s e l e c t i o n f o r each mosaic image l i n e i s based on the minimum l o c a l d i f f e r e n c e i n the overlap area. A d i f f e r e n c e i s formed at each p i x e l l o c a t i o n i n the overlap area, and the p i x e l at which the minimum di f f e r e n c e occurs i s chosen as the seam point. The f i n a l step i s seam smoothing. This i s again l i n e a r i n t e r p o l a t i o n with a chosen ramp function width. Using t h i s method of seam point s e l e c t i o n , one finds a succession of points whose h o r i z o n t a l positions are unrelated to one another. This random positioning of the seam i s known as "feathering", since i t seems that one image i s "feathered" i n t o the other i n the mosaic. This method has the advantage of being able to track around clouds or other features present i n one but not i n the other. This method has the disadvantage of introducing d i s c o n t i n u i t i e s i n the v e r t i c a l d i r e c t i o n . A refinement would be to r e s t r i c t the range of candidate seam points depending on the magnitude of the minimum l o c a l d i f f e r e n c e on the previous l i n e . Milgram also proposes that the range of candidate seam points should be r e s t r i c t e d from l i n e to l i n e . An extension to r e f i n e the seam point locations has also been proposed by Milgram [MILG77]. He defines a cost function as the sum of minimum l o c a l differences f o r a l l l i n e s . The best seam through the mosaic i s the least cost path. 231 This algorithm suffers from the fact that the operator has to input the value of overlap width where the minimum l o c a l d i f f e r e n c e i s to be obtained, an expression for range of a candidate seam point and the smoothing ramp function width. An added disadvantage i s that p i x e l i n t e n s i t i e s of one scene have to be adjusted f i r s t ( f o r gray l e v e l balance) and t h i s sometimes may not be desirable.. 6.3.3 F i l t e r i n g J A f i l t e r i n g method has been employed by the USGS for several years. The steps for a two-frame mosaic, depicted i n Figure 6.11a,..are: 1. Mosaic the two images with a v e r t i c a l seam. 2. Low pass f i l t e r the mosaicked image. 3. High pass f i l t e r the o r i g i n a l two images and mosaick them with the same v e r t i c a l seam as i n Step 1. Size of the high pass f i l t e r should be equal to that of the low pass f i l t e r . 4. Add the outputs from Steps 2 and 3. A geometric i n t e r p r e t a t i o n of t h i s method i s shown i n Figure 6.11b. The low pass f i l t e r gives the effect of a ramp function and the high pass f i l t e r restores the high frequency components which are removed by the low pass f i l t e r . T his technique y i e l d s a mosaic s i m i l a r to that of l i n e a r i n t e r p o l a t i o n , but i s computationally more expensive. METHOD I M A C E 1 I M A O E 2 H I G H P A S S F I L T E R E D H I G H P A S S F I L T E R E D M O S A I C W I T H V E R T I C A L S E A M F I N A L P R O D U C T M O S A I C W I T H V E R T I C A L S E A M L O W F I L T I P A S S E R E D GEOMETRIC INTERPRETATION I M A G E 1 I M A G E 2 M O S A I C O F L O W P A S S F I L T E R M O S A I C O F H I G H P A S S F I L T E R ' S U M O F L A S T T W O P R O F I L E S FIGURE 6.11 FILTERING METHOD OF IMAGE MOSAICKING 233 6.3.4 Experiments The linear interpolation and f i l t e r i n g methods are considered in the experiment. The feathering method is not considered because of i t s "feathering" artifact. The test data used were two scenes imaged over an ice area of Cambridge Bay, Canada. In order to highlight the mosaicking effect, one image has been adjusted to have a higher mean than the other. The result i s shown in Figure 6.12 where the two scenes were mosaicked with a vertical seam. Figure 6.13 shows mosaicking by linear interpolation while Figure 6.14 shows mosaicking by f i l t e r i n g . The mosaic frame is 240 pixels wide and the overlap area is 60 pixels in the middle of the mosaic frame. The mean of the left frame is 20 intensity levels (out of 256) higher than that of the right frame. The results show that the smoothing effect of these two methods is quite similar (ignoring the effect of the f i l t e r width). However, in the f i l t e r i n g method there appears to be a vertical seam through the lake in the overlap area. This is caused by the difference in contrast of land to lake in the frames. Such difference in contrast does not have any effect in the linear interpolation method. For this reason, linear interpolation is recommended and futhermore, i t is computationally faster than the f i l t e r i n g method. 234 235 FIGURE 6.13 MOSAIC BY LINEAR INTERPOLATION FIGURE 6.H MOSAIC BY FILTERING 237 6.4 A Future Planetary Mission 6.4.1 Mission Description The Venus Or b i t i n g Imaging Radar (VOIR) due to be launched i n the future under the auspices of NASA/JPL w i l l carry a SAR System on board to chart the planet's surface. The s a t e l l i t e and i t s SAR system are designed to have the parameters presented i n Table 6.1. Complementary to the SAR system, the s a t e l l i t e w i l l also carry an altimeter on board. In the r e c t i f i c a t i o n of SAR imagery obtained from the mission, s p e c i a l attention must be paid to: - the navigation phase since no GCPs are a v a i l a b l e , - the remapping transform which should also include the e f f e c t of atmospheric r e f r a c t i o n i n the thick Venus atmosphere, and the use of altimeter data. 6.4.2 Navigation Due to the lack of GCPs, navigation w i l l have to be performed '• by the automatic focusing method. The method has been outlined i n Section 5.5. The s e n s i t i v i t y of FM rate to various parameters i s shown i n Table 6.1. TABLE 6.1 SENSITIVITY OF FM RATE TO PARAMETER VARIATIONS IN VOIR IMAGERY % Change i n FM Rate U n i t s T a r g e t L o c a t i o n 0°N 71°N T a r g e t L a t i t u d e 0.51 x 10"* 0.35 x 10"* /deg. O r b i t Semi-major A x i s 0.016 0.016 /m O r b i t E c c e n t r i c i t y 2.0 2.0 S a t e l l i t e A n g u l a r Frequency 3117 3055 / ( d £ 2 ) KsecJ O r b i t I n c l i n a t i o n 0.055 0.095 /deg Venus Semi-Major A x i s 0.017 0.017 /m Venus R o t a t i o n Rate 15.7 29.1 /deg 239 6.4.3 Venus Atmospheric Refraction Since an electromagnetic wave r e f r a c t s on entering Venusian atmosphere at nonnormal incidence, the SAR beam geometry departs from the l i n e a r . This e f f e c t changes the l o c a t i o n of the swath by about 0.5 km, and causes the low r e s o l u t i o n swath to be about 16 metres narrower than i n the absence of r e f r a c t i o n . T h e o r e t i c a l l y , the effect w i l l contribute to inaccuracies i n the o v e r a l l target l o c a t i o n equations, and can be taken i n t o account by introducing a s h i f t y , z g (see Figure 6.15) i n y and z In the s a t e l l i t e coordinate system y points i n the s s • s along track d i r e c t i o n and z points from the s a t e l l i t e to the planet s cent re. A d e s c r i p t i o n of the technique for computing these deviations follows. The index of r e f r a c t i o n i s calculated from the temperature, pressure density, and composition as a function of a l t i t u d e . The dependence of index of r e f r a c t i o n on height i s shown i n Figure 6.16 [MCCA80], where a s t a t i c s p h e r i c a l symmetric multilayer model i s used f o r the Venusian atmosphere. To trace a ray path from the s a t e l l i t e to the surface, the atmosphere i s assumed to consist of many concentric s h e l l s of f i n i t e thickness and constant r e f r a c t i o n index. When the ray i n t e r s e c t s a boundary between two l a y e r s , S n e l l ' s law i s applied to calculate the change i n d i r e c t i o n . The ray then propagates i n a straight l i n e to the next boundary, where the process repeats. The mathematical treatment i s presented i n Appendix H. 240 FIGURE 6.16 INDEX OF REFRACTION VERSUS ALTITUDE 241 The procedures to compute the deviations Ay g and Az g are then (refer to Figure 6.15): (a) Trace the ray launched at the i n i t i a l angle for free space until i t intersects the planet surface at point C. Tracing is through the refracting atmosphere (Appendix H). Store the y g and z coordinates of C. s (b) Compute the two-way time T taken by the beam illuminating the target at C. (c) Using the value of x from Step (b) and free space value for the speed of light, solve for x and y . This gives point B. s s (d) Compute Ay , Az which are the differences between the results • s s of Steps (a) and (c). Simulations for the model atmosphere of Venus give the following results: Near range Ay s = 0.539 km A z = s 0.026 km Mid range Ay s = 0.546 km A z = s 0.028 km Far range A y s = 0.555 km A z = s 0.030 km Since the range of variation for the correction is small (16 metres maximum), the dependence of Ayg and Az^ on T can be assumed to be negligible in the neighbourhood of a chosen target. 242 6.4.4 Use of Altimeter Data In the rectification transformation (Appendix G) , the sa t e l l i t e coordinates (x ,y z ) of the target are a function of s a t e l l i t e distance S S b H from the centre of the planet and planet radius r at target location. The accuracies of these two variables are: H = 6301.4 ± 1 km r = 6051.4 ± 0.1 km Errors in H and r would introduce errors in the computation of x g, y g and z s The VOIR altimeter measures the s a t e l l i t e altitude h from nadir, and i t s error is ±0.24 km. Because the error in h is less that the error in H and since H can be written as h + r, the s a t e l l i t e coordinates , i f expressed as a function of h and r, would be more accurately determined using altimeter data. Referring to Appendix G, z g is given by: 2 2 2 (CT /2) + H - r z = m (6.3) 's 2 H and 3z 8z * zs " H * + 77 « ( 6- A ) Substituting H = h + r, z g is given by: 243 ( c T / 2 ) 2 + H 2 + 2hr ID Z s 2 (h + r) (6.5) and 3z 9z (6.6) Table 6.2 compares the maximum error, using the nominal s a t e l l i t e o r b i t parameters, introduced i n the computation of s a t e l l i t e coordinates of both representations. It i s seen that the h representation (Equation 6.4), using the altimeter measurement data, i s better than the H representation (Equation 6.3) i n terms of r e c t i f i c a t i o n accuracy. Also, i n the h representation, ( 3 z g / 3r) Ar i s n e g l i g i b l e . 6.4.5 D i g i t a l T e r r a i n Model D i g i t a l t e r r a i n data can be incorporated i n the remapping process to correct f o r r e l i e f displacement as described i n Section 2.6.2. The only DTM a v a i l a b l e for Venus has been compiled by P e t t e n g i l l et a l [PETT79] with a 200 km gri d spacing. Note that even i f r e l i e f displacement i s not corrected, reasonable swath-to-swath r e g i s t r a t i o n can s t i l l be achieved because of the VOIR geometry, parallax of a target i n two adjacent s t r i p s being n e g l i g i b l e f o r most of the planet. In imaging from one s t r i p to the next the o r b i t i s displaced about 10 km, and i t can be shown that the elevation needs to be about 7 km above the datum to produce a 200 metre parallax. TABLE 6.2 H REPRESENTATION VERSUS h REPRESENTATION E r r o r AD = V ( A x s ) 2 + ( A y s)z + (A2 g ) 2 H R e p r e s e n t a t i o n AD Near range Mid range F a r range 1.299 ( A H - A r ) 1.272 ( A H - A r ) 1.248 ( A H - A r ) h R e p r e s e n t a t i o n Near range Mid range F a r range AD 1.299 Ah 1.272 Ah 1.248 Ah Worst AD AH = 1 km Ar = 0.1 km 1.4289 km 1.3992 km 1.3728 km Worst AD Ah = 0.24 km 0.3118 km 0.3053 km 0.2995 km 245 CHAPTER 7 CONCLUSIONS AND RECOMMENDATIONS 7.1 Conclusions The u n i f i e d approach to the r e c t i f i c a t i o n of remotely sensed d i g i t a l imagery discussed i n t h i s t hesis provides the basis f o r a system to r e c t i f y various s a t e l l i t e imagery to a common data base, y i e l d i n g a high throughput while r e t a i n i n g geometric and radiometric accuracy. (For a 5000 x 5000 Seasat image, remapping and using nearest neighbour i n t e r p o l a t i o n requires about 6 hours of computer time. On the contrary, the suggested method takes only 1/2 hour.) Hence imagery with d i f f e r e n c t sensor geometry and c h a r a c t e r i s t i c s can be overlaid as has been i l l u s t r a t e d by the Landsat/Seasat composite. The thesis has investigated the two phases of the r e c t i f i c a t i o n process — geometric transformation and i n t e r p o l a t i o n , with geometric transformation further subdivided i n t o navigation and remapping. The i n v e s t i g a t i o n was performed f o r various s a t e l l i t e s with d i f f e r e n t imaging geometries (e.g., scanner imagery and SAR imagery). Regardless of the type of imagery, high throughput i s achieved with one-dimensional i n t e r p o l a t i o n of imagery data. Based on the foregoing i n v e s t i g a t i o n s , the following conclusions can be drawn. 1. A DTM i s required f o r the r e c t i f i c a t i o n of high resolution imagery. This has been v e r i f i e d i n experiments with Seasat imagery and with synthesized TM and "SPOT imagery. T e r r a i n information i s required i n both the navigation phase and the remapping phase. 246 2. In the navigation phase, unknown parameters in the rectification transformation can be expanded as time series, the coefficients of which can be determined with the aid of GCPs and a recursive estimation technique. The number of GCPs required to attain a certain rectification accuracy depends upon the quality of the GCPs. This method has been investigated from high resolution and narrow swath width (Landsat MSS) to low resolution and wide swath width imagery (TIROS-N). It was shown that about 10 GCPs per scene were necessary for subpixel rectification accuracy. 3. For SAR imagery, another way of navigation is via automatic focusing. This method has been discussed. However, investigation of this method is left for future research. 4. To retain radiometric accuracy and high throughput using one-dimensional interpolation, three one-dimensional passes have to be applied. In each pass, a 16-point sine interpolation kernel weighted by a Kaiser window w i l l yield a maximum of 1 level RMS error out of 256 levels. This resulting accuracy is far better than cubic convolution which has been proposed by a few researchers. 5. The presence of scan gaps in TM imagery has led to the investigation of a new interpolation kernel. This kernel is able to yield an acceptable 1 level RMS error out of 256 levels in the case of overlap. The error in the case of underlap depends on the scan gap size. 247 6. The example of a Landsat/Seasat composite has succinctly demonstrated the benefits of multi-sensor i n t e g r a t i o n . The a r t i f i c i a l l y produced image composite combines the data acquired by two d i f f e r e n t sensors and s a t e l l i t e s , one a passive system emphasizing tonal q u a l i t i e s of ground cover and the other an active radar system accentuating surface texture. 7.2 Recommendations f o r Future Work As a result of the i n v e s t i g a t i o n s i n the t h e s i s , the following recommendations are made for further research work: 1. In the navigation phase, i t i s desirable to reduce the number of GCPs. It i s a very time consuming process to mark GCPs and i n many areas of the world, high q u a l i t y maps are not a v a i l a b l e . One way to reduce the number of GCPs i s to model the s a t e l l i t e a t t i t u d e by a p h y s i c a l law of s a t e l l i t e motion: rate of change of angular momentum (a function of rate of change i n atti t u d e angles) equals the t o t a l torque acting on the s a t e l l i t e . Hence i f the t o t a l torque acting on the s a t e l l i t e i s known, the v a r i a t i o n s i n the a t t i t u d e angles Can be evaluated. This method requires an extensive knowledge of the s a t e l l i t e ' s physical parameters (such as moment of i n e r t i a , mass and dimension of solar panels, control parameters, etc.) [FRIE83], 2. Again on the topic of navigation, the process of marking GCPs can be automated with image to image c o r r e l a t i o n techiques [BARN72]. Here, a small section of an image around each GCP i s stored for a p a r t i c u l a r 248 sensor. When another image with the same sensor or a d i f f e r e n t sensor i s obtained over the same area, automatic c o r r e l a t i o n can be evoked. However, the following problems must be solved: - How w i l l geometric d i s t o r t i o n s such as those present i n TM imagery a f f e c t the c o r r e l a t i o n technique? - How should sensor to sensor (with d i f f e r e n t sensor resolutions) c o r r e l a t i o n be performed? Furthermore what should be the or i e n t a t i o n of the GCP reference subimage? - How would seasonal and sun p o s i t i o n changes af f e c t the correlation? - Automatic r e g i s t r a t i o n for two' radar scenes over the same area i s d i f f i c u l t , l e t alone r e g i s t r a t i o n of radar to scanner imagery. The problems l i e i n the i n t e n s i t y migration, presence of speckles and p o l a r i z a t i o n e f f e c t s i n radar imagery. 3. Automatic focusing appears to be promising f o r the navigation of SAR imagery. This method should be investigated with Seasat SAR imagery and simulated imagery. The following should be studied: What geometric parameters can be updated? - What i s the effect of "matched f i l t e r control point" d i s t r i b u t i o n ? - Performance and l i m i t a t i o n s . 4. In the processing of SAR data, a one-dimensional remapping and i n t e r p o l a t i o n pass has been evoked to correct f o r range c e l l migration. Can t h i s be coupled with the three one-dimensional passes method to reduce 249 the t o t a l number of remapping and i n t e r p o l a t i o n passes? 5. A r a d i o m e t r i c a l l y corrected or destriped Landsat scene always shows some residual s t r i p i n g . How does t h i s residual s t r i p i n g a f f e c t the i n t e r p o l a t i o n accuracy? 6. Is i t possible to develop an i n t e r p o l a t i o n kernel the s i z e of which i s adaptive to the frequency content of the subarea i t i s operating on? Then t h i s f i l t e r must have b u i l t - i n i n t e l l i g e n c e to know or predict the gray l e v e l v a r i a t i o n s i n the subarea. For r e l a t i v e l y f l a t areas ( i n terms of i n t e n s i t y ) , a kernel s i z e of less than 16 i s adequate to r e t a i n the required radiometric accuracy of 1 l e v e l RMS out of 256; hence i n t e r p o l a t i o n time i s reduced. 7. The following i s the general problem of using t e r r a i n data to perform radiometric c o r r e c t i o n . In the correction of SAR imagery, the slope foreshortening e f f e c t causes a narrow band of high p i x e l i n t e n s i t i e s to be mapped to a wider area. While i t i s correct from the si g n a l processing point of view to use a sine function for i n t e r p o l a t i o n , a d i f f e r e n t method which also accounts for the slope effect may y i e l d better r e s u l t s for image i n t e r p r e t a t i o n purposes. An example i s to average out the i n t e n s i t y of the narrow band to the mapped wider area. 8. The r e l i e f displacement f o r Mapsat and a possible high squint SAR imagery w i l l also have an across scan component. How can the present approach be modified to accommodate the correction of such d i s t o r t i o n ? 250 9. How can the Landsat/Seasat overlay be used? An example w i l l be to show how such overlay can improve classification accuracy. Another example w i l l be to detect features which cannot be detected in only one component. 10. The resulting Landsat/Seasat composite has the high resolution of the Seasat image while retaining the spectral information of the Landsat image. This leads one to believe that i n the design of a high resolution MSS system, not a l l the spectral bands need to be of the same resolution. This is for data reduction purposes. For a set of correlated spectral bands, i t may be required only to have one band with the high resolution and a l l the other bands in the set with a relatively lower resolution. When the lower resolution bands are modulated by the corresponding high resolution band for a l l sets of spectral bands, the result w i l l be a high resolution colour image. Hence data are reduced without significantly degrading the amount of information. Experiments can be performed with Landsat MSS or TM imagery. Simulation of the low resolution bands can be done by low pass f i l e r i n g the imagery data in the corresponding bands. 11. A major area of research definitely l i e s in the production of DTMs. The following work should be investigated: The side looking mode of SPOT w i l l allow a stereo pair to be obtained and hence terrain information can be extracted of the imaged area. However the relief displacement vectors in two neighbouring passes are not parallel (except at the equator). How would this affect the correlation process to obtain 251 t e r r a i n height and the height accuracy? One technique' to obtain a contoured image i s by Moire contourography [PIR082]. In t h i s technique, an object has to be much closer to the interference grating than the i l l u m i n a t i o n to the grating. Can t h i s technique be modified to s a t e l l i t e borne imagery? Since i n t h i s a p p l i c a t i o n the t e r r a i n i s much farther away from the grating than the i l l u m i n a t i o n i s to the grating, the problem of d i f f r a c t i o n has to be overcome. A SAR system i s a suitable choice of i l l u m i n a t i o n , but w i l l the SAR data processing a n n i h i l a t e the desired contours which are formed by the p r i n c i p l e of interference? Another problem i s that the contours thus obtained do not show concavity or convexity. Another technique of obtaining a contoured image i s mentioned by Graham [GRAH74], It uses two antennas separated by a v e r t i c a l distance on the spacecraft. Hence the two images obtained w i l l create an interference pattern to produce contour l i n e s . The pos i t i o n i n g of the contour depends on the imaging geometry of the two antennas. However, the concavity and convexity of the contours are not immediately known and have to be i n f e r r e d . Also the contours are displaced according to contour elevation. '•> Vast a p p l i c a t i o n of remotely sensed imagery has already been found i n geology, a g r i c u l t u r e , f o r e s t r y , oceanography, cartography and many other earth resources management areas. With the increased res o l u t i o n and 252 number of spectral bands in the present and future generations of s a t e l l i t e sensors, the rectified imagery w i l l undoubtedly be more invaluable than ever in these application areas. 253 BIBLIOGRAPHY [ALLA78] [AVER8.1] [BAKE75] [BARN72] [BENN78] [BENN79] [BENS80] [BERN76] [CAR075] [CNES81] [COLV79] [CUMM79] Al l a n , M. M. "DTM A p p l i c a t i o n i n Topographic Mapping", Photogrammetric Engineering and Remote Sensing, Vol.44, No.12, Dec.1978, pp.1513-1520. Avery, J . E. and J . S. Hsieh, "Geometric Correction Resampling for the Landsat-D Thematic Mapper", International Geoscience and Remote Sensing Symposium, Washington, D.C., June 1981. Baker, J . R. and E. M. M i k h a i l , "Geometric Analysis and R e s t i t u t i o n of D i g i t a l M u l t i s p e c t r a l Scanner Data Arrays", LARS Information Notes 052875, Purdue Un i v e r s i t y , 1975. Barnea, D.T. and H. Silverman, "A cl a s s of Algorithms for Fast D i g i t a l Image Regi s t r a t i o n " , IEEE Transactions on Computers, Vol.C-21, Feb.1972, pp.179-186. Bennett, J.R. and I.G. Cumming, " D i g i t a l Technique f o r the Multilook Processing of SAR Data With A p p l i c a t i o n to Seasat-A", F i f t h Canadian Symposium on Remote Sensing, V i c t o r i a , B.C., Aug.1978, pp.506-516. Bennett, J.R. and I.G. Cumming, "A D i g i t a l Processor for the Production of Seasat Synthetic Aperture Imagery", Surge Workshop, Esrim, F r a s c a t i , I t a l y , July 1979. Benson, M., " D i g i t a l Processing of Seasat-A SAR Data using L i n e a r Approximations to Range C e l l Migration Curves", IEEE International Radar Conference, 1980, pp.176-181. Bernstein, R., " D i g i t a l Observation Sensor Data", Image Processing of Earth IBM Journal of Research arid Development. Jan.1976, pp.40-57. Caron, R.H. and K.W. Simon, "Attitude Time Series Estimation f or R e c t i f i c a t i o n of Spaceborne Imagery", Journal of Spacecraft and Rockets, Vol.12, No.l, Jan.1975, pp.146-151. "SPOT to Ground Stat i o n France, May, 1981. Interface Document' CNES , Colvocoresses, A.P., "Proposed Parameters f o r Mapsat", Photogrammetric Engineering and Remote Sensing , Vol.45, No.4, Apr.1979, pp.501-505.. Cumming, I.G. and J.R. Bennett, " D i g i t a l Processing of Seasat SAR Data", IEEE International Conference on Acoustics, Speech and Si g n a l Processing , Washington, D.C., Apr.1979, pp.710-718. 254 [DEB079] [DOYL78] [DUFL67] [EBNE76] DeBoor, C , "A Practical Guide Springer-Verlag, 1979, Chapter IV. to Splines", [EKLU72] [FRIE81] [FRIE83] [ G 0 L D 8 1 ] [GR/H74] [GSFC82] Doyle, F.J., "Digital Terrain Models, An Overview", Photogrammetric Engineering and Remote Sensing, Vol.44, No.12, Dec.1978, pp.1481-1485. Du Plessis, R.M., "Poor Man's Explanation of Kalman Fi l t e r i n g " , Northern American Rockwell Corporation, June 1967. Ebner, H., "A Mathematical Model for Digital Rectification of Remote Sensing Data", XIII Congress of the International Society for Photogrammetry, Commission III, Helsinki, Finland, July 1976. Eklundh, J.O., "A Fast Computer Method for Matrix Transposing", IEEE Transactions on Computers, Vol.C-21, June 1972, pp.801-803. Friedmann, D.E., '*Two-Dimensional Resampling of Line Scan Imagery by One-Dimensional Processing", Photogrammetrl c Engineering and Remote Sensing. Vol.47, No.10, Oct.1981, pp.1459-1467. Friedmann, D.E., J.P. Friedel, K.L. Magnussen, R. Kwok and S. Richardson, "Multiple Scene Precision Rectification of Space Borne Imagery with Very Few Ground Control Points", Photogrammetric Engineerine and Remote Sensing. Vol.49, No.12, Dec.83, pp.1657-1667. Goldbogan, G.C., "PRIM: A Fast Matrix Transpose Method", IEEE Transactions on Software Engineering, Vol.SE-7, No.2, Mar.1981, pp.255-257. Graham, L.C., "Synthetic Interferometer Radar for Topographic Mapping", Proceedings of IEEE, Vol.62, No.6, June 1974, pp.763-768. Also a selected paper in [KOVA76], pp.275-280. "Landsat-D to Ground Station Interface Description", Revision 5, Goddard Space Flight Centre, Maryland, Aug.1982. [HARD71] Hardy, R.L., "Multiquadratic Equation of Topography and Other Irregular Surfaces", Journal of Geographical Research". Vol.76, No.8, 1971. [HEIL78] H e l l , R.V. and S.M. Brych, "An Approach for Consistent Topographic Representation of Varying Terrain", DTM Symposiurn, St. Louis, Missouri, May 1978, pp.397-410. 255 (HORN78 ]. Horn, B.K.P. and B.K. Bachman, "Using Synthetic Images with Surface Models", Communications of the ACM. Vol.21, No.11, Nov.1978, pp.914-924. [HORN79] Horn, B.K.P. and B.J. Woodham, "Landsat MSS Coordinate Transformations", Proceedings of the Fifth International Symposium on Machine Processing of Remotely Sensed Data. 1979, pp.59-68. [H0RN81] Horn, B.K.P., ' H i l l Shading and the Reflectance Map", Proceedings of the IEEE, Vol.69, No.l, Jan.1981, pp.14-47. [HOVA79] Hovanessian, S.A., "Introduction to Synthetic Array and Imaging Radars", Hughes Aircraft, Mar.1979. [JAME82] James, W., "Unveiling Venus with VOIR", Sky and Telescope, Feb.1982, pp.141-144. [JERR77] J e r r i , A.J., 'The Shannon Sampling Theorem - Its Various Extensions and Applications: A Tutorial Review", Proceedings of the IEEE, Vol.65, No.11, Nov.1977, pp.1565^1596. [JPL77] "Seasat-A Instrument Data Processing System Detail Functional Specification", Vol.1, JPL, Pasadena, California, June 1977. [KELL77] Kelly, R.E., P.R.H. McConnell and S.J. Mi ldenberger, 'The Gestalt Photomapper System", Photogrammetric Engineering and Remote Sensing, Vol.43, No.11, Nov.1977, pp.1407-1417. [KONE76] Konecny, G., "Mathematical Models and Procedures for the Geometric Restitution of Remote Sensing Imagery", XIII Congress of the International Society for Photogrammetrv. Commission III, Helsinki, Finland, July 1976. [K0RN61] Korn, G. and T. Korn, "Mathematical Handbook for Scientists and Engineers", McGraw-Hill, 1961, Sections 21.4 and 21.8. [KOVA76] Kovaly, J.J., "Synthetic Aperture Radar", Artech, 1976. [KRAT72] Kratky, V., "Photogrammetric Solution for Precision Processing of ERTS Images", XII Congress of the International Society of Photogrammetry, Ottawa, June 1972. [KRAU78] Kraus, K., "Rectification of Multispectral Scanner Imagery", Photogrammetric Engineering and Remote Sensing. Vol.44, No.4, Apr.1978, pp.453-457. [LARS80] Larson, J., "Rectification of Digital Images for Remote Sensing Analysis", Department of Photogrammetry, Royal Institute of Technology, Stockholm, Sweden, 1980. 256 [LATH65] Lathi, B.P., "Signals, Systems and Communications", John Wiley and Sons, 1965, Chapter 11. [LECK80] Leckie, D.G., "Use of Polynomial Transformation for Registration of Airborne Digital Line Scan Images", Sixth Canadian Symposium on Remote Sensing, Halifax, N.S., May 1980, pp.635-641. [LEON70] Leondes, C.T., "Theory and Applications of Kalman Fil t e r i n g " , University of California, Los Angeles, Feb.1970. [LINF39] Linfoot, E.H. and W.M. Shepherd, "On a Set of Linear Equations (II) " . Journal of Mathematics, Vol.10, 1939, pp.84-98. [LITT80] L i t t l e , J.J., "Automatic Rectification of Landsat Images Using Features Derived from Digital Terrain Models", M.Sc. Thesis, U.B.C., Vancouver, 1980. [MALI73] Maling, D.H., "Coordinate Systems and Map Projections", George Philip and Son Limited, London, England, 1973, p.194. [MCCA80] McCarthy, J.F., "Summary of Atmospheric Refraction Study", Hughes Aircraft internal memo, June 1980. [MDA76] "Feasibility Study Signal Processing System for Satellite-Borne Synthetic Aperture Radar", MDA Technical Report for CCRS and CRC, Ottawa, Feb.1976. [MDA77a] "Considerations in the Implementation of the Range/Doppler Domain Approach to Seasat-A SAR Data Processing", MDA, Dec. 1977. [MDA77b] "A Theoretical Foundation for the Design of a Ground Data Processor for a Sa t e l l i t e Borne Synthetic Aperture Radar", MDA, Aug.1977. [MDA78] '*The Development of an Ice Status System", MDA, Mar. 1978. [MDA79] "SAR Seasat Processor Functional Description", MDA, Apr.1979. [MILG75] Milgram, D., "Computer Methods for Creating Photomosaics", IEEE Transactions on Computers, Vol.C-24, No.ll, Nov.1975' pp.1113-1119. [MILG77] Milgram, D., "Adaptive Techniques for Photomosaics", IEEE Transactions on Computers, Vol.C-26, No.ll, Nov.1977, pp.1175-1180. [MILL71] Millar, R.A., "The Principles of Kalman Fil t e r i n g " , CRC, Ottawa, Sept.1971. 257 [NASA76a] "Real-Time Data Systems for the TIROS-N Spacecraft Series", NASA Document, Mar.1976. [NASA76b] "ERTS Data Users' Handbook", NASA Document No.71SD4249, Sept.1976. [ORTH77] Orth, R. and D.S. Sloan, "A Self-Contained Landsat Data Reception and Precision Cartographic Image Production System", Proceedings of the International Symposium on Image Processing, Interactions with Photogrammetry and Remote Sensing, Graz, Austria, Oct.1977, pp.189-196. [ORTH78a] Orth, R., K. Brydges and F. Wong, "Photomaps from Precision Rectified Landsat Imagery", Fifth Canadian Symposium on Remote Sensing, Victoria, B.C., Aug.1978, pp.292-295. [ORTH78b] Orth, R., F. Wong and J.S. MacDonald, 'The Production of 1:250,000 Maps of Precision Rectified and Registered Landsat Imagery using the MDA Image Analysis System: I n i t i a l Results", Twelfth Symposium on Remote Sensing of Environment. Manila, Phillipines, Apr.1978, pp.2163-2176. [OTEP78] Otepka, G., "Practical Experience in the Rectification of MSS Images", Photogrammetric Engineering and Remote Sensing. Vol.44, No.4, Apr.1978, pp.459-467. [PETT79] Pettengill, G.J., E. Eliason, P.G. Ford, G.B. Loriot, H. Masursky and G.E. McGill, "Pioneer Venus Radar Results: Altimetry and Surface Properties", Journal of Geophysical Research, Dec.1979. [PEUC78] Peucker, T.K. , R.J. Fowler, J.J. L i t t l e and D.M. Mark, "The Triangulated Irregular Network", DTM Symposium, St. Louis, Missouri, May 1978, pp.515-540. [PIR082] Pirodda, L., "Shadow and Projection Moire Techniques in Absolute and Relative Mapping of Surface Shapes", Optical Engineering, Vol.21, No.4, Aug.1982, pp.640-649. [PRAK81] Prakash, A., and E.P. Beyer, "Landsat-D Thematic Mapper Image Resampling for Scan Geometry Correction", Proceedings of the Seventh International Symposium on Machine Processing of Remotely Sensed Data. June 1981, pp.189-200. [RABI75] Rabiner, L.R. and B. Gold, 'Theory and Application of Digital Signal Processing", Prentice-Hall, 1975, Section 3.11. [RAMA77] Ramapriyan, H.K., "Data Handling .for the Geometric Correction of Large Images", IEEE Transactions on Computers, Vol.C-26, No.ll, Nov.1977, pp.1163-1167. 258 [ RICH72] Richar.dus, P., and R.K. Adler, "Map Projections for Geodesists, Cartographers and Geographers", North-ttolland 1972, Section 7.5. [SAKR68] Sakrison, D.J., "Communication Theory: Transmission of Waveforms and Digital Information", Wiley, 1968, Section 2.7. [SCHU76a] Schumaker, L., "Fitting Surfaces to Scattered Data", Approximation Theory II, ed. G. Lorentz, Academic Press, 1976. [SCHU76b] Schut, G.H., "Review,of Interpolation Methods for Digital Terrain Models", XIII Congress of the International Society for Photogrammetry, Commission III, Helsinki, Finland, July 1976. [SHLI79] Shlien, S., "Geometric Correction, Registration and Resampling of Landsat Imagery", Canadian Journal of Remote Sensing, Vol.5, No.l, May 1979, pp.74-89. [SIM075a] Simon, K.W., S.S. Rifman and R.H. Caron, "Precision Digital Processing of Landsat MSS Data", Eleventh International Symposium on Space Technology and Science. Tokyo, Japan, July 1975. [SIM075b] Simon, K.W., "Digital Image Reconstruction and Resampling for Geometric Manipulation", Proceedings IEEE Symposium on Machine Processing of Remotely Sensed Data, June 1975, pp. 3A-1 - 3A-11. [UCLA79] [WOOD80] [YEN56] [ZOBR82] "UCLA Kalman Filtering Lecture Notes", UCLA, California, Mar.1979. Woodham, R.J., "Using Digital Terrain Data to Model Image Formation i n Remote Sensing", SPIE. Image Processing for Missile Guidance. Vol.238, July 1980, pp.361-369. Yen, J.L., "On Non-Uniform Sampling of Bandlimited Signals", IRE Transactions on Circuit Theory. Dec.1956, pp.251-257. Zobrist, A., "Computational Aspects Imagery", Proceedings of the of Remapping Digital NASA Workshop on Registration and Rectification, NASA, July pp.358-366. 1982, 259 APPENDIX A LANDSAT-4 SENSORS A.l General Description The Landsat-D s a t e l l i t e has been successfully launched in late 1982 to become Landsat-4. The sensors on board consist of a multispectral scanner (MSS) and a thematic mapper (TM). The s a t e l l i t e and the sensors have the following characteristics: Orbit inclination: 99° Nodal period: 99 minutes S a t e l l i t e altitude: varies from 685 km to 740 km, 705 km at 40 N°latitude Swath width: 185 km Pixel resolution: 58 m x 79 m MSS 30 m x 30 m TM Bands 1,2,3,4,5 and 7 120 m x 120 m TM Band 6 The MSS sensor i s similar, both in spectral bands and resolution, to the MSS sensors in the previous Landsat s a t e l l i t e s . Only the TM high resolution bands are considered in this thesis due to the complex imaging geometry of the detectors and the high susceptibility to geometric distortions as a result of the high resolution. The TM is a bidirectional scanner, imaging the earth in the forward and reverse scanning directions across the ground track as shown in Figure A . l . The scanning i s accomplished by an oscillating mirror. A scan line corrector (SLC) helps to align the ground track projections together. Without the SLC, 260 DETECTORS NCI " • « O 0 « 0 ° •AND 7 DETECTOR ORIENTATION AT PRIMARY FOCAL F-LANE DETECTOR ORIENTATION AT COOLED fOCAL PLANE SPACECRAFT COORDINATES DETECTORS tiZ. 1 SCAN DIRECTION LANDSAT AROUND TRACK SOUTH FIGURE A . l DETECTOR ARRAY PROJECTIONS ON GROUND TRACK [GSFC82] 261 t h e r e e x i s t s a s i g n i f i c a n t amount of o v e r l a p and u n d e r l a p f r o m swath t o swa t h . The f u n c t i o n s of t h e SLC a r e shown i n F i g u r e A.2 w h i c h shows t h e ground p r o j e c t i o n o f t h e s c a n w i t h t h e s a t e l l i t e s t a t i o n a r y . W i t h t h e v e l o c i t y v e c t o r o f t h e s a t e l l i t e , t h e f i n a l and r e q u i r e d p r o j e c t i o n i s shown i n F i g u r e A.2c. The c h a r a c t e r i s t i c o f TM imagery i s t h e p r e s e n c e of s c a n gaps t o be d e s c r i b e d i n what f o l l o w s . The a c t u a l imagery o b t a i n e d w i l l not be as i d e a l and i s shown i n F i g u r e A.3, because t h e SLC i s s e t f o r n o m i n a l s a t e l l i t e ground speed. But t h e ground speed d e v i a t e s f r o m n o m i n a l because o f : s a t e l l i t e a l t i t u d e v a r i a t i o n , and - s a t e l l i t e v e l o c i t y v a r i a t i o n . F u r t h e r m o r e , t h e TM s e n s o r s s u f f e r f r o m h i g h f r e q u e n c y v i b r a t i o n j i t t e r and f r o m t h e b o w t i e e f f e c t ( p i x e l p r o j e c t i o n on t h e ground i s l a r g e r f o r a l a r g e r s c a n a n g l e ) • These d i s t o r t i o n s i n t r o d u c e u n e v e n l y spaced samples when p r o j e c t e d on t h e ground. T h i s i s shown i n F i g u r e A.3c i n w h i c h A^ d e n o t e s t h e i n t e r l i n e s p a c i n g w i t h i n a swath and A j d e n o t e s t h e s p a c i n g between t h e two end d e t e c t o r s i n two swaths. I d e a l l y , AQ s h o u l d be t h e same as A^. But due t o t h e v a r i o u s d i s t o r t i o n s , AQ may not be e q u a l t o A^. N o t e t h a t A^ i s c o n s t a n t o v e r t h e t i m e i n t e r v a l of i m a g i n g a s c e n e . As a r e s u l t , t h e r e a r e 15 e q u a l sample s p a c i n g s f o l l o w e d by one u n e q u a l sample s p a c i n g , s i n c e each swath c o n t a i n s 16 l i n e s o f i m a g e r y . a) Ground T r a c k Without SLC (t>) SLC M o t i o n 1 1 1 1 1 1 M 1 I I 1 (c) I d e a l Ground Track FIGURE A.2 TM GROUND TRACK 263 A B l i n e s ] nth scan }l n * l ) t h _ scan Ideal Scans 1 :) nth seen ( n + l ) t h scan (b) D i s t o r t e d Scans A (G>0) B (G=0) C (G<0) G >0 fl0 A-0 ( c ) S e c t i o n s FIGURE A.3 TM SCAN GAPS 264 A.2 Scan Gap Analysis A characteristic of TM imagery is the presence of scan gaps. This section analyzes the scan gap magnitude. Figure A.3 shows that the ground coverage contains a certain amount of underlap (missing data) and overlap (redundant data). The underlap/overlap is sometimes referred to as scan gap and is defined as: G = Aj - A Q (A.l) Therefore, when G i s 0, a l l the samples are evenly spaced: when G is positive, underlap occurs and when G i s negative, overlap occurs. For TM imagery, G l i e s between +2 lines. The distortions that contribute to the scan gap are discussed in what follows. 1. J i t t e r J i t t e r i s caused by vibration in the mirror assembly. Anticipated RMS j i t t e r error, referenced to the spacecraft coordinate system, is shown in Table A.l [GSFC82]. Significant error i s not expected to occur above 77 Hertz. The non-zero scan gap is partly caused by j i t t e r of frequencies higher than 0.4 Hertz. For example, j i t t e r of 7 Hertz or more corresponds to one cycle in less than 32 lines, which is two swaths of data. The across scan distortion, caused by j i t t e r (3a level) i s thus between ±6 metres. Hence, the maximum scan gap is 12 metres at the 3o level in the worst case i f the j i t t e r s for two consecutive swaths are exactly out of TABLE A . l ERROR MAGNITUDE CAUSED BY JITTER A c r o s s - S c a n Frequency Number o f Number o f Swaths E r r o r Magnitude E r r o r Magnitude Range ( H e r t z ) L i n e s p e r C y c l e p e r C y c l e (3a) i n a r c - s e c (3a) i n metres 0.4 t o 7 560 t o 32 35 t o 2 0.9 a l l axes 3.1 G r e a t e r than Less than 32 Less than 2 2.79 r o l l 9.5 7 0.60 p i t c h 2.1 0.90 yaw 0.3 266 phase. J i t t e r also introduces d i s t o r t i o n i n the along scan d i r e c t i o n , but t h i s i s slow varying as shown below. Taking the highest s i g n i f i c a n t frequency of j i t t e r to be 77 Hertz and the p i x e l sampling rate 9.6 x 10 ^ sec, one cycle corresponds to at least -6 1/(77 x 9.6 x 10 ) = 1352 p i x e l s 2. S a t e l l i t e A l t i t u d e V a r i a t i o n The a l t i t u d e of the Landsat-4 o r b i t , considering both o r b i t e c c e n t r i c i t y and the earth shape, w i l l vary between 685 and 740 km. At 705.3 km, the swath width i n the along scan d i r e c t i o n i s 480 metres. V a r i a t i o n i n the o r b i t a l t i t u d e w i l l produce a scale change i n the l i n e spacing within a swath but the distance between swath centres w i l l not change (assuming constant s a t e l l i t e v e l o c i t y ) . The contribution to the scan gap can be calculated by the r e l a t i o n : G = 480 (1 - TQI^) (A.2) where h i s the o r b i t a l t i t u d e i n km. For the s p e c i f i e d a l t i t u d e range, G varies from 13.8 to -23.6 metres. 3. S a t e l l i t e V e l o c i t y V a r i a t i o n V e l o c i t y v a r i a t i o n causes a change i n the spacing between swath centres but not within the swath and hence, a f f e c t s the scan gap. 267 Assume t h a t t h e major c o n t r i b u t i o n t o v e l o c i t y v a r i a t i o n i s a l t i t u d e v a r i a t i o n , t h e n t h e s a t e l l i t e v e l o c i t y f o l l o w s f r o m K e p l e r ' s Law: 0.5 (A.3) where K i s a p r o p o r t i o n a l i t y c o n s t a n t , h i s o r b i t a l t i t u d e , R i s t h e l o c a l e a r t h r a d i u s a t t h e n a d i r p o i n t , and R + h i s t h e l o c a l s a t e l l i t e r a d i u s . The g round speed i s t h e n 1.5 v G = KR ( z~TT ] (A.4) and t h e swath c e n t r e s p a c i n g i s W = KRT { — i - r \ (A.5) ( R i h ) where T i s t h e swath t i m e . S i n c e W i s about 480 met re s w i t h an o r b i t a l t i t u d e o f 705.3 km, t h e v a l u e o f W a t any o t h e r o r b i t a l t i t u d e h i s t h e n g i v e n by: R / R„ + 705.3 X 1* 5 - r0[ — W = 480.0 - " u (A.6) where RQ i s t h e l o c a l e a r t h r a d i u s a t t h e s u b s a t e l l i t e p o i n t when t h e o r b i t a l t i t u d e i s 705.3 km. The maximum d i f f e r e n c e between R and R^ i s 22 km on t h e e a r t h s p h e r o i d , and t h i s d i f f e r e n c e g i v e s an e r r o r o f l e s s t h a n one met re i n W i f R i s assumed t o be e q u a l t o R^ i n t h e above e q u a t i o n . 268 When the s a t e l l i t e o r b i t i s highest at 740 kmj, the value of W i s 476.5 metres and when i t i s lowest at 685 km, the value of W i s 482.1 metres. Hence, the contribution to the scan gap varies from -3.5 metres to 2.1 metres. 4. Scan Skew This i s a combined effect of the SLC rate, s a t e l l i t e a l t i t u d e , ground v e l o c i t y and earth r o t a t i o n . To produce the i d e a l ground track as shown i n Figure A.2c, the s a t e l l i t e a l t i t u d e and ground v e l o c i t y roust exactly compensate f or the SLC. This occurs only at 40°N i n l a t i t u d e where the a l t i t u d e i s 712.8 km and ground v e l o c i t y i s 6.821 km/sec [PRAK81]. At the other a l t i t u d e s , the ground v e l o c i t y changes according to Equation A . 4 . The scan skew components are shown i n Figure A.4 i n which A = ground v e l o c i t y x active scan time, B = SLC rewind ground distance, C = distance t r a v e l l e d during rewind, E = earth r o t a t i o n i n along-track d i r e c t i o n i n one swath time, and dp dj = ground distance between swath centres at both ends. From the fi g u r e d 0 = B + C and d 1 = 2(A + E) - B + C At the extreme a l t i t u d e s , the values are A l t i t u d e km d^ metres d^ metres (d^ - d^) metres 685 473.3 516.7 -42.9 740 504.6 474.0 30.6 269 Maximum scan angle magnitude Zero scan angle Coverage of f i r s t detector Coverage of last detector Maximum scan angle magnitude FIGURE A.5 BOWTIE EFFECT 270 If scan skew d is defined as the skew measured from the swath centre, then d = (d 0 - dj) / 2 and i t varies between ±21.5 metres at 685 km altitude and between ±15.3 metres at 7AO km. 5. Bowtie Effect The bowtie effect is caused by the earth curvature and the way the detectors are positioned in the focal plane. This effect is illustrated in Figure A.5 which shows the positions of the f i r s t and last detectors of a particular band. The swath width in the along track direction is narrower in the centre than at both ends, hence causing a line displacement. The line displacement d, as shown in the figure, is 2.3 metres in the worst case. Hence, this gives a scan gap G of 0 to -A.6 metres, depending on the pixel number. Including the effect of altitude change, the range of gap G i s from 0 to -A.8 metres. This distortion i s of high frequency which is obvious from the fact that the line displacements are in opposite directions for the f i r s t . and last detectors. 6. Earth Rotation Earth rotation causes the line spacing between two swaths to widen. The magnitude of scan gap G at the swath centre (nadir) is 271 G = T u cos 8. (A.7) e 1 where T Is the swath time, to is the earth rotation velocity at the e equator and 6^ i s the s a t e l l i t e inclination. For the northbound and southbound passes, G is always positive. 7. Swath Time Variation The time taken to sweep one swath of imagery plus repositioning for the sensor to start imaging the next swath of data may vary. The magnitude of variation is not precisely known yet. The effect of this variation is that swaths of data are distorted in the across scan direction. Such distortion is of high frequency since the variation is different between swaths. 8. Summary Table A.2 shows that the scan gap varies between -2 A^ and 2 AQ where A^ is the line spacing within a swath. In the table, the worst case for each scan gap contribution factor is assumed so that the factors a l l add constructively whereverpossible, The total scan gap size in the table does not include swath time variation, the magnitude of which is not known at this moment. 272 TABLE A.2 WORST CASES OF SCAN GAP SIZE A l t i t u d e 685 km 740 km L i n e s p a c i n g w i t h i n a swath A 0 29.14 m 31.48 m SCAN GAP J i t t e r (3a) 12.0 m -12.0 m A l t i t u d e V a r i a t i o n 13.8 m -23.6 m V e l o c i t y V a r i a t i o n 5.6 m -7.3 m Scan Skew 21.5 m -15.3 m Bowtie E f f e c t 0 m (swath c e n t r e ) -4.8 m (swath edge) E a r t h R o t a t i o n 5.2 m 5.2 m Swath Time V a r i a t i o n not known not known T o t a l ( n o t i n c l u d i n g swath time v a r i a t i o n ) 58.1 m = 2.0 A 0 -57.8 m = -1.8 A 0 273 APPENDIX B SPOT SENSORS B.l General Description The SPOT s a t e l l i t e is designed to have two modes of operations : multispectral and panchromatic. The multlspectral mode employs a multilinear array (MLA) of detectors. The panchromatic mode also employs a linear array of detectors and sometimes i s referred to as panchromatic linear array (PLA). A key feature of the sensors i s the provision for off nadir viewing, with the off nadir angle adjustable to +26°. This feature provides stereo viewing capability. The s a t e l l i t e and sensors have Orbit inclination : Nodal period : Sa t e l l i t e altitude : Swath width : MLA pixel resolution : PLA pixel resolution : The sensors have no mechanical using the "push-broom" scanning modi a linear array of 6000 detectors i n of the other three bands; and line the following characteristics: 98.7° 101.5 minutes 820 km nominal 60 km nadir pointing 80 km off nadir angle = 26° 20 m nadir pointing 26 m off nadir angle = 26° 3 bands 10 m nadir pointing 13 m off nadir angle = 26° moving parts; images are obtained !. Each line of the image is formed by the panchromatic band, 300 in each scanning is performed electronically. 274 Successive l i n e s of the image are produced as a r e s u l t of the s a t e l l i t e motion along i t s o r b i t . B.2 Scan Angle The detectors are spaced, except f o r between groups of detectors, equally apart. The scan angle can be expressed mathematically with the aid of Figure B.l as follows: Let £Q = angle of centre p i x e l made with the v e r t i c a l , 6 = angle of imaging p i x e l made with the v e r t i c a l , p = imaging p i x e l number, the f i r s t p i x e l being 1, 6 = angle subtended by sensor from the f i r s t p i x e l to the l a s t p i x e l a f t e r removing the detector gaps, r = spacing between two detectors i n f o c a l plane, f = f o c a l length, and N = t o t a l number of detectors. Then, the angle g i s given by: 8 = e Q + tan" 1 { [tan(0/2)] - ( p - l ) r / f } (B.l) B.3 Resolution The r e s o l u t i o n i n the along scan d i r e c t i o n i s then a function of o f f nadir pointing angle and i s given by: d 0 Resolution = R -r-r AB a p where 0 and R are defined i n Figure B . l . The res o l u t i o n thus obtained i s plotted i n Figure B.2. 275 FIGURE B . l IMAGING PIXEL AND ITS SCAN ANGLE RESOLUTION IN METRES -t-.5' OFF-AXIS POINT ANGLE 5'* 10° 1 ° 20° 25* FIGURE B.2 RESOLUTION VERSUS OFF-AXIS POINTING ANGLE IN SPOT IMAGERY 277 APPENDIX C RECURSIVE ESTIMATION FILTER C.l Background In 1960, R.E. Kalman introduced a set of matrix recursion relations for optimally estimating noise corrupted measurement data. This set of equations i s known as the discrete time Kalman f i l t e r [LEON70]. The parameters to be estimated are collected in a state vector. If these parameters are time independent, the Kalman f i l t e r degenerates to a recursive estimator. A recursive estimation procedure yields an estimate following each measurement. Specifically, i t enables one to estimate an unknown quantity based on the last estimate and the present measurement. In general, the procedure can be formulated as: x(i) = k y y(i) + k x x(i - 1) (C.l) where x is the estimate, y the measurement, i the incremental index number, and k and k are predetermined coefficient matrices. Such a procedure x y eliminates the need for storing, and performing calculations on, any amount of past data. Each new measurement is utilized as i t is received. This appendix f i r s t describes briefly the underlying principles of recursive f i l t e r i n g using the recursive estimation procedure. Then i t shows how the recursive estimator can be applied to the state vector made up of the attitude time series coefficients using GCPs as measurement data. The method can be generalized to update SAR orbit trajectory data 278 using GCPs or possibly the FM rate as measurement data. C.2 Recursive Estimator Equations F i r s t the one-dimensional case i s considered. Let x^ and be two 2 2 independent estimates of x with variances and o^ resp e c t i v e l y . Then the best estimate, or minimum variance estimate, can be written as: x = (1 - w) Xj + w = Xj + w(x£ - Xj) (C.2) where the weight w can be obtained by minimizing the variance of x. The variance of x i s given by: a 2 = E[(x - E ( x ) 2 ] 2 where E(x) denotes the expected value of x. By se t t i n g da /9w = 0, the weight w i s given by: 2 °1 w= (C.3) °1 + C T 2 2 ^ 2 2 * °2 X l °1 X2 a l Hence x = ^ \~ = x ] + 2 2 ^ X2 ~ X l ^ ( c«4) a i + °2 °1 + C T 2 2 2 2 01 a2 2 and a z = — j — 1— = a * (1 .- w) (C.5) °1 + a2 These equations can e a s i l y be extended to a multi-dimensional case Furthermore, the two estimates are not of equal dimensions. Let x and _y denote two variables related by: 279 v = M x (C6) Let x, and _y_„ be the corresponding estimates with covariance P. and P respectively. Then the optimal estimate of x i s given by: x = x + W'M (y - x (C7) By minimizing the covariance of x, i t can be shown that W ' = Pj M [M P M + P 2] -1 (C8) and the covariance of x_ i s given by: P = (I - W'M)P 1 (C9) Equations C.7, C.8, and C.9 are the key recursive estimator equations. C.3 Application to Attitude Time Series Coefficients Recursive estimation is well suited to the rectification problem with the ECR coordinates of GCPs as measurement data, and the attitude time series and s a t e l l i t e altitude coefficients as the state vector. The nominal values of pitch, r o l l , yaw and s a t e l l i t e altitude are f i r s t used i n the transformation for the f i r s t GCP selected. The ECR-coordinates of the GCP are computed and compared with the measured (from mapsheets or surveys) values. The difference is then fed to the recursive estimater which refines the coefficients. Figure C.l shows how the equations are used in the process. After each GCP is selected, the state vector and i t s covariance matrix are updated. Here x denotes the state vector, P i t s covariance matrix, y 280 I n i t i a l i z e X(0) and P(0) b a s e d on a p r i o r i knowledge o f system s t a t i s t i c s i=0 i = i + l Compute optimum w e i g h t i n g c o e f f i c i e n t m a t r i x W ' = P ( i ) M T [ M P ( i ) M T+ E J " 1 13x3 13x3 3x13 13x3 3x3 13x3 13x13 Make optimum e s t i m a t e x ( i ) = X l i - l ) + b [ y ( i ) - M X.(i-1) ] 13x1 13x1 13x3 3x13 3x13 A d j u s t P t o a c c o u n t f o r optimum e s t i m a t e P ( i ) = ( I - W M) P ( i ) 13x1 13X1 1 3 x 3 3x13 B x l NOTE: THE NUMBER BELOW EACH MATRIX DENOTES ITS DIMENSION. THERE ARE FOUR COEFFICIENTS (CUBIC POLYNOMIAL) FOR EACH OF THE THREE ATTITUDE ANGLES AND ONE FOR SATELLITE ALTITUDE FIGURE C . l LEAST SQUARES RECURSIVE FILTER EQUATIONS FOR ATTITUDE TIME SERIES MODEL 281 the measured ECR coordinates of the current GCP, E the measurement error covariance matrix (the level of confidence in identifying the GCP in the image) and I the identity matrix. 282 APPENDIX D BILINEAR AND AFFINE TRANSFORMATIONS This appendix derives the coefficients of the bilinear and affine transformations. D.l Bilinear Transformation Let (x,y) be the output grid and (u,v) be the input grid coordinates; Let the output grid be segmented into rectangular blocks. The bilinear transformation for each block can be represented by: u = a^ + a 1 Ax + a 2 Ay + a^ Ax Ay (D.la) v = b Q + b1 Ax + b 2 Ay + b 3 Ax Ay (D.lb) where Ax and Ay are measured from the block's top left corner as shown in Figure D.l. If the mapping of the 4 corners of the block is known, then the coefficients a^ and b^ ( i = 0 to 3) can be determined as follows: a0 = U l aj= (u 2 - up /X a 2 = (u^ - u ^ / Y a 3 = (u L + u 3 - u 2 - u A) / (XY) b0 = V l V ( V2 " V 1 X b 2 = (v 4 - Vj) / Y b 3 = (v x + v 3 - v 2 - v 4) / (XY) where X and Y are the block dimensions and (u , v i)» 1 = 1 to 4, are 283 MAPPING (X., Y.) < • (U , V ) , i = 1, 2, 3, 4 i i i i FIGURE D . l REMAPPING BLOCK 284 coordinates of the mapped corner points in the input grid as shown in Figure D.l. The mapping can easily be proved to be continuous from one block to the surrounding blocks. This can also be deduced from the fact that the same two corner points, common to two neighbouring blocks, are used to evaluate their respective bilinear coefficients. D.2 Affine Transformation If a^ and b^ are set equal to 0 in Equation D.l, the resulting transformation i s affine. Then the other six coefficients are given by a least squares f i t : and T - I T (M M) M (D.3a) where T - I T (M M) M (D.3b) M = 0 0 1 • X 0 1 0 Y 1 X Y 1 (D.4) 285 Upon simplification, a0 = (" U1 + U2 " U3 + U 4 ) / ( 2 X ) a1 = (-Uj - u 2 + u 3 + u A) / (2Y) a_ = u + (-u + u„ + u„ - u.) / 4 Z 1 1 Z J * (D.5) b0 = (~ V1 + V2 " V3 + V / ( 2 X ) bj = (-Vj -v 2 + v 3 + v 4) / (2Y) b 3 = vj + (-vj + v 2 + v 3 - v 4) / 4 This transformation has the restriction that parallel lines are. always mapped into parallel lines, but the mapping i s not continuous across blocks. 286 APPENDIX E CLOSED FORM TRANSFORMATION FOR SCANNER IMAGERY This appendix i s to derive a transformation between the sensor instantaneous viewing direction and the geographic coordinates of any image point. The derivation starts by postulating a unit scan vector viewing an image in the platform coordinates. Then i t is transformed into the sa t e l l i t e coordinates and fi n a l l y into the ECR coordinates. A l l vectors are in ECR coordinates unless otherwise stated. Relative orientations of various coordinate systems are shown in Figure E . l . E.l Basic Formulation Let u be the unit vector pointing from the sa t e l l i t e to a point on the earth, R the ground point vector and R the sat e l l i t e position vector -g -s as shown in Figure 4.6. Then R = R + Du (E.l) -g -s where D i s the distance between the sa t e l l i t e and the ground point. The f i r s t step is to solve for D. Let x be the geocentric latitude of the ground point, a the earth's semi-major axis and e i t s eccentricity. Define e' by: 2 e' = 6 . (E.2) 1-e It can be shown that 287 ECR c o o r d i n a t e system: ( x g , y e» z f i) S a t e l l i t e c o o r d i n a t e system: (x , y . z ) 5 d d P l a t f o r m c o o r d i n a t e system: (x , y , z ) FIGURE E . l RELATIVE POSITIONS OF COORDINATE SYSTEMS and R "g R -s 2 + 2 D R 2 - u + D 2 288 (E.3) 2 2 = a / (1 + e' s i n x ) (E.4) Let R = (R , R , R ) -g gx' gy' gz then s i n Y = R / |R I A gz '-g' (E.5) S u b s t i t u t i n g s i n x i n t o Equation E.4 and sim p l i f y i n g gives 2 2 |R I = a - e' R l-g l gz But from Equation E . l R = R + D u gz sz z T T where R = (R , R , R ) and u = (u , u , u ) -s sx sy sz - x y z (E.6) (E.7) S u b s t i t u t i n g Equation E.7 into E.6 and si m p l i f y i n g gives 9 9 9 9 9 |R | = a - e' R - 2 e ' D R u - e ' D u -g gz z z z z (E.8) Combining Equations E.8 and E.3 gives a quadratic equation i n D. The roots to t h i s equation are summarized i n Equations 4.3 to 4.5. E.2 Scan Vector i n Platform Coordinates The two instantaneous viewing angles of a detector are deployment angle oc and scan angle 3 as i l l u s t r a t e d i n Figure E.2. In platform coordinates, the scan vector Du i s = D 1 = sin a cos a s i n 3 cos a cos 3 (E.9) 289 FIGURE E . 3 SENSE OF ATTITUDE ROTATION 290 where 1 is the unit scan vector in platform coordinates. E.3 Scan Vector in Satellite Coordinates The platform has to be rotated through a yaw (K ), pitch ($) and r o l l (to), not necessarily in this order, to align i t s axes with those of the sa t e l l i t e coordinate system. This serves to define the platform coordinates with respect to the sat e l l i t e coordinate system. To conform with photogrammetric convention the order of rotation is taken to be K , $ and and the sense of rotation is shown in Figure E.3. The scan vector Du in terms of s a t e l l i t e coordinates is then X s y s _ z s . = M ps V y p > _ = D ps - (E.10) where u -M p s -1 0 0 Cu> 0 -Su> C<{> 0 0 1 S<)> 0 CK SK -SK CK 0 0 (E.ll) 0 Sw Cm -s<t> 0 C(|> 0 0 1 and where C0 = cos 9 and S9 = sin 6 Since the magnitude of K, <J> and to are small (less than 1° is sufficient for the following approximation), M can be expressed as : 291 M ps CKC$ K ~K CIUCK u> -U) C<t>Cuj (E.12) E.4 Scan Vector In ECR Coordinates Before the scan vector can be expressed in ECR coordinates the sa t e l l i t e coordinate system has to be defined with respect to the ECR coordinate system. The following assumptions are made: (i) The z-axis points towards the earth's centre, ( i i ) The x-axis points in the direction of the sa t e l l i t e velocity vector which is given by the sa t e l l i t e heading, including earth rotation effect, with respect to the geographic north. The transformation between ECR (x , y , z ) and sa t e l l i t e coordinates e e e o (x g, y , z g) is given by: - — - -X X e s = M y J e se z e s + B where M can be derived with the aid of Figure E.4 to se FIGURE E.4 SATELLITE COORDINATES TO ECR COORDINATES 293 M s e (1,1) M s e (1.2) M s e (1,3) M s e (2,1) M s e (2,2) M s e (3,1) M s e (3,2) M s e (3,3) = -SA Sy - Sx C\ Cy = M s e (2,3) M s e (3,1) = -c-x CA = CX Sy - SX Sk Cy = M s e (3,3) M s e (1,1) = -CX SA = CX Cy = M s e (1,3) M s e (2,1) = -Sx M s e (3,3) M s e (2,1) M s e ( 1 ' 3 ) M s e W M s e ( 2 ' 3 ) M s e (E.13) where Y is the sa t e l l i t e heading angle, x is the geometric latitude and X is the longitude of the sa t e l l i t e as shown in Figure E.4. The vector JJ is the displacement vector between the satellite coordinate system and and the ECR coordinate system. 294 APPENDIX F THEORY OF SAR AND SAR DATA PROCESSING This appendix summarizes p r i n c i p l e s of SAR and SAR data processing. F . l Formation of SAR Figure F . l depicts a radar system carried by a spacecraft. The function of the radar i s to obtain a high r e s o l u t i o n map of the t e r r a i n as shown i n the f i g u r e . High range r e s o l u t i o n can be obtained by transmitting pulsed chirp radar signals and matched f i l t e r i n g the return s i g n a l s . Resolution i n azimuth i s i n p r i n c i p l e l i m i t e d by the half power beamwidth (HPBW) of the antenna i n the azimuth d i r e c t i o n . The HPBW at range R i s roughly XR/D where A i s the wavelength of operation and D i s the antenna aperture. Hence the r e s o l u t i o n improves with increase i n aperture. However the physical si z e of the antenna that can be carried on board the spacecraft i s l i m i t e d , and hence l i m i t i n g the azimuth r e s o l u t i o n . Azimuth resolution, higher than the conventional l i m i t of AR/D can be attained by the synthetic aperture technique. If radar pulses are transmitted from a sequence of positions along the f l i g h t path, the return s i g n a l from a target w i l l be a chirp waveform i n the azimuth d i r e c t i o n . If the chirp's frequency modulation (FM) rate i s known, the s i g n a l can be matched f i l t e r e d to achieve high azimuth r e s o l u t i o n . The resolution i s equivalent to synthesizing an aperture that can be thousands of metres long. FIGURE F . l SAR GEOMETRY 296 F.2 Resolution i n Range Consider the case i n which an unmodulated pulse of length T and fixed frequency i s transmitted. The range r e s o l u t i o n i s given by where c i s the v e l o c i t y of propagation. An improvement i n range resolution can be obtained by transmitting a short pulse which requires more transmitter peak power. Equipment peak power l i m i t a t i o n s therefore l i m i t the obtainable range r e s o l u t i o n . The r e s o l u t i o n l i m i t f o r an unmodulated pulse i s limited to the value given by Equation F . l because the beginning and the end of the pulse are i n d i s t i n g u i s h a b l e . If the frequency across the pulse i s made to vary, then the returned pulse at various positions can be distinguished by v i r t u e of the frequency received. This provides the motivation for studying frequency modulated pulses. For a chirp radar, the transmitted s i g n a l i s a l i n e a r FM pulse: q r = c T / 2 ( F . l ) f . ( t ) = r e c t ( t / x) (F.2) where t = time c a r r i e r frequency K = FM rate rect (x) = elsewhere 297 The instantaneous frequency of the pulse is f^ + Kx . Ignoring a multiplicative constant and the rect function the return signal of a target with slant range r is f (t) = cos 2 TT r (F.3) The baseband signal, obtained by beating f r ( t ) with cos ^irf^t and followed by low pass f i l t e r i n g i s f „(t) = cos 2 IT rB (F.4) The matched f i l t e r output, ignoring a multiplicative factor, is a sampling function f = sine [ TT K T (t - 2r/c)] (F.5) r m where sine x = sin x/x. The function f is plotted in Figure F.2. The rm effective pulse width is the -3db width of the main lobe and is given by x = 0.886/(KT) (F.6) err Range resolution i s therefore p v _ ° T e f f _ c 0.886 r ^ 2 T 7 ( F * 7 ) 9 2 For Seasat, K = 562 x 10 /sec and T = 33.9 ysec, hence x = 0.0465 ysec and p r = 7 metres (compare to q^ = 5085 metres). Therefore by transmitting the linear FM pulse, the gain is 730 in range resolution. 298 Amplitude FIGURE F.2 MATCHED FILTER OUTPUT 299 F.3 Azimuth Resolution For s i m p l i c i t y , l e t the transmitted wave form be: f x ( t ) = cos 2 tr f Q t (F.8) then the return s i g n a l i s f (t) = cos 2 77 f _ ( t - 2r/c) (F.9) a u and the baseband s i g n a l i s f a B ( t ) = cos 2 TT f 0 ( 2 r / c ) (F.10) Figure F . l shows the SAR geometry i n which r = / R s 2+ x 2 ( F . l l ) For R » x, r ^ R + x 2 / (2R ) (F.12) s s s Also x = vt (F.13) where v i s the v e l o c i t y of spacecraft. Hence f a B ( t ) = cos 27rf 0[2R s + v 2 t 2 / ( c R g ) ] = cos [ 0 + 2 T r f Q v 2 t 2 / ( c R g ) ] (F.14) where 0 i s a constant phase angle which i s i r r e l e v a n t to the subsequent development. Hence the return i s a l i n e a r FM s i g n a l by v i r t u e of ve h i c l e motion. Comparing t h i s with Equation F.4, the azimuth FM rate i s 2f v 2 K a = -ah (F-13) 2 For Seasat K = 520/sec . The time duration i n which the target stays i n the antenna beam (within the HPBW) i s T = R X/ (D v) (F.16) s s The matched f i l t e r output i s again a sampling function f = sine (TTK T t) (F.17) am a s 300 The plot of t h i s function i s s i m i l a r to that of f i n Figure F.2. The rm function f then shows that a high resolution can be obtained i f the rm target remains illuminated for a longer period of time. The azimuth r e s o l u t i o n , s i m i l i a r i n form to the range resolution given by Equation F.7 i s Pa - SfeSSSj ( F. 1 8 ) a s S u b s t i t u t i n g Equation F.15 and F.16 into F.17 and noting that c =^q » the azimuth r e s o l u t i o n i s P = 0.886 D / 2 (F.19) 3. For the Seasat antenna, D = 10.7 metres. Hence p = 4.74 metres. The a conventional r e s o l u t i o n ( q = X R / D ) i s 18.7 km. The gain i s 3. S therefore 3943 times i n azimuth r e s o l u t i o n . The length of the array synthesized i s 43.2 km. C h a r a c t e r i s t i c s of SAR can be summarized by Equation F.19. They are: (a) Azimuth r e s o l u t i o n i s higher with a shorter antenna. This i s due to the fact that the target stays illuminated for a longer period of time for a shorter antenna (larger HPBW). (b) The azimuth res o l u t i o n i s range independent. Normally the f a r t h e r the sensor i s from the target, the lower the r e s o l u t i o n . However, t h i s i s exactly counter balanced by the fact that the target stays illuminated for a longer period of time when the sensor i s further away. 301 F.4 Two-Dimensional S i g n a l If the transmitted pulse i s f g ( t ) as given by Equation F.2 and i s repeated once every T seconds, then the return w i l l be a two-dimensional (nT and r,. where n i s an integer) signal given by: f r ( n T ) = cos 2 (F.20) Again a m u l t i p l i c a t i v e constant has been ignored and i n the equation assumes an Integral value. The baseband signal i s g B(n) = cos 2 TT 2 V ( T l ) + K/ _ gjrCTO c 2 1 c (F.21) where n has replaced nT which represents time i n the azimuth d i r e c t i o n . T h i s i s a l i n e a r FM s i g n a l i n both d i r e c t i o n s . The matched f i l t e r output i s the product of f and f , i . e . , a two-dimensional sampling function am rirr ' v 6 g (t,n) = sine t TT K x (t - 2r/c)] sine [ ir K T n ] (F.22) in 3 s where K i s again given by Equation F.15. In the d e r i v a t i o n of Equation 3 F.2 i t has been assumed that r ( n ) = y R c 2 + x 2 * R„ + v 2 n 2 /(2 R j s s s i n which the closest approach of the target i s assumed to be at n. = 0. (F.23) An equal phase plot of g„ and g i s shown i n Figure F.3. Therefore D m the fundamental operations to obtain SAR imagery from raw data ( o r 302 FIGURE F.3 TWO DIMENSIONAL SIGNAL AND MATCHED FILTER OUTPUT 303 baseband signal) i s to match f i l t e r the data i n range and azimuth. These operations are c a l l e d range compression and azimuth compression r e s p e c t i v e l y . F.5 The P r i n c i p a l Problem i n S a t e l l i t e SAR Processing [MDA77b] 1 2 2 2 In Equation F . l l , the Pythagorus theorem r = R + x has been invoked, thus assuming that the ground has no curvature. Also, x has been expressed as v and t h i s assumes that the ground i s stationary. However, both of these assumptions are incorrect for any planet. Taking the shape of the earth and i t s r o t a t i o n i n t o account r ^ R + B n 2 / ( 2 R ) (F.24) s s where B i s a function of t r a j e c t o r y data, earth curvature and r o t a t i o n . The azimuth FM rate i s then given by: K = 2 B /(X R ) a s 2 Therefore for a f l a t and non-rotating earth, B = v . The target t r a j e c t o r y , by Equation F.24, i s a parabola and i s plotted i n Figure F.4. The equal phase plot of the baseband s i g n a l i s also shown i n Figure F.4. The range compression i s performed along l i n e s of constant n and produces an output on each n l i n e , a sampling function given by Equation F.5. The sampling function f o r a l l n l i n e s , as shown i n the equation, i s d i s t r i b u t e d along the locus t - 2r(n)/c. As a r e s u l t , the azimuth 1 Sections F.5 and F.6 are abridged from t h i s reference. 304 Azimuth FIGURE F.4 TARGET TRAJECTORY AND EQUAL PHASE PLOT 305 compression must be performed along l i n e s of constant t - 2r(n)/c as shown by the t r a j e c t o r y i n Figure F.4. While t h i s poses l i t t l e problem mathematically, i t can impose considerable implementation d i f f i c u l t y depending upon the desired r e s o l u t i o n and the v e h i c l e / t e r r a i n geometry. In the case of airborne SAR, the maximum duration of t -2r(rj)/c i s normally less than 2 r e s o l u t i o n c e l l s and azimuth compression can, i n p r a c t i c e , be performed along l i n e s of constant range. However, i n the s a t e l l i t e borne case, t h i s deviation can be up to 100 resolution c e l l s . The azimuth compression must then be performed along l i n e s that "wander" through the two-dimensional data. This problem i s referred to as the range c e l l migration (RCM) problem and i s the p r i n c i p a l problem i n s a t e l l i t e borne SAR data processing. F.6 Matched F i l t e r i n g The procedures to process SAR data are range compression, RCM c o r r e c t i o n and azimuth compression. The data i s received a range l i n e at a time and i s immediately range compressed. This can r e a d i l y be done since the FM rate of the transmitted pulse i s a function of system design. Once range compressed, data i s then corner turned so that i t can be read an azimuth l i n e ,at a time. -. The RCM co r r e c t i o n i s then performed by an i n t e r p o l a t i o n process as shown i n Figure F.5. Azimuth matched f i l t e r i n g i s performed on the interpolated data array which e f f e c t i v e l y represents a target point locus. 306 FIGURE F.5 RANGE CELL MIGRATION CORRECTION [MDA77b] 307 The azimuth FM rate i s a function of target p o s i t i o n , s a t e l l i t e p o s i tion and v e l o c i t y . The antenna i l l u m i n a t i o n pattern i s not nec e s s a r i l y exactly broad-side ( i . e . , not exactly right angle to the s a t e l l i t e v e l o c i t y vector) due to a combination of earth r o t a t i o n v e l o c i t y and s a t e l l i t e p i t c h , r o l l and yaw. The beam centre l i n e l i e s along the n = 0 l i n e (time c l o s e s t approach to target) only i n the f i c t i o n a l case of a stationary earth and zero p i t c h and yaw or i n the rare case where the three e f f e c t s cancel. A more t y p i c a l case i s depicted i n Figure F.5. It i s necessary to r e s t r i c t the impulse response duration of the matched f i l t e r as much as possible, confining i t to the illuminated p o s i t i o n only. As shown i n Figure F.6, targets (a) and (b) are matched f i l t e r e d at the illuminated portions of t h e i r t r a j e c t o r i e s and the outputs are placed at locations A and B r e s p e c t i v e l y . This i s frequently referred to as compression to zero Doppler since at n = 0 a l l targets exhibit zero Doppler v e l o c i t y . Benson [BENS80] at OCRS, Ottawa, has approximated the target t r a j e c t o r y by a s t r a i g h t l i n e . This reduces the number of arithmetic operations i n processing the SAR data^ but at the expense of degrading the f i n a l r e s o l u t i o n . The beam centre p o s i t i o n can be obtained by Doppler centroid tracking [MDA79] which estimates the centroid of the antenna pattern from Doppler s p e c t r a l a n a l y s i s . It can be shown that the peak of the mean sp e c t r a l d i s t r i b u t i o n occurs at the Doppler frequency corresponding to the peak of the antenna azimuth pattern. The Doppler frequency i s a l i n e a r function of azimuth distance. 308 FIGURE F.6 A TYPICAL BEAM ORIENTATION 309 F.7 Multi-Look Processing High r e s o l u t i o n SAR imagery su f f e r s from a ' s c i n t i l l a t i o n ' problem. A reduction of the s c i n t i l l a t i o n can be achieved by superimposing images taken from d i f f e r e n t angles to the r e f l e c t i n g surface ( d i f f e r e n t "look" angles). In e f f e c t , a number of coherently generated subimages are incoherently summed to reduce the s c i n t i l l a t i o n i n the f i n a l image. A system of t h i s nature i s c a l l e d a "multi-look" system. The multiple looks can be achieved by separately azimuth compressing the returns from d i f f e r e n t portions of the antenna h o r i z o n t a l beamwidth and incoherently summing the r e s u l t s i n a "post detection i n t e g r a t o r " . This reduces the Doppler bandwidth a v a i l a b l e f o r each look by a m u l t i p l i c a t i v e f actor equal to the inverse of the number of looks. This degrades the f i n a l r e s o l u t i o n but reduces s c i n t i l l a t i o n by a f a c t o r i n the order of /N where N i s the number of looks. Figure F.7 shows the antenna beam i s divided into subbeams i n a 4-look s i t u a t i o n and the method can be generalized to any number of looks. Matched f i l t e r i n g i s done for each subbeam. Let t h e i r outputs be 1^, I^, I^j, where N i s again the number of looks. Then the looks are summed incoherently to give an output: I ( t ) = I j + I 2 + • * ' + 1^ (F.26) If the absolute magnitude of each output i s not taken, the effect i s coherent summation corresponding to 1-look processing. It can be shown that the incoherent summation gives 3 1 0 FIGURE F.7 FOUR-LOOK PROCESSING 311 I(t) i N sine ( T f K T t / N ) (F.27) 3 S Hence the resolution for N looks decreases to N times that for 1-look. For Seasat imagery, the resolution for 1-look is 4.74 metres. Hence for 4-looks i t i s 4.74 x 4 = 19 metres. Due to a sidelobe suppression technique which degrades the resolution 1 somewhat, the fi n a l azimuth resolution is about 25 metres. 1 Range resolution i s 7 metres in slant range and this corresponds to 19 metres in ground range. Sidelobe suppression in the range direction degrades the range resolution to about 25 metres. 312 APPENDIX G CLOSED FORM RECTIFICATION TRANSFORMATION OF SAR IMAGERY This appendix summarizes the geometric rectification transformation of SAR imagery from line, pixel coordinates to map projection coordinates. G.l Slant Range/Azimuth to Map Projection Transformation The process of transforming from line/pixel coordinates to the output map projection coordinates i s performed by a series of discrete transformations, beginning with the line/pixel coordinates in the input slant range image, as follows: u, v Line/pixel in input image coordinates Slant range, azimuth coordinates x , y , z s J s s Sat e l l i t e coordinates Earth centred rotating (ECR) coordinates x, y Map projection coordinates Note that T and t are in units of time. m r The transformations except for input image discrete transformations are the same as coordinates to are specified i n those presented In Appendix E sa t e l l i t e coordinates. These order below. 313 G.2 Input Image to Sa t e l l i t e Coordinate System Given a target's line and pixel coordinates in image memory, i t s slant range and azimuth coordinates can be computed from sampling rates in both directions. If this target has a Doppler frequency f^, the target's coordinates (x g, y , z g) in the sa t e l l i t e coordinate system can be solved from the intersection of three surfaces: Doppler surface m x + m y + m z = f ( f _ ) x s y s z s D Wave front surface (free space, c = constant) (G.l) 2 9 2 2 x + y + z = ( c T 12) s s s m (G.2) Planet surface 2./ us2 2 y s + ( z s -H) = r (G.3) where f ( V m , m , x y m = a function of Doppler frequency, = s a t e l l i t e distance from centre of earth, = planet radius at target location, = velocity of propagation in free space, and = known functions of spacecraft velocity and position. It can be shown that m m = - M dR -s se dt 314 Here M is the rotation matrix from s a t e l l i t e coordinates to ECR se coordinates presented in Equation E . l l and R is the s a t e l l i t e position -s vector in ECR coordinates (see Section E . l ) . In the s a t e l l i t e coordinate system, the z g axis points to the centre of the earth and x g axis lies along the s a t e l l i t e velocity vector. By solving the equations representing the three surfaces, i t can be shown that [f(f_) -mx - m z ] / m D x s z s x 7(CT / v m , v 2 2 . 2) - y - z m s £ t(cx /2) 2 + H 2 - r 2 ] / (2H) m (G.4) The value of Y g i s positive for a radar clock angle of 270° and negative O for a clock angle of 90 . The value of r is target position dependent and requires the use of a planetary datum, plus the planetary coordinates of the target, and incorporates a DTM. G.3 S a t e l l i t e Coordinates to Earth Centred Rotating Coordinate System The coordinates of the target in the ECR coordinate system are given by: x = M se _ s . + B (G.5) where the rotation matrix M shown in Equation E.12 is a function of se s a t e l l i t e position. The vector B is the displacement vector between the s a t e l l i t e coordinate system and the ECR coordinate system. G.4 ECR to Map Projection Transformation The map projection coordinates (x, y) are functions of (x x = g1 (x g, y e, x g) y " g 2 ( V V Z e } The functions g. and g 0 depend on the map projection specified 316 APPENDIX H REFRACTION OF ELECTROMAGNETIC WAVE H.l Introduction This appendix describes a method to trace an electromagnetic wave path when entering a planet's atmosphere, taking refraction into account Refraction i s caused by the difference i n speed of propagation at different altitudes due to change in atmospheric density. The velocity of propagation i s modelled by c' = c Q exp (0z g) (H.l) where c^ is the velocity of propagation at s a t e l l i t e altitude and 3 is the velocity attenuation coefficient due to the planet's atmosphere. Referring to Figure H.l, expressions are derived for z g, x g and P. for a given 6,., and the relation between 9~ and 9 is determined. The values of u u n x and z are used in Section 6.4.3. s s H.2 Expression for z s If the atmosphere is divided into infinitesimally thin layers, the law of refraction states that (Figure H.2): ; sin 9~ sin 6, sin 9., = _ _ i = . . . » _ _ i = . . . (H.2) c 0 c l c i Therefore at any layer i sin 9 i = (Cj / c Q) sin QQ (H.3) and vertical displacement is given by: dz = (dr) cos 9. 317 FIGURE H.2 REFRACTION IN A THIN LAYER 318 But dr = c^dt, ther e f o r e dz = c, cos 0 , dt s i i From Equations H.l and H.2, c o s 6 1 =Jl - [ci I c Q ) s i n 6 Q ] 2 Therefore, =Jl - exp(-'2Bz g) s i n 2 6 Q (H.4) g z g ) ^ 1 - exp(-2gz g) s i n 2 e Q dt d z g = c Q exp(- ^) and dt = exp(Bz ) dz s s so that t = C0J 1 ~ exp(-2gz g) s i n 0 Q 2 s exp(3 z ) dz s s Q c Q J l - exp(-23z g) s i n 2 e 0 = exp(Bz g) Jl - exp (-2Bz g) s i n e Q - cos e Q (H.5) 1 2 2 Hence z = — In [ ( t c. + cos fl_) + s i n o.] (H.6) s 28 0 0 - 0 H.3 Expression f o r x g R e f e r r i n g t o Figure H.2, the l a t e r a l displacement i n l a y e r i i s dx = (dz ) tan 6 s s 1 From Equations H.l and H.2 and Figure H.2, i t can be shown that dx = (dz ) tan { s i n * [exp(-gz ) s i n 8 n]} S S S \) = WO exp(-Bz g) s i n e Q s -1 cos { s i n [exp(-Bz ) s i n 6 n]} s u 319 exp(-gz ) sin Jl - exp(-Bzs) sin 0 Q Therefore, x = s s exp(-gzg) sin 0^ 0 ^1 - exp ( - 3 z g ) sin 6Q = | {eQ - s i n " 1 (exp(-Bz g) sin 0Q]} (H.6) H.A Expression for R Referring to Figure H.2, the path length is dR = dz / cos 6. s 1 dz dR = Jl - exp(-26zg) s i n y 0 Q Therefore, s dz R = • Q y i - exp ( - 2 3 z g ) sin QQ i i + JTZ 2 exp(-2gz ) sin 6. = j l n S., exp(Bz ) (H.7) 3 / T + 1 - sin* e 0 s H.5 Relation Between 0 and 0. n U Referring to Figure H.2, sin 0 = (c /c n) sin 0 n = exp(~3z ) sin 0_ n n (J u s u From Equation B.7, 2 2 exp(gz g) = (t c 0 + cos e Q) + sin 0 Q 2 2 = 2 t 3 c Q cos 0 Q + 1 + t 3 c Q Therefore, 2 2 7 sin 0 n =(2 t 3 c Q cos 0 ( ) + 1 + t g c Q) s i n ^ 0 q PUBLICATIONS R. Orth, F.Wong and J.S. MacDonald, "The P r o d u c t i o n o f 1:250,000 Maps o f P r e c i s i o n R e c t i f i e d and R e g i s t e r e d L a n d s a t Imagery U s i n g t h e MDA Image A n a l y s i s System: I n i t i a l R e s u l t s " , T w e l f t h I n t e r n a t i o n a l Symposium on Remote S e n s i n g o f Envir o n m e n t , M a n i l a , P h i l i p p i n e s , A p r i l 1978. F.Wong and R.Orth, " R e g i s t r a t i o n o f S e a s a t / L a n d s a t Composite Images t o UTM C o o r d i n a t e s " , S i x t h Canadian Symposium on Remote S e n s i n g , H a l i f a x , N.S., May 1980. F. Wong, D.E. Friedmann and R. Orth, "The Use o f D i g i t a l T e r r a i n Model i n t h e R e c t i f i c a t i o n o f S a t e l l i t e Borne Imagery", F i f t e e n t h I n t e r n a t i o n a l Symposium on Remote S e n s i n g o f Envi r o n m e n t , Ann A r b o r , M i c h i g a n , May 1981. F. Wong and R. Orth, " O r b i t Refinement f o r t h e R e c t i f i c a t i o n o f SAR Imagery", ISPRS, Symposium o f Commission I I I , H e l s i n k i , F i n l a n d , June 1982. F. Wong and R. Orth, " R e c t i f i c a t i o n o f TIROS-N Imagery w i t h R e c u r s i v e F i l t e r i n g T e c h n i q u e f o r T i m e - S e r i e s E s t i m a t o r " , ISPRS, Symposium o f Commission I I I , H e l s i n k i , F i n l a n d , June 1982. F. Wong, "Resampling of Landsat-D Thematic Mapper Imagery in the Presence of Scan Gaps" In Preparation
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A unified approach to the geometric rectification of...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A unified approach to the geometric rectification of remotely sensed imagery Wong, Frank Hay-chee 1984
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | A unified approach to the geometric rectification of remotely sensed imagery |
Creator |
Wong, Frank Hay-chee |
Publisher | University of British Columbia |
Date Issued | 1984 |
Description | Many applications of remotely sensed digital imagery require images that are corrected for geometric distortions. It is often desirable to rectify different types of satellite imagery to a common data base. A high throughput rectification system is required for commercial application. High geometric and radiometric precision must be maintained. The thesis has accomplished the following tasks: 1. The sensors used to obtain remotely sensed imagery have been investigated and the associated geometric distortions inherent with each sensor are identified. 2. The transformation between image coordinates and datum coordinates has been determined and the values of the parameters in the transformation are estimated. 3. A unified rectification approach has been developed, for all types of remotely sensed digital imagery, which yields a high system throughput. 4. Use of digital terrain models in the rectification process to correct for relief displacement has been incorporated. 5. An efficient image interpolation algorithm has been developed. This algorithm takes into account the fact that imagery does not always correspond to sampling on a uniform grid. 6. The applications of rectified imagery such as mosaicking and multisensor integration have been studied. 7. Extension of the rectification algorithm to a future planetary mission has been investigated. The sensors studied include TIROS-N, Landsat-1, -2 and -3 multispectral scanners, Seasat synthetic aperture radar, Landsat-4 thematic mapper and SPOT linear array detectors. Imagery from the last two sensors is simulated. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-06-14 |
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.0051914 |
URI | http://hdl.handle.net/2429/25697 |
Degree |
Doctor of Philosophy - PhD |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1984_A1 W65.pdf [ 36.03MB ]
- Metadata
- JSON: 831-1.0051914.json
- JSON-LD: 831-1.0051914-ld.json
- RDF/XML (Pretty): 831-1.0051914-rdf.xml
- RDF/JSON: 831-1.0051914-rdf.json
- Turtle: 831-1.0051914-turtle.txt
- N-Triples: 831-1.0051914-rdf-ntriples.txt
- Original Record: 831-1.0051914-source.json
- Full Text
- 831-1.0051914-fulltext.txt
- Citation
- 831-1.0051914.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}]}"
data-media="{[{embed.selectedMedia}]}"
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-0051914/manifest