Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Motion compensation for airborne interferometric synthetic aperture radar Stevens, David Robert 1994

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1994-0148.pdf [ 3.41MB ]
Metadata
JSON: 831-1.0065119.json
JSON-LD: 831-1.0065119-ld.json
RDF/XML (Pretty): 831-1.0065119-rdf.xml
RDF/JSON: 831-1.0065119-rdf.json
Turtle: 831-1.0065119-turtle.txt
N-Triples: 831-1.0065119-rdf-ntriples.txt
Original Record: 831-1.0065119-source.json
Full Text
831-1.0065119-fulltext.txt
Citation
831-1.0065119.ris

Full Text

Motion Compensation for AirborneInterferometric Synthetic Aperture RadarbyDAVID ROBERT STEVENSB.A.Sc., Engineering Physics, University of British Columbia, 1989A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE STUDIESTHE DEPARTMENT OF ELECTRICAL ENGINEERINGWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAMarch 1994© David Robert Stevens, 1994In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)__________________________Department of l.QC1Y)ttA-( EL?eers49The University of British ColumbiaVancouver, CanadaDate 1rO4 3DE-6 (2188)AbstractAirborne SAR interferometry has the potential to provide topographic data with a precisionof the order of one meter. However, to generate data accurate to this level it is essential tomeasure and compensate for the antenna baseline motion. Conventional motion compensationtechniques and their errors are analyzed and extended to the two channel simultaneous imagingscenario of InSAR. An evaluation of the modelling is made using point target simulation andreal motion and InSAR data. Phase compensation of both channels to the same reference trackand compensation to two separate tracks are considered. The single track approach allowstrack segmentation to follow aircraft drifts without causing discontinuities in the differentialphase, but is sensitive to range cell migration effects. The dual track approach is not sensitiveto this but suffers from discontinuous differential phase at segmentation boundaries, whichcomplicates the phase unwrapping process. A new formulation for each approach is presentedthat compensates for unknown terrain coupled with low frequency aircraft motion. In addition,a new approach that uses the dual track approach initially and then converts to a single referencetrack after compression is proposed. This realizes the benefits of both approaches with onlya small increase in computation.11ContentsAbstract . .List of Tables. . .List of Figures . . .. . . . . . . ...Acronyms, Abbreviations and Symbols .AcknowledgmentsChapter 1 Introduction1.1 Topographic Mapping With Radar1.2 Topographic Mapping with Interterometric1.2.1 Spaceborne Interferometric SAR .1.2.2 Airborne Interferometric SAR1.3 Thesis Outline and ContributionsChapter 2 Airborne InSAR Background2.1 Conventional SAR2.1.1 Conventional Motion Compensation2.2 Outline of Airborne lnterferometric SAR2.2.1 Basic Processing Steps2.2.2 Airborne InSAR Motion Compensation• . . . . viiixixiv11Radar . . 223477810101114141515Chapter 3 Single Channel Airborne Motion Compensation.3.1 Assumptions3.2 Ideal Motion Compensation3.2.1 Point Target FormulationifiIi• • vii3.3 Non-Ideal Motion Compensation.3.3.1 Zero Doppler Pulse Processing3.3.2 Effects of Unknown Terrain3.3.3 Effects of No Motion Compensation Resampling . .3.3.4 Extension to Clutter Data3.3.5 Effects of Inertial Data Errors3.3.6 Summary of Single Channel Errors and ConstraintsChapter 4 Motion Compensation for Interferometry4.1 Assumptions4.2 InSAR Motion Compensation Requirements4.3 Ideal Motion Compensation4.3.1 Dual Reference Track Assuming Ideal Compensation4.4 Non-Ideal Motion Compensation4.4.1 Single Reference Track with Unknown Terrain4.4.2 Dual Reference Track with Unknown Terrain4.4.3 Effects of Unknown Terrain4.4.4 FM Rate Induced Differential Phase Errors4.4.5 RCMC Induced Differential Phase Errors4.4.6 Effects of No Motion Compensation Resampling4.4,7 Effects of Inertial Data Errors4.4.8 Summary of Errors and Constraints for InSAR4.5 Proposed Motion Compensation Approaches4.5.1 Realization of Requirementsiv22222833343540434344474849505456646668697074744.5.2 Combination of Single and Dual Reference TracksChapter 5 Evaluation of Non-Ideal InSAR Motion Compensation.767979798181879699100105108110111. . . . . 112Appendix B Single Channel Phase Errors due to Unknown TerrainC Effect of Linear Phase Errors on CompressionD Directional Random Walk ModelE Inter-channel Azimuth BroadeningF Single Reference Track Phase ErrorsG Dual Reference Track Phase Errors5.1 Point Target Simulation5.1.1 The Point Target Simulator5.1.2 Modelling of Typical Aircraft Motions.5.1.3 Simulation Results of Modelled Motions.5.1.4 Simulation Results of Actual Motion Data5.1.5 Discussion of Point Target Results . .5.2 Processing CCRS InSAR Data. . . . .5.2.1 Three Hills Data Results5.2.2 Inertial Data ErrorsChapter 6 Conclusions6.1 Recommended Approach for CCRS System6.2 Recommended General ApproachChapter 7 Suggestions for Future Work . . . . .Bibliography . . . . . . . . . . . . . . .Appendix A Effect of Quadratic Phase Errors on CompressionAppendixAppendixAppendixAppendixAppendix113116119121122125127129VAppendix H Unknown Terrain and Translational Motion 132Appendix I Motion Compensation Resampling. . . . . . . . . . . 133Appendix J MATLAB Program Listings 135viList of TablesTable 1 Typical CCRS Parameters 21Table 2 Single Channel Phase Errors 41Table 3 Single Channel Motion Compensation Limitations (CCRSTypical) 42Table 4 InSAR Motion Compensation Errors and Limitations #1 . . 72Table 5 InSAR Motion Compensation Errors and Limitations #2 . . . . 73Table 6 Comparison of Single and Dual Reference Track MotionCompensation 74Table 7 Typical Flight Motion for CCRS Convair 580 81Table 8 Dual Reference Track Motion Compensation with RCMC . . . 83Table 9 Effects of Angular Acceleration, aroll = 0.3 , h = 1 km, andno antenna pattern 85ViiList of FiguresFigure 1Figure 2Figure 3Figure 4Figure 5Figure 6Figure 7Figure 8Figure 9Figure 10Figure 11Figure 12Figure 13Figure 14Figure 15Figure 16Figure 17Figure 18Figure 19Figure 20• . 12• . 17• . 19• . 2325293435363649508282Synthetic Aperture Radar Imaging GeometryCCRS InSAR Processing StepsSingle Target Motion Compensation GeometryZero Doppler Applied Motion Compensation Cross-sectionMultiple Target Motion Compensation GeometryConstant Offset Flight MotionMotion Compensation to a Reference Track.Flight Profile Without ResamplingFlight Profile With ResamplingNeglecting Motion Compensation ResamplingCross-section of Neglecting Resampling GeometryDouble Reference Track GeometrySingle Reference Track GeometryConvair 580 Flight Displacement DataConvair 580 Roll VariationUnknown Terrain with = 0.5 , and R = 10 km 84Unknown Terrain with ajj03 = O.Olg, and R = 10 km 85Effect of Sinusoidal Roll on Differential Phase (0.15 deg peak,period = 3 see, h = 1 km, R = 10 km) 86Effect of Sinusoidal Roll on Height Estimate (0.15 deg peak,period=3see, h=lkm, R=lOkm) 87Single Track Bias Errors, R = 10 km, and no antenna pattern. 88VuFigure 21 Measured Versus Expected Differential Phase Dual Track(R 10 km, h = 1 km) 89Figure 22 Height Estimation Errors and Roll (R = 10 km, h 1 kin) . . . 90Figure 23 Height Estimation Errors (R = 10 km, h = 1 km) 91Figure 24 Height Estimation Errors (R = 15 kin, h = 1 km) 91Figure 25 Height Estimation Errors (R = 20 km, h = 1 km) 92Figure 26 Azimuth Broadening with Unknown Terrain (R = 10 km,h=lkm) 93Figure 27 Azimuth Broadening with Unknown Terrain (R = 10 km,h=500m) 94Figure 28 Azimuth Shift with Unknown Terrain (h = 1 km) 95Figure 29 Measured Versus Expected Differential Phase Single Track(R = 10 km, h = 1 km) 96Figure 30 Height Estimation Errors Single Track no RCMC (R = 10 kin,h=lkm) 97Figure 31 Differential Phase with Segmented Reference Tracks(R=lOkm,h=lkm) 98Figure 32 Height Estimates with Segmented Reference Tracks(R=lOkm,h=lkm) 98Figure 33 Square Root of Interferogram Magnitude Image 101Figure 34 Dual Reference Track Interferogram Phase 102Figure 35 Dual/Single Reference Track Interterogram Phase 103Figure 36 Single Reference Track Interferogram Phase 104Figure 37 Elevation Estimates 105ixFigure 38 1:50 000 Contour Map of 3 Hills Area (51°45’ lat. 113°30’ long.),25 ft contours (82 P/11/12 Energy, Mines, and ResourcesCanada, 1990) 106Figure 39 Differential Resampling Geometry 133xAcronyms, Abbreviations and SymbolsRADAR RAdio Detection And RangingSAR Synthetic Aperture RadarInSAR Interferometric Synthetic Aperture RadarDEM Digital Elevation ModelCCRS Canada Centre for Remote SensingRCM Range Cell MigrationRCMC Range Cell Migration CorrectionGPS Global Positioning SystemJPL Jet Propulsion LaboratoryRVPC Range Varying Phase CorrectionINS Inertial Navigation SystemFM Frequency ModulationI In-phase ChannelQ Quadrature ChannelPSF Point Spread FunctionPDF Probability Density FunctionFIR Finite Impulse ResponseChannel A Transmit and Receive Radar ChannelChannel B Receive Only Radar ChannelRT Reference TrackISLR Integrated Side Lobe RatioR Slant rangeV,v Velocityt Azimuth time variableT Processed aperture time) WavelengthPRF Pulse Repitition Frequencyç& PhaseF Differential phased Displacement from reference trackxi(9 Off nadir angleH Altitude of reference track (channel A)los Line of sight componentI los Perpendicular to los componentf0 Range spectral shiftc Speed of lightFr Range sampling rateBWr Range Bandwidth°r Range Oversampling FactorK FM rateYa Altitude of transmit/receive antenna/Rmf Error in range used for matched filter generationh Error in assume terrain elevationLP Azimuth shift in position of compressed peaka accelerationg acceleration due to gravityLRINS Inertial Data position erroro- Phase variance of point target peako- Variance of phase noiseN Number of independent samplesy Inter-channel correlation coefficient4’ Differential phaseVariance of differential phaseaab Antenna baseline angleaRT Reference track baseline angleb Antenna baseline separationbRT Reference track baseline separationA Antenna A aperture center positionB Antenna B aperture center positionE(A) Error in applied phase correction for channel AE(B) Error in applied phase correction for channel BLaab Variation of roll angle about mean valuearoll Angular acceleration of antenna baselinexuD Length of antenna along trackM Displacement from aperture centerbias Differential Phase Bias/hbias Bias error in height estimate/&Rab Inter-channel difference in slant range matched filter errorszd Inter-channel difference in dxliiAcknowledgmentsI would like to extend my sincere thanks to my supervisors Dr. Ian Cumming and Dr. M. Itofor their enthusiastic and instructive supervision and financial support. A special thanks is alsoextended to Dr. Laurence Gray of CCRS for his many enlightening contributions, discussions,financial support, and for supplying the data used in this work which was prepared by PeterFarris-Manning. In addition, I would like to thank my fellow students Gordon Davidson,Michael Seymour, and Yong Luo for their moral and technical support. I would also like tothank Paul Lim of MacDonald Dettwiler for getting me interested in SAR in the first place.Last but not least, I would like to thank my wife Kathy and my daughter Michelle for theirconstant support and understanding. Financial support for this project was generously providedby scholarships from the Natural Sciences and Engineering Council, MacDonald Dettwiler, theB.C. Science Council, and UBC. Additional support came from contracts with CCRS. All dataused in this study was provided generously by CCRS.xivChapter 1 IntroductionThe knowledge of the Earth’s topography has been of practical importance for all civilization; initially locally and then gradually to the global scale. For example, the design of anancient Roman aqueduct would need to consider the local topography in some detail. Untilrecently, the consolidation of terrain features and topography into a map has been a laboriousprocess of collating and analyzing information from ground surveys. With the advent of theaviation age aerial photographs could be used in a stereoscopic mode to estimate terrain elevation and identify terrain features. Since the introduction of the computer, attempts have beenmade to automate this arduous process leading to many new techniques of mapping terrainfeatures and topography.1.1 Topographic Mapping With RadarOne such new technique uses an active microwave sensor (radar) mounted on an airplane oron a satellite. Radar maps are now available in which high-resolution Synthetic Aperture Radar(SAR) images are used to portray the terrain features [1—3]. The radar sensor provides coveragein all weather, and images the terrain in a unique way which has certain advantages to users suchas geologists and hydrologists [4]. In addition, this airborne and spaceborne imaging technologyhas been extended to topographic mapping through the processing of multiple images.In the initial method of stereo radargrammeiry, two passes of airborne SAR data areobtained in an aerial survey, and a correlation process is used to extract terrain height bymanipulating the parallax between the two generated images [5]. However, radar specklenoise, caused by coherent scattering effects, and geometric distortion, limit the accuracy ofthe correlation.11.2 Topographic Mapping with Interferometric RadarA new technology, Interferometric SAR (InSAR), is being developed in which elevationdata is obtained by processing a pair of stereo images in a unique way [6, 7]. This processutilizes the fact that SAR is a coherent imaging system that is able to produce a complex valuedimage. If the processing is carefully done, the phase of a given image pixel is a function of thepath length travelled by the radar pulse coupled with scene scattering effects. The magnitudeof each image pixel is a function of the reflectivity of the corresponding patch in the scene.When the two receiving channels are processed together, interferometric phase patterns areobtained that are caused by the height differences in the imaged terrain. The phase patternsobtained in this “interferogram” are then analyzed by a computer to obtain the elevation ofeach pixel in the scene by triangulation. The resulting digital elevation model (DEM) can bemore accurate than those previously available with SAR.12.1 Spaceborne Interferometric SARThe spaceborne implementation of this approach combines two SAR images of the sameregion from two separate satellite passes called satellite InSAR [8]. The spaceborne implementation has certain advantages due to the predictable platform motion and large groundcoverage. The disadvantages are:• The image signal to noise ratio is not very good.• There may be changes in the ground reflectivity between the particular orbits used, leadingto temporal decorrelation between the two images.• The two images from different orbits must be accurately registered to each other (the orbitpositions are difficult to determine accurately).• The inter-channel baseline separation, which determines the parallax, is a function of theoffset between the two orbits and therefore is not easily adjustable.2The baseline separation, or parallax, between the two images, and the radar frequencydetermine the interferometer’ s sensitivity to terrain elevation and also determines the amountof baseline speckle noise (decorrelation caused by the parallax). For typical satellite altitudes,radar frequencies, and for adequate terrain elevation sensitivity, the baseline needs to be severalhundred meters [9]. It is therefore impractical to produce a DEM from a single pass of a satelliteby mounting two typical radar antennas on a single satellite.Several alternate approaches to obtain single pass InSAR for the satellite case have beenproposed. One approach is to use a very short wavelength allowing a much smaller baseline[10], and another is to used a tethered satellite system [11]. With the recent success of theGlobal Positioning System (GPS) it has been proposed that very accurate sateffite interferometrycould be done by having two satellites travelling in nearly parallel orbits making simultaneousimagery [12]. By using GPS in a differential mode it should be possible to estimate the relativeorbits very accurately. The simultaneous imaging removes the temporal decorrelation inherentto conventional satellite techniques. However, each of these possible improvements to satelliteinterferometry pose severe technical challenges, and remain impractical for the present time.1.2.2 Airborne Interferometric SARPerhaps the most promising approach to radar topographic mapping in the short term isin using a single flight pass of a dual antenna airborne SAR [7]. The InSAR system uses aconventional SAR for transmitting and receiving the radar signal, but in addition, a secondantenna and receiving channel is used to obtain another registered image displaced acrosstrack. This single pass approach is possible because the flight altitude is much smaller thanthat required for the satellite orbit, thus requiring only a several meter baseline. There areseveral advantages to the airborne approach:• The altitude of the flight can be adjusted to provide the desired elevation estimation3sensitivity.• The signal to noise ratio is typically quite large (10—30 dB).• The inter-channel registration is simpler due to the rigid structure connecting the twoantennas.• There is great flexibility in the direction and number of passes made of a particular region,which can be useful if radar layover and radar shadow are present in rugged terrain regionsleading to regions of insufficient sampling of the terrain,• There is no temporal decorrelation between the two images due to changing groundscattering conditions because the two images are made simultaneously.The main disadvantage of the airborne approach is that the trajectory of the antennas issubject to perturbations due to air turbulence and aircraft flexure. The effect of undesiredaircraft flight motion on the InSAR system is very complex and not well understood.1.3 Thesis Outline and ContributionsThe primary focus of this thesis is to understand the role that motion compensation plays inairborne InSAR. More specifically, we wish to establish the fundamental limitations that motioncompensation places on InSAR. It is necessary to first understand the mechanisms involvedin conventional airborne SAR before an analysis of interferometric SAR can be made. Giventhat the purpose of InSAR is to estimate topography, the effect of unknown terrain on motioncompensation will be dealt with in some detail. It will be shown that there are fundamentallimitations on the accuracy of motion compensation even when the aircraft motion is knownprecisely. Fortunately, most errors due to motion compensation are similar in both InSARchannels and thus cancel out, leading to small height estimation errors. Some of the residualerrors can be minimized by following new approaches, for instance, by combining the best4aspects of two existing methods of motion compensation for InS AR. The thesis report is laidout as follows:• Chapter 2 will be devoted to the background for airborne InSAR. This will include anoverview of Synthetic Aperture Radar (SAR) as well as the basic processing involvedin airborne InS AR, The chapter concludes by introducing the issues involved in motioncompensation for airborne InSAR.• Chapter 3 will deal with single channel airborne motion compensation. Initially, idealmotion compensation will be discussed which assumes that all geometric parametersneeded for the compensation are known and infinite processing resources are available.Following this, the realistic situation of non-ideal motion compensation will be considered.Several new formulations will be presented such as the coupling effect between the rangevarying phase correction (RVPC) and range cell migration correction (RCMC), the effectsof unknown terrain, and the directional random walk model for phase smoothing of pointtargets by azimuth compression.• This single channel analysis is then extended to a comprehensive chapter on motioncompensation for interferometry. The bulk of Chapter 4 is original work as little has beenpublished in this area. As in Chapter 3, ideal compensation and non-ideal compensationwill be considered. The highlights of the new work include a requirements analysis forInSAR motion compensation, a detailed comparison of the single and dual reference trackapproaches, closed form solutions for estimating terrain elevation for each approach whichremoves some terrain induced bias errors, an analysis of differential effects caused byunknown terrain and differential errors caused by the coupling between the range varyingphase correction (RVPC) and range cell migration correction (RCMC), and the proposal ofa new approach that combines the best aspects of the single and dual track methods.• The validity of the presented theory and modelling is verified in Chapter 5 by experimen5tation. An InS AR point target data generator and processor is used with modelled andreal aircraft flight motions. Good agreement between the theory and simulation results isdemonstrated. Real InSAR data from CCRS is also processed by a basic InSAR processorusing the single and dual track approaches. The results illustrate the differences in theeffects of the various motion compensation approaches on the differential phase of theinterferograms. A slant range topographic map is generated and compared visually witha standard 1:50 000 contour map.• Chapter 6 presents a summary of the conclusions from this work. This includes recommendations for the CCRS InSAR motion compensation processing as well as a more generalapproach.• Chapter 7 lists suggestions for further study such as removing some of the assumptionsmade in the current work and extending the analysis to other operational scenarios.6Chapter 2 Airborne InSAR BackgroundIn this chapter the general background required to understand airborne InSAR will bepresented. This will begin with an overview of conventional SAR followed by the extensionto airborne InSAR.2.1 Conventional SARConventional Synthetic Aperture Radar (SAR) is itself a very interesting and challengingfield [13, 3]. For the purposes of interferometry, only topics in conventional SAR that arerequired to understand the topographic mapping problem will be considered. The basic idea ofSAR is that a long antenna is synthesized in the flight direction by transmitting and receivingpulses of microwaves as the much smaller real antenna passes over a region (Figure 1). Thebackscatter echo from each pulse, incident on a scene patch, is processed coherently to obtaina high resolution reflectivity value in the along flight or “azimuth” direction. In the acrosstrack or “range” direction pulse compression techniques are used to obtain good resolution inthis coordinate. The range to the target is obtained from the radar timing.A simplified view of conventional SAR is that the 2—D raw data received in range andazimuth space can be focussed into a reflectivity image by a 2—D linear compression operation.The range coordinate is actually the radar line-of-sight direction and therefore a conversionto “ground” range is necessary for mapping purposes. One complication in the processingcomes from the fact that as the synthesized aperture is formed the slant range to the targetchanges. This causes the echo data from a given target to migrate across range called RangeCell Migration (RCM). In order to perform the compression in the azimuth direction the echodata must be aligned with the flight direction. This operation is called Range Cell MigrationCorrection (RCMC). For many airborne applications RCMC can be ignored because only the7central part of the aperture is used so the range migration is small. But, when very highresolution in azimuth is required the entire aperture must be processed and often RCMC isnecessary. Another complication which is specific to the airborne implementation is the needfor motion compensation.Figure 1 Synthetic Aperture Radar Imaging Geometry2.1.1 Conventional Motion CompensationThe need for motion compensation to achieve accurate focusing of airborne SAR images haslong been understood [14]. Although, efficient approximations to the ideal compensation arestill being developed and evaluated. The basic approach involves using the displacement historyof the antenna phase center with respect to a selected reference track to calculate a range gate(or time) delay adjustment for each pulse and a slant range dependent phase correction for eachbeam footprint8range-compressed pixel. This correction extends along the entire synthesized aperture and isapplied to the range compressed data prior to azimuth compression. This can be accomplishedin real-time on the aircraft or as a post-ifight processing step. The compensated data thenemulates the transmission and reception from an antenna travelling in a straight line along thereference track provided the compensation has been successful [15].Accurate phase correction requires knowledge of the complete geometry of the imagingsituation. This requires knowledge of the elevation angle to the scene patch of each range-compressed pixel as well as the precise position of the antenna. Some work has been doneon estimating the flight trajectory from analyzing the data post-flight, for example [16], butmany approaches rely on Inertial Navigation Systems (INS) to provide the antenna positioninformation [17]. Estimating the elevation angle is a more difficult problem.Calculation of the elevation angle for each image pixel requires a prior estimate of theterrain height. Without access to a Digital Elevation Model (DEM) of the imaged region, someassumptions or estimates must be made. Typically, the terrain is assumed to be flat at somereference elevation, but the effect of this assumption has not been studied in detail. Mostmotion compensation error analysis has considered uncompensated motion, such as INS errorsand their effect on image focus [17], without addressing the effect of terrain height on motioncompensation,When the residual phase errors after motion compensation, due to INS errors and unknown target elevation, are small, subsequent azimuth compression can be accomplished ina straightforward manner and defocus and image distortion are minimized. As conventionalSAR is typically concerned with a magnitude image not the phase, little is known about of theeffect of aircraft motion on the image phase. In order to evaluate the importance of variousmotion compensation errors for terrain mapping an overall understanding of airborne Interferometric SAR processing is required. In addition, it is possible to compensate for some motion9compensation errors at later processing stages.2.2 Outline of Airborne Interferometric SARAirborne InSAR is basically the interferometric processing of two channels of airborneSAR data. The individual channel SAR processing algorithms and subsequent interferometricprocessing will depend upon the radar hardware used to obtain and possibly process the receiveddata. In this section, the basic processing steps for a typical C-band airborne InSAR systemwill be discussed followed by an overview of InSAR motion compensation.2.2.1 Basic Processing StepsThere are many processing steps involved in a typical InSAR system (Figure 2). Thesteps described below follow those used in the Canada Centre for Remote Sensing (CCRS)InSAR system [18]. Other approaches change the order of some steps and may use alternatealgorithms to carry out a particular function.Before the data is digitized the received echoes are first compressed in the range directionby applying a matched filter operation. This step is often performed in the analog electronicsof the receiver system for each of the two channels. The following steps are performed bydigital computation within a computer:• The two channels are digitized and then registered to each other across range to ensurethat each corresponding pixel in the two channels are correlated.• Motion compensation is then applied to both channels.• The motion compensated range compressed data is then compressed in the azimuth directionby application of a matched filter.• The interferogram is then formed by applying a complex conjugate multiplication, pixel bypixel, on the two conventional SAR images. The phase of the interferogram is thereforethe difference in the phases of the two images and is often called the differential phase.10• The differential phase of the interferogram can be smoothed by sample averaging in theazimuth and range directions to reduce phase noise.• The modulo 2ir differential phase must be unwrapped by adding some multiple of 27r toeach sample in order to obtain the absolute differential phase which represents the differencein the path length phase between the two antennas.• This absolute phase can then be converted into an elevation estimate from the imaginggeometry.• The scene patch location estimate in the SAR reference frame must finally be transformedinto a convenient mapping coordinate system.For airborne InSAR the two steps that pose the greatest challenges are motion compensationand phase unwrapping. The phase unwrapping problem is common to satellite InSAR, as wellas other applications in signal and optical processing, and has been studied in detail but arobust solution for InSAR has yet to be found [19—21]. The motion compensation problemwill now be considered.2.2.2 Airborne InSAR Motion CompensationThe InSAR motion compensation situation is more complicated than in conventional SARbecause two antennas are involved and terrain height estimation requires high phase accuracyin the processed images. Inter-channel phase bias errors lead directly to biases in the derivedterrain model and must therefore be minimized.Two approaches to InSAR motion compensation have been demonstrated, but neither hasbeen analyzed in detail. The first approach is that adopted by the Jet Propulsion Laboratory(JPL) group. This method involves defining two reference tracks, one for each antenna, andthen applying motion compensation for each channel separately [7, 22, 21, 23]. In the JPLcase the reference tracks are segmented to keep the displacements small to allow for various11Channel A (main antenna) Channel B (insar)Resample to Channel ACell CentresMotion Compensation Motion CompensationAzimuth Compression Azimuth Compressionlnteerogram Generationlnterterogram SmoothingAbsolute Phase UnwrappingPhase to ElevationConversionSAR Geometry toMapping GridFigure 2 CCRS InSAR Processing Stepsapproximations to be made. The absolute phase of each segment is estimated by using theSplit-Frequency Approach [21]. In addition, the antenna is not yaw steered so processing tozero Doppler can not be used.The second approach is used by the Canada Centre for Remote Sensing (CCRS). Thismethod motion compensates both antennas to the same reference track [18]. This approach12basically attempts to undo the relative mapping between the two images caused by the parallax.Any residual mapping is attributed to the unknown terrain effects. The basic processing to datehas not used motion compensation resampling and relies upon control points to resolve theabsolute phase ambiguity. In addition, the antenna is yaw steered, so zero Doppler processingcan be used.Motion compensation is a complex component of InSAR processing due to the difficultyin predicting the effects of inertial data errors, unknown terrain elevation effects, and possibleside-effects from the motion compensation operation. In the next chapter a detailed analysisof single channel motion compensation will be made, specifically considering effects that areimportant to interferometry. Next, the InSAR motion compensation problem will be definedand studied in detail. The primary errors sources will be identified and new approaches tominimizing their adverse effects will be presented. The developed theory wifi be verified bysimulation of the InSAR imaging and processing situation using modelled and real motiondata. A simple combination of the single and dual reference track approaches will be proposedthat combines the benefits of both. This new processing approach will be demonstrated byprocessing actual C-band InSAR data from the CCRS InSAR system.13Chapter 3 Single Channel AirborneMotion CompensationMotion compensation for single channel SAR systems is a complex problem of attemptingto compensate for non-linear motion of the aircraft which occurred during the aperture synthesis.Complete compensation considers the yaw, pitch, roll, and translational motion of the aircraft.In order to focus on particular effects that are important for InSAR a number of reasonableassumptions can be made which will simplify the analysis and allow certain error sources tobe isolated.31 AssumptionsThe following is a list of the main assumptions made for the subsequent analysis:1. The antenna is yaw steered to zero Doppler. This means that the antenna is assumed to bepointing perpendicular to the velocity vector of the aircraft.2. The SAR processing is to zero Doppler.3. The pitch of the aircraft is fixed.4. The spacing between transmitted pulses on the ground is kept fixed by adjusting the PulseRepetition Frequency (PRF) to be proportional to the aircraft velocity.5. The inertial data will be assumed to be correct, for the most part.6. The motion compensation to be applied will be post-flight, apart from those mentionedin the above list.7. Translational motion and aircraft roll must be compensated through motion compensation.For a typical airborne system, such as the CCRS system, most of these assumptions arereasonable. Errors in the yaw and pitch of the aircraft will lead to radiometric problems but notphase problems if processing is to zero Doppler. The main limiting factor is the assumption14that the inertial data is correct. When inertial data errors occur this can lead to defocussingand mis-registration. These effects have been studied before [24, 15] so will not be consideredhere except for their effect on the image phase.With these assumptions in mind the formulation of motion compensation in the ideal casewill be investigated. This necessitates the development of the echo response from a pointtarget which depends upon the range history between the antenna and target, called the rangeequation. Following this, non-ideal compensation will be investigated considering such issuesas the effect of unknown terrain.3.2 Ideal Motion CompensationThe ideal case of motion compensation assumes that all the necessary geometrical information is available for compensation and that there are no processing restrictions. The necessarycompensation can be analyzed by considering how the echo from a point target (or impulsereflectivity) in the scene is affected by non-linear flight motion. The required compensationwill be described followed by an analysis of an unavoidable side-effect of the phase correctionoperation.3.2.1 Point Target FormulationThe point target formulation considers motion compensation to correctly image a singlepoint scatterer in a perfectly absorbing background. The radar transmits a linear FM chirpand the receiver compresses the echo using a matched filter. The time delay of the rangecompressed sine-like point spread function (PSF) depends upon the two-way distance the pulsetravelled to the point target [3]. This is called the range equation.The Range EquationThe ideal range equation is obtained for a straight flight and a flat earth (Figure 3):R(t) = + v2t (3.1)15where R0 is the slant range of closest approach, v is the aircraft ground velocity, and t is time.The phase of the range compressed point spread function for each pulse is [3]:= 4rR(t) (3.2)where ) is the carrier wavelength and the phase component due to delays in the radar systemhave been removed. The target energy in azimuth can be compressed by using this phasehistory to form a matched filter. The slant range R can be expanded in a binomial series to4th order in t:v2t v4tR(t) R0 + — (3.3)Therefore, the instantaneous Doppler frequency as a function of time is:1 ô 2vt v4t3 (3.4)It is clear that the position of the range compressed PSF will vary along the aperture, calledrange migration (RCM), according to the range equation. Range cell migration correction(RCMC) must be performed in order to extract the correct energy. A phase correction is notapplied to compensate for this because it is this Doppler phase history that is required in orderto compress the energy in azimuth. For instance, the quadratic term in the range equationproduces the linear FM part of the Doppler history. When the flight is non-linear the rangeequation becomes more complicated.16delv(t)Non-linear Flight MotionWhen the flight is non-linear the expression for R has a component that represents thedisplacement of the antenna from the nominal linear track along the direction of the antennato the target. It will be assumed that the along track velocity is constant:R1(t) JR10(t) + v2t (3.5)whereR10(t) = fn(O, H, R0, delh(t),delv(t)) (3.6)is the projected distance in the plane of the particular pulse and is a function of many parameters.When processing a single target only, there are two options:1. Form a matched filter using R1(t) and process. This involves resampling the range datato follow the target energy as it migrates due to R1(t). The matched filter will not bea linear FM chirp.VtRno(t)HFigure 3 Single Target Motion Compensation Geometry172, Modify the received data to appear as though the flight were linear, then compress usingnormal matched filtering and normal RCMC.In the second approach, the following calculation is required:LRe(t) = R1(t) — R(t)v2t 1’ v2t’\ (3.7)R10(t)+ 2R10tj 2RThe following phase correction can then be applied:= 4TrLRe(t) (3.8)as well as a resampling operation to extract the correct target energy along the aperture [15].The resampling operation is referred to as the motion compensation resampling operation.Typically, the range compressed data is at baseband and processing after motion compensation is based on this assumption. It will now be shown how the required phase correction,which varies across range, can invalidate this assumption even when the phase history alongazimuth has been adequately corrected [21].Range Varying Phase Correction (RVPC)An approximation can be made to the phase correction applied at the zero Doppler pulse(Equation 3.8) by defining (Figure 4):d(t) = delv(t)2 + delh(t)2 (3.9)a(t) = atari2(delh(t), delv(t))By assuming parallel rays this yields:R110(t)— R0 d(t) cos (a(t) + 0) (3.10)giving a phase correction of [21]:4irLq(t) —,-d(t)cos(a(t)+0) (3.11)18AntennaFor a given pulse the phase correction varies as a function of the off-nadir angle (0) to thetarget. Differentiating Equation 3.11 yields:A change of variables using:andRcosO = H(3.13)(3.14)RTRnIo(t)Figure 4 Zero Doppler Applied Motion Compensation Cross-sectionwhere1/H0=cos (—\ Baoand H is defined in Figure 3.(3.12)9047r—--d(t) sin (a(t) + 0)ãzq(t)8R80 18R Rtan0yields______4ir d(t) sin (a(t) + 0)This describes the range varying phase correction (RVPC) which has been noted before [211.RtanO (3.16)19This range varying phase correction, which is locally linear, will alter the range spectrawhich can invalidate the baseband assumption that is usually made by subsequent processingsteps. By using the following Fourier Transform property an estimate of the resulting rangespectral shift can be made:.F{g(t)}= (f)(3.17){g(t)e_i2t}= (f + f°)The magnitude of the frequency shift can be calculated from the range time varying phasecorrection, yielding:c 1 47rdjj03=r )RtanO (Hz) (3.18)If the baseband assumption is to be valid all subsequent processing must have transfer functionswith a constant spectral response as well as a sufficiently large range oversampling factor. Thelimit on the oversampling factor becomes:BWr Fr2 +If0I<-- (3.19)where BWr is the bandwidth of the range chirp and Fr is the range sampling rate. Thissimplifies to:1____DTTT rJ__) V Vor1 JJos 321\RtanOBW < rwhere O = - is the oversampling factor. For typical CCRS parameters (Table 1):1 + O.0321d±1 3 < O. (3.22)This provides a limit to the allowed displacement from the reference track. Given that thephase modulation varies non-linearly across range, significant distortions to the range spectracan occur even if the oversampling factor is sufficient.20Table 1 Typical CCRS ParametersParameter Typical ValueWavelength ). 56.56 mmAltitude H 6 kmBaseline length b 2.8 mBaseline angle c 40 degreesSlant range Rao 10 kmPulse Repetition Frequency FRF 337 HzProcessed Aperture T 3 sec 2.23 degreesRange Bandwidth BW 25 MHzRange Sampling Rate F 37.5 MHzVelocity v 130 m/sOne example of subsequent processing that is sensitive to the spectral shift is RCMC,which is normally accomplished by baseband interpolation. Normal baseband interpolation (Iand Q channels interpolated separately with a baseband kernel) can be performed provided therange oversampling factor is large enough and the spectrum of the interpolation kernel if flatand real valued, When the shift becomes large non-baseband interpolation is required. Giventhat the shift is range dependent this implies performing interpolation with a range varyingkernel which is undesirable.Since the range varying phase correction depends upon the displacement from the referencetrack, the magnitude of the range spectral shift can be minimized by segmenting the referencetrack to follow the aircraft drifts. The consequence of this is that the data becomes discontinuousalong azimuth across the boundaries. Alternately, one might expect that a real-time range gatedelay adjustment can be used (if available in the hardware) to maintain a constant swath centerposition. This will help the alignment of the target energy along azimuth but does not help therange varying phase correction. The range gate delay adjustment does not change the phaseof the echo energy from a given ground scatterer.21It is important to note that because the displacement from the reference track varies alongtrack the motion compensation required at each pixel of the raw data depends upon the locationof each target on the ground. This means that in order to perform ideal compensation theentire aperture for each target on the ground must be processed uniquely and separately. Thisleads to prohibitively large computational loads. For this reason non-ideal approximate motioncompensation approaches are usually used, The errors resulting from this will now be analyzed.3.3 Non-Ideal Motion CompensationTypical of most engineering problems is the trade-off between the quality of the resultand the resources and time required to generate the result. Given an unlimited amount oftime and computer resources and complete knowledge of the geometry of the situation, idealmotion compensation could be used. In reality, much more efficient and therefore approximatesolutions are required. This leads to the zero Doppler pulse processing approach that handlesmultiple targets at once.3.3J Zero Doppler Pulse ProcessingA second point target will now be added to the previous single target ideal motioncompensation approach to examine how multiple targets can be processed simultaneously.Multiple Point TargetsBy considering Figure 5 it is clear that there are different phase corrections to be appliedfor a given range cell of a given pulse depending upon which target you are correcting for, evenif the targets are at the same elevation. As the echoes from each target are overlapped alongthe aperture, only one correction can be applied for the pulse at this range bin if the targets areprocessed together. The only way to process each target correctly is to process the apertureseparately for each target. But, this is inefficient and is often not necessary. Therefore, one canuse R10(t) for each target [251. In other words, the center of aperture (or zero Doppler) pulse22is corrected properly for each target but the correction has a small error across its aperture foreach target. This can be called zero Doppler pulse processing.PulseFigure 5 Multiple Target Motion Compensation GeometryThe slant range correction for zero Doppler pulse processing is:yielding a phase correction of:zS.R(t) = R10(t)— R0 (3.23)zq5(t) = 4ir/.iR(t) (3.24)In order to extract the correct target energy the range data must be interpolated to obtainR10(t) instead of the normal sampling range of R0.An analysis of the error in the applied motion compensation as the pulse position variesfrom the zero Doppler pulse can be made. The difference between the desired phase correctionfor target 1 at a given pulse Lqe(t) and the phase correction actually applied at the slant rangedelv(t)VtRnlo(t)HTarget 2/et123where the echo energy resides q(t) is:47r / v2t v2tLS4e(t) — t)— —A \nlot) 2 (3 25)2rvt (R0 — R10(t)A RBy incorporating the parallel ray approximation one gets:A2irv2td(t) cos (a(t) + 0(t))Lk/err ) ——_____________________________(3.26)2irvtd103(t)A Rwhere d108 is the projected line-of-sight displacement. Another way to view this effect is toconsider this error as a mismatch between the FM rates of the corrected data and the matchedfilter. This approach is described below.FM Rate ErrorsMotion compensation has attempted to modify the data to appear as thought it was generatedfrom an ideal trajectory. One important case to consider is when the flight is parallel but offsetfrom the reference track and the terrain is flat (Figure 6). For this case there is a constant phasecorrection applied along the aperture for each azimuth line. The consequence of this is thatthe FM rate of the data is unchanged. If the motion compensation had been able to considerthe RCM along the aperture for a single target (ideal compensation) then the FM rate wouldhave been changed. It is therefore important to use an FM rate for the matched filter thatmatches the data, not the reference track, If there is a difference between the effective slantrange of the data and that of the matched filter it is easy to calculate the phase error resultingfrom the different FM rates. The FM rate is:=(3.27)so, the phase error is therefore:iqerr == 2t(,,,— “) (3.28)24HFigure 6 Constant Offset Flight MotionNotice that this error is of the same form as the multiple target scenario. The effect of thisinherent problem in non-ideal motion compensation can be reduced by adjusting the matchedfilter to match the data from the true antenna position rather than from the reference track.Given that this adaptive matched filter will vary along the aperture as the antenna deviatesfrom the reference track, it will be necessary to define a nominal filter if compression isimplemented using fast convolution [26]. The residual phase error will be due to an effectivedifference in the slant range used to generate the matched filter /Rmf.In Appendix A the effect of quadratic phase errors on the phase and shape of the compressedpeak has been analyzed. It was found that for quadratic phase errors of less than about ir atthe edge of the processed band, the peak phase error was about one third of what was applied.This is in agreement with results from [271. From simulation it was found that if the appliedquadratic phase error is less than about ir/2 then the broadening is less than about 5% and thesidelobe level increases about 4 dB. This is in agreement with results from [24]. Ignoring thepeak phase errors for the moment, a limit on the allowable quadratic phase errors can be made to25keep the peak deformation reasonable, For instance, one could (somewhat arbitrarily) say thatthe phase error must not vary more than - over the aperture. Therefore, from Equation 3.28:21rVLRmf (T)2 (3,29)and\R2/Rmf< v2T(m) (3.30)where LRmf is the displacement of the true antenna position from the one used to generate thematched filter and T is the duration of the processed aperture. For typical CCRS parameters(Table 1) this corresponds to a maximum displacement of 37 m between the antenna and thereference track. This is not normally a problem for the single channel case,Effects of Range Cell Migration Correction (RCMC)Another consequence of zero Doppler pulse processing is the coupling effect between therange varying phase correction (RVPC) and RCMC [21]. Consider again the case of a trajectorywhich is offset but parallel to the reference track (Figure 6). RCMC extracts the target energyfollowing the range migration. Due to the zero Doppler pulse processing and the RVPC theeffective phase correction that is applied for a given target varies along the aperture. Theconstant component is the zero Doppler phase correction and the varying component is equalto the product of the slope of the RVPC and the range cell migration (RCM):app (t) — R20) + t)1RRcM(t)-d(t) cos (c(t) + 0)— (t) + ) (-) (3.31)47r 2r d(t) sin (c(t) + O)v2t—-d(t) cos ((t) + 0) — R2 tan 0There are two possible ideal cases to compare this result with. The first is if the data isexpected to be perfectly corrected to the reference track such that the reference track is used26to define the FM rate of the matched filter, In this case, the ideal applied phase correctionwould be:I4ideal1(t) = —(Ri(t) — R2(t))47r7 v2t” 1 1—— (,(R10—R2o)+—(3.32)2 R10 R2047r 2r d(t) cos (a(t) + O)v2t-.--d(t) cos (a(t) + 0) — R2Therefore, the resulting quadratic phase error along the aperture, by usingsin (a— b) sin a cos b — cos a sin b (3.33)is:err12vtd(t)(((t) + 0) — sin (o(t) + 0))R tan (3 34)2Trvtd(t) sin (t)sin0The alternative is to assume that effectively a constant phase correction has been appliedalong the aperture that has not altered the FM rate. In this case, the filter is matched to theoriginal position of the antennas. The ideal correction for this approach is:zq5ideal2(t) = R10 — R20) (3.35)Therefore, the resulting quadratic phase error along the aperture is:_____2vtd(t) (sin (a(t) +0)err2‘ 9t/.R(t— tan i (3.36)This quadratic phase error is slightly smaller, on average, than the previous one, therefore, theantenna positions should be used to define the filters.If a peak error is estimated using the perpendicular to the line-of-sight displacementcomponent, one obtains:A__________2 (.AR tan027Furthermore, if a limit of ir/2 is placed on the quadratic phase error, the maximum alloweddisplacement from the reference track is:,\R2tanOd103< 4v2t(m) (3.38)For typical CCRS parameters (Table 1) this yields:djj08 <49 m (3.39)The easiest way to reduce this effect is to segment the reference track to follow the antennaalong track so that the displacement is small enough to produce insignificant range varyingspectral shifts. But, this leads to discontinuous phase data.3.3.2 Effects of Unknown TerrainAnother important source of error is the unknown terrain elevation, When the elevationof the scene patch (or point target) is unknown it is impossible to correctly calculate the slantrange shift that occurred [15]. An approximation is to resample to the R,10(t) points and applythe phase correction assuming flat terrain at the reference level (Figure 7). It should be notedthat the resampling will only extract the correct part of the echo when the target is actuallyat the reference level.28AThe applied phase correction assuming flat terrain at the reference level is:qmc(applied) =— Rao) (3.40)for the zero Doppler pulse of a particular target. The proper correction is:cbmc(ideal) —f(R4h— Baa) (3.41)The error in the applied phase correction is therefore:err = mc(app1ied) — mc(ideal) = —(R0— Rjh) (3.42)When the slant range distances are expanded in a binomial series and parallel rays are assumed,one obtains (Appendix B):err(t)4hd103(t) (343))Rao sin 0(t)29dRTA: antenna positionRT: reference trackRaoReference LevelFigure 7 Motion Compensation to a Reference Trackwhere 8(t) is the assumed angle to the target which will vary according to the geometry, andd±10(t) = d(t) sin (a(t) + 0(t) + (t)) (3.44)is the component of the displacement that is perpendicular to the line-of-sight direction, andz0(t) 0h — Oao is the error in the assumed angle to the target.It is clear that as the flight deviation varies with time (d(t) and a(t) change) and so willthe error in the applied phase correction. If the displacement has a constant offset term with asmaller time varying term the constant term will produce a phase bias along the aperture whichwill be unaltered by the compression operation.The form of the phase error is very nearly the same as the form of the flight motion. Forinstance, a linear deviation from the reference track along the aperture leads to approximatelya linear phase error. The same is true for quadratic and higher order motions. This can beseen by noting that the sine term will be approximately linear in a(t) provided the variationis small, When the angular variation is large d(t) must be small (close to the reference track)resulting in a very small non-linear component. One can therefore proceed to analyze the effectof linear and higher order cross-track flight motions by examining the effects on compressionof linear and higher order phase errors along the aperture. This has been described in detailin Appendix A and C.Cross-track VelocityAppendix C examines the effect of linear phase errors on the compression of a linear FMpulse. It was found that the peak is shifted in position, identical in shape, and has a phaseerror at the compressed peak. More specifically, a linear phase error = 7rL produces a shiftin the compressed peak position of:=—& = 4rv2ôt (sec.) (3.45)30where K is the FM rate. The corresponding phase error at the peak is:— L2— R0 (N2 346err— 4K 8v2 at)This shift and phase error can be put in terms of changes to the motion compensationphase errors along the aperture:4irh ad108(t) (3 47)at at .\RsinO OtThe derivatives involving the off-nadir angle are ignored as they do not contribute appreciably.The unknown parameter required to predict the shift is the target elevation (or equivalently thecorrect off-nadir angle). If an estimate of this were available one could predict the shift andcorrect for it. The same is true of the phase error at the peak. The resulting shift is:h ad8—. (sec.) (3.48)v2srnO atIf this is limited to less than one sample then:2‘-‘Llos V SIll 3 49hFRFwhere FRF is the Pulse Repetition Frequency. For typical CCRS parameters (Table 1) thisyields:d10 40 m (3 50)at <hsWhen the phase error at the peak is evaluated one obtains:2Trh (ad10524peak ARv2 sin2 at ) (3.51)The limits on the phase error of the compressed data depends upon the use for the data. Forinterferometry this effect is negligible.31Cross-track AccelerationIt is difficult to predict analytically the effect of non-linear phase errors on the compressedpeak PSF. In Appendix A an approximate analysis is given for the effect of quadratic phaseerrors on the phase and shape of the compressed peak. This corresponds to uniform translationalacceleration coupled with unknown terrain. Ignoring the peak phase errors for the momentone can place a limit on the allowable quadratic phase errors to keep the peak deformationreasonable (to be less than . over the aperture). It is desirable to extend this limit to representmaximum flight motions given a particular imaging situation.It is shown in Appendix B that the largest variation in the phase error occurs for motionperpendicular to the line-of-sight of the radar. This is interesting as the phase of the receivedsignal is sensitive to displacements in the line-of-sight direction. But, for terrain induced phaseerrors of concern is the variation of this error as a function of flight motion. The approximatemaximum phase error from Equation 3,43 is:4K h zd08352Yerrmax Rao sinOwhere dJjø8max is the maximum flight motion displacement perpendicular to the line-of-sight. If one limits this motion to be such that the motion compensation phase variation overthe aperture is less than from uniform acceleration, this yields:/ \2-ajj08 = ..ajjosT2 (3.53)The resulting upper limit on uniform acceleration is:,\RaosiriO fm’ajj03 < 1 (3.54)niFor typical CCRS parameters (Table 1) this is:ajj05 < () (3.55)where g is the acceleration due to gravity.323.3.3 Effects of No Motion Compensation ResamplingUp to this point in the analysis the range gate delay adjustment of the echo energy, whichis the motion compensation resampling operation, has been mostly ignored. SAR processingis most sensitive to phase but still relies on the signal strength. Some simplifications to thephase processing result if the range gate delay adjustment is ignored because the slant rangesample spacing is constant. The effect of this is now considered.If resampling had not been done and the processing had simply sampled the range data atfixed slant range values and phase compensated these samples then the effect of turbulent flightis to produce some uncompensated range cell migration (RCM). Uncompensated RCM leadsto broadening in both the range and azimuth directions. Figures 8 and 9 show flight profilesfor the two cases of resampling and no resampling, indicating from where on flat terrain thesamples would be taken, Figure 10 shows an example of the no-resampling case.The error in the sampled slant range for a particular pulse from Figure 10 is:/.R(t) R8 — R0 (3.56)One can approximate the difference in the slant range values by considering Figure 11. Thedifference is approximately the projected line-of-sight displacement from the center of apertureposition (An):— R0 x(t) sin(t) i4,(t) (3.57)Given that typical flight motion over an aperture is small compared with the slant rangeresolution, this effect is not expected to produce significant broadening. This will be confirmedby simulation in Chapter 5.An alternative to motion compensation resampling is to have a real-time range gate delayadjustment which maintains a constant ground range at the swath center. This range gate33Horizontal Motion //V/VV V/ FFigure 8 Flight Profile Without Resamplingdelay adjustment alters the start time of the echo receive window. This will provide a roughcompensation for aircraft drift but only correctly compensates for one slant range (assumingthe terrain assumption is correct as well).3.3.4 Extension to Clutter DataSo far the analysis has considered the effect of motion compensation errors on point targets.Since typical imaged terrain has many contributing scatterers per resolution cell one must extendthe results to consider this. Given that the aircraft flight motions are essentially constant duringthe time corresponding to a single resolution cell (of order 30 ms) it is reasonable to expect thatall contributing scatterers for a given sample value are affected by the motion in an identicalVConstant R varying 934VFigure 9 Flight Profile With Resamplingway. For instance, all the point spread functions corresponding to the contributing scattererswill be broadened or shifted by the same amount because each scatterer will experience thesame phase error along the aperture. This means that all scatterers will be altered by the sameamount leading to an altered clutter value.3.3.5 Effects of Inertial Data ErrorsIn the preceding analysis the position of the antennas was assumed to be correct. In reality,there will be errors in the inertial data. Inertial data errors basically fall into two categories:high frequency errors and low frequency drifts. The high frequency errors have periods thatMotion Compensation With Resampling////Horizontal MotionV.V/V///35detv(t)RTFigure 11 Cross-section of Neglecting Resampling Geometryare short compared to the processed aperture. The low frequency drifts have periods that arelong compared with the processed aperture.VtReHFigure 10 Neglecting Motion Compensation Resamplingaperture center position AAcdelv(O)ReRts36The low frequency drifts will yield effectively constant antenna position errors during anaperture synthesis. These errors will behave much like the errors due to the unknown terrainand will depend upon the motion during the aperture. For instance, a roll bias will behave verysimilarly to an error in the assumed target elevation because it is basically a rotation error in thegeometry. Most biases, such as a roll bias, are small and should not have a significant impactat the motion compensation stage. Some work has been done to determine upper bounds oninertial data errors for a specific system by simulation [17]. In addition, a more general rulehas been proposed that for a maximum broadening of 10% the phase errors must be less thanK!2 when the period is greater than the processed aperture [24]. In terms of INS errors thiscorresponds to about 7 mm for the CCRS C-band system.The high frequency errors translate into small high frequency errors in the antenna positionestimates, which in turn become high frequency phase correction errors along the aperture.Limits on the allowable high frequency phase errors in terms of the impulse response of thesystem have been analyzed [24, 28]. To maintain an Integrated Side Lobe Ratio (ISLR) of atleast —20 dB, one standard deviation of the high frequency phase error must be less than about 6degrees. In terms of INS errors this corresponds to about 0.5 mm for the CCRS C-band system.A preliminary investigation into the effect on the resulting phase of the impulse response hasbeen made [29]. This will be extended further in the following section.Directional Random Walk for Point TargetsThe effect of high frequency (or random) phase errors on the phase of the azimuthcompression result can be modelled. Azimuth compression can be thought of as the repeatedsummation of complex samples. When there are errors in the phases this will lead to peakdeformation and phase distortion. To simplify the problem, consider only the resultant of thesummation that yields the complex peak value. Furthermore, assume that the expected phase37of the compressed peak is zero. Each sample that contributes to this summation is the resultof the complex multiplication of the data sample with an estimate of the complex conjugateof that sample. The matched filter attempts to remove any imaginary component in the datasample leaving a purely real result. When there is no noise and the filter is perfectly matchedthen the peak value will be purely real valued because every vector summed will be directedalong the positive real axis in the complex plane. When noise is present or the filter is notcorrectly matched in phase, the contributing sample will have an imaginary component. Thecomplex peak value can be modelled as the result of what could be called a directional randomwalk. Each sample in the summation can be thought of as taking a step in the complex planein the general direction of the positive real axis but with a small random direction error.The statistics of the complex summation result for the case of communication theoryproblems has been studied [30]. The summation can be characterized by:Re = X + jY=(3.58)where each contribution Ae’ is independent and the expectation of Y is zero (expect zerophase). The A’s can be regarded as the modulated scattering strength along azimuth includingany compression windowing. The calculation of phase variance at the peak is carried out inAppendix D for Gaussian phase noise. The variance to be solved is:2 <(x—)> <y2>- =- (3.59)X2The resulting phase variance using integrals from [31] is:=(1 — e2) A(3.60)2e° (N 2L i=138where cr is the variance of the Gaussian phase noise and the A ‘s represent the compressionweighting. When there is uniform weighting this reduces to:/ 2(1 —=2 (3.61)2N ed’The validity of this model to describe the phase smoothing ability of compression of a pointtarget was verified using simulation with typical parameters for SAR. The simulation resultsagreed to very high precision with the theory. The intuitive explanation for the significantsmoothing observed (‘-.‘ -) is that the compression operation has little effect on the distributionof the phase noise within the data but collapses most of the signal energy into a small region.The phase noise at the compressed peak is very small due to the significant increase in thesignal to noise ratio at that point.INS Errors for C’utter DataWhen the resulting phase errors for the clutter data case is considered very different resultsare obtained. Previous work [29] and repeated simulations and real data processing in thecurrent work have demonstrated that random phase errors are not smoothed when the signal isclutter data. The real data experiments added random phase noise to uncompressed SAR dataand compared the focussed results with the noiseless compressed image. There was no evidenceof smoothing. The resulting phase errors depended upon the specific processing conditions,such as the processed beamwidth. The intuitive understanding of this is that the clutter to noiseratio at a particular pixel is unchanged by compression. The individual point target responsesare focussed but the overall signal strength at any one point is, on average, unchanged. Thismeans that high frequency phase errors are transferred to the compressed image result withoutthe effect of smoothing.3933.6 Summary of Single Channel Errors and ConstraintsIn ideal motion compensation all the required geometric parameters are known and thereis no restrictions placed on the processing. This leads to accurate compensation except thatthe range varying phase correction induces a range spectral shift. This can be minimized bykeeping the reference track close to the antenna position by segmenting the track.In realistic processing situations non-ideal motion compensation must be used. In thisscenario multiple targets must be processed together and certain important information, such asthe terrain elevation, is not known. This leads to various errors (Table 2). Given the assumptionsmade at the beginning of the chapter the following tables summarize the phase errors occurringin single channel non-ideal motion compensation and the limitations these errors place on theimaging configuration and flight motion. The zero Doppler pulse processing leads to phaseerrors along the uncompressed azimuth data due to multiple target processing, possible FMrate errors, and the coupling between the range varying phase correction and RCMC. Giventhat the target elevation is not typically known, assuming flat terrain also leads to phase errorsalong the aperture. Inertial data errors translate into phase correction errors along the apertureproportional to the error in the antenna position. This leads to peak deformation and phaseerrors.Table 3 indicates how these errors translate into limitations on the allowed flight motionfor a given imaging configuration. Specific numbers are included for the nominal CCRS case.• Errors in the slant range value used to generate the matched filters will cause peakbroadening.• The coupling between the range varying phase correction (RVPC) and RCMC producespeak broadening that depends upon the displacement from the reference track perpendicularto the line-of-sight direction.40Table 2 Single Channel Phase ErrorsPhysical Effect Phase Error Along ApertureMultiple Target Consideration2irvt2d03 (t)cerr(t)FM Rate Errors27rvt/Rmf/-‘4err(t) AR2Range Cell Migration Correction2irvtdj103(t)Lcerr(t) AR2 tanUnknown Terrain Elevation4ir h dj0(t)err() ARsinOInertial Data ErrorsLqerr(t)• The range varying phase correction (RVPC) causes a non-linear range spectral shift whichdepends upon the displacement from the reference track perpendicular to the line-of-sightdirection. This can be an important factor when subsequent processing is considered.• The effect of unknown terrain is to produce peak broadening for acceleration perpendicularto the line-of-sight direction and a peak shift for velocity perpendicular to the line-of-sightdirection.• The effect of inertial data errors is to produce peak broadening and phase errors. For highfrequency errors there is significant smoothing for point targets but not for clutter data.In addition, neglecting to resample for motion compensation will lead to some PSF deformation.41Table 3 Single Channel Motion Compensation Limitations (CCRS Typical)Physical Effect LimitationFM Rate Errors 5% broadening:\fl2/Rmf< vT<37mRange Cell Migration Correction 5% broadening:,\R2tand±103 <v2T<49mDisplacement Offset Range Oversampling Factor Or:(range varying phase correction)2cdjj08\RtanOBWr <1 + 1.032 k&uosI < OrUnknown Terrain Elevation 5% broadening:AR sin 0ajj05 < hT25g m<——h s2one sample shift:öd±108 v2 sin0öt hPRF40 m<——hsInertial Data Errors High Frequency INS ErrorsClutter Data: no smoothingPoint Target:/ 2I 1 — e°2 \°9 22N42Chapter 4 Motion Compensation for InterferometryNow that the foundations of motion compensation for a single channel SAR have beenformed, it is now possible to consider the specific problem of motion compensation for InSAR.This chapter will begin by defining the assumptions made and the requirements of motioncompensation for InSAR. Following this, ideal motion compensation will be investigated andthe corresponding interpretation of the differential phase of the interferogram. The requiredinterpretation for non-ideal motion compensation will then be analyzed for the case of singleand dual reference tracks. Finally, InSAR errors resulting from non-ideal compensation willbe investigated by extending the preceding single channel error analysis to the two channelcase under the same assumptions.4.1 AssumptionsThe assumptions made are the same as those presented in the previous chapter, namely:1. Both antennas are yaw steered to zero Doppler. This means that both antennas are assumedto be pointing perpendicular to the velocity vector of the aircraft.2. The SAR processing is to zero Doppler for both channels.3. The pitch of the aircraft is fixed.4. The spacing between transmitted pulses on the ground is kept fixed by adjusting the PulseRepetition Frequency (PRF) to be proportional to the aircraft velocity.5. The inertial data will be assumed to be correct, for the most part.6. The motion compensation to be applied will be post-flight, apart from those mentionedin the above list.7. Translational motion and aircraft roll must be compensated through motion compensation.43These basic assumptions are consistent with the CCRS Convair 580 InSAR system [18], and theEnvironmental Research Institute of Michigan (ERIM) IFSARE system in strip-map mode [32].4.2 InSAR Motion Compensation RequirementsSAR Interferometry is concerned with extracting the phase difference between two processed images, pixel by pixel, in order to estimate the terrain elevation, InSAR is no longerconcerned with the “focus” of the image in terms of its radiometric appearance, but rather howthe focus affects the estimate of the phase, and in turn the estimate of elevation. This requiresviewing the motion compensation and azimuth compression stages of conventional airborneSAR in a different way.For example, InSAR is not affected by large phase bias errors that exist in each channelprovided they are common to both and thus cancel out in the differential phase. As the twoantennas are rigidly connected on the aircraft, much of the motion experienced by each channelwill be common to both. Therefore, phase errors resulting from the common flight motion maynot degrade the interferometer. What affects the height estimate is differential phase biasesand differential phase variances, which translate into height estimate biases and height estimatevariances.In a recent paper by Bamler and Just [27], a frequency domain analysis was conductedto evaluate the effect on the differential phase of differences in the transfer functions between channels. For distributed homogenous targets the signals for the two channels can becharacterized by:S1 = (4.1)82 = R2e38 (4.2)44where each is a circular Gaussian random variable [101. The expected differential phase biaswas shown to be the phase of the inter-channel correlation coefficient [27]:bias = arg (‘y) (4.3)* f HA(f)H(f)df<s1s2>—007= =______________(4.4)vSi9’><S2> 00HA(f)2IHB(dwhere < > is the expectation operator, and HA(f) is the transfer function of channel A, andHB(f) is the transfer function of channel B. The probability density function of the differentialphase can be shown to be (page 46 of [331):() = i _2) (45)where/3= 7I°() (4.6)The corresponding variance has no closed form solution, but can be evaluated numerically:4 =f p()d (4.7)This work formalizes the idea that differences in the transfer functions lead to differential phasebiases or increased differential phase variances. For instance, if the motion compensation errorslead to differences in the PSF’ s between channels this will lead to a decrease in the inter-channelcorrelation and therefore an increased differential phase variance. It is shown in Appendix Ethat the inter-channel correlation coefficient factor due to inter-channel azimuth broadening ofj > 1 is:7CC— (4.8)for sinc-like point spread functions. Therefore, inter-channel broadening of about 10 percentleads to a correlation factor of about 0.95.45But, this does not mean that defocussing common to both channels has no effect on theInSAR system. The amount of phase noise from certain sources, such as baseline decorrelation[10], depends upon the transfer function in each channel. In addition, a more local estimateof the terrain elevation is achieved with a narrower point spread function (PSF). Given thatthe height estimate is very sensitive to the differential phase, significant spatial averaging isrequired in InSAR to facilitate phase unwrapping and reduce the variance of the estimate. Thisrelaxes the localization requirements because the averaging will typically take place over manyresolution lengths. But, the narrower the PSF the more uncorrelated the samples are and thusbetter phase smoothing will be accomplished. The amount of spatial averaging may need to bereduced when the local elevation variation is large in order to avoid differential phase aliasing,This will increase the sensitivity to focussing problems.Once the differential phase has been calculated it must be unwrapped into the absolute phaseby adding some multiple of 2w. This is a two dimensional problem which requires continuity inthe differential phase [191. If the differential phase is discontinuous, for instance by referencetrack segmentation, then the unwrapping process becomes very complicated. Discontinuitieswill be inevitable in rugged terrain due to shadow and layover effects but should not be causedby motion compensation mechanisms.In summary, identical, high resolution, geometrically correct, PSFs are desired for InSARmapping. Phase errors in the compressed images that are common to both channels cancelout in the interferogram’s differential phase. InSAR motion compensation has some of theconventional requirements, such as resolution and sidelobe level control, that exists in thesingle channel case. But, the primary concern is to maintain the similarity in the transferfunctions between the two channels. Differential phase errors should be minimized at all costs.In addition, the differential phase should be continuous for phase unwrapping. The ability tomeet these requirements of motion compensation for InSAR will now be examined for the46single and dual reference track approaches.43 Ideal Motion CompensationIdeal motion compensation for InSAR is very similar to the single channel scenario. Twoparallel reference tracks can be defined, one for each channel. Given precise knowledge ofthe entire imaging geometry, including inertial data and terrain information, each target onthe ground is processed independently as prescribed in the ideal single channel case. Thecompensated range compressed data then emulates transmission and reception from the twoparallel tracks. The transmit/receive channel is processed as in the single channel case, but,the receive-only channel calculates corrections based on the transmission from the two-wayantenna and reception at the one-way antenna.As was the case in the single channel scenario there is a range varying phase correctionthat depends upon the displacement of each channel from its reference track. The possible side-effects from subsequent processing were outlined in the previous chapter, but, for InSAR theconcern is what will lead to differential phase errors. If the displacements from the referencetracks are almost identical in each channel one would expect that the phase errors wouldessentially cancel out. This will be investigated further in Section 4.4.5.Another approach could be to motion compensate both channels to the same referencetrack. If ideal compensation was carried out then the resulting differential phase would be zeroeverywhere because the phase would be the same in both channels (ignoring other effects suchas noise). The usefulness of such an approach seems remote at present but it highlights themain differences between the two main approaches to motion compensation for InSAR, namely,using single and dual reference tracks. The single reference track approach attempts to undothe relative mapping between the two channels. The dual reference track, on the other hand,47attempts to make the parallax be with respect to the two reference tracks. The two approachesrequire different interpretations of the resulting differential phase.4.3.1 Dual Reference Track Assuming Ideal CompensationThe approach taken by the Jet Propulsion Laboratory (JPL) [21, 22] is to assume that idealmotion compensation has been performed and they proceed to interpret the differential phasefrom the two reference tracks (Figure 12). The interpretation of the differential phase used isbased purely on triangulation from the reference tracks (RTa, RTb):= (Rtah — (4.9)The target coordinates are measured from the channel B reference track, requiring 0b and Rtbh.Quoting their results [21, 22]:= cos1 (cog /i1 sin (8,—— sin sin (Ob—(4.10)whereID ).DP’2 122 z2kibh— E)— ‘tbh — VRTSlfliub—c)hURT “tbhand=—cxab (4.12)This approach does not consider the possible biases that result when motion compensationis performed assuming flat terrain at some reference level. This assumption is required forInSAR because the terrain is not known apriori.If there are large drifts from the reference tracks then segmentation of the reference trackswill be required due to the range varying phase correction (RVPC) problems. This will resultin discontinuous differential phase because the interpretation of the differential phase is fromthe two reference track positions. This will complicate the phase unwrapping process.48Bc4.4 Non-Ideal Motion CompensationGiven that InSAR would be superfluous if the terrain elevation was initially known idealmotion compensation is not possible for InSAR. At best, each individual target could beprocessed as in the single channel ideal case but with flat terrain assumed. As in the singlechannel case, this is very computationally intensive and therefore zero Doppler pulse processingis most often used. In this section, two new approaches to interpreting the differential phaseby considering the biases caused by unknown terrain will be presented. Following this, theFigure 12 Double Reference Track Geometry49effects on the InSAR height estimate, due to the same sources of errors treated in the singlechannel case, will be investigated.4.4.1 Single Reference Track with Unknown TerrainThe single reference track approach to motion compensation was first used for InS AR bythe Canada Centre for Remote Sensing (CCRS) [25, 18]. Both channels are phase compensatedto the same reference track using the geometry of Figure 13, assuming the terrain is flat atsome reference level. For the following analysis the arc on which the target lies is assumedto be centered at the channel A antenna. An equally valid view centered on the channel Bantenna could be used. The fact that there is a choice highlights the fact that the two channelsare not perfectly registered due to the unknown terrain. When clutter data is processed the best3D position estimate for the pixel is essentially midway between the two results.BcFigure 13 Single Reference Track GeometryAc: channel A aperture center positionBc: channel B aperture center positionbYaReference Track RbhReference Level50The applied phase correction for each channel assuming channel A is transmit/receive andchannel B is receive-only is (Figure 13):cbA(applied) = —(R0— Rao)2ir (4.13)qB(applzed)= _T((Rt0 — Rao) + (R0 — Rbo))If ideal compensation had been applied (by virtue of 0 being known):qA(ideal) = —-(Rh— Rah)(4.14)2ircbB(zdeal) =—T((Rth — Rah) + (Rth — Rbh))the differential phase would be zero everywhere because the relative phase mapping due to theparallax would have been removed,Due to the realities of non-ideal compensation it is necessary to make a “center of apertureassumption”. The assumption is that motion compensation and the resulting differential phaseof a compressed image sample of the interferogram is due only to the relative geometry ofthe InSAR system at the zero Doppler pulse. If the applied phase correction for a givenpixel is “backed” out after compression, the resulting differential phase is due simply to theinstantaneous relative geometry of the antennas and target for that pulse [181. In addition, theinter-channel registration resampling that is performed on the range compressed data is assumedto carry over into the azimuth compressed data unchanged. The center of aperture assumptionis in fact correct for trajectories which are parallel but offset from the reference track(s). Thisis the case because the applied phase correction is constant for a given slant range cell alongthe aperture and azimuth compression therefore does not alter the constant phase correction.By invoking the center of aperture assumption the differential phase will be zero for scenepatches that have the same elevation as the reference level because the motion compensationfor these patches will have been done properly. The differential phase will be non-zero for anyscene patch whose elevation differs from the reference level. The value of this phase is thedifference in the errors in the phase correction applied between the two channels due to the51reference level assumption. These errors depend upon the imaging geometry and can be relatedexplicitly to the elevation of the scene patch relative to the reference level. This interpretationapplies for slowly varying drifts of the aircraft from the reference track. Higher frequencymotions coupled with the unknown terrain will cause phase errors that deviate from the centerof aperture assumption. These errors will be considered in Section 4.4.3.The resulting single reference track differential phase 4) can easily be shown to be(Appendix F):4) = (qA(applied) — cA(zdeal))— (qB(applied) — gB(ideal))2ir (4.15)= _T(R_)The off-nadir angle 0ah’ measured from the instantaneous antenna position A, can be obtainedin terms of known and measured parameters from simple trigonometry:Rh b2 + R0 + 2bRao COS (cxab + 0ah)(4.16)R0 = b2 + R0 + 2bRao COS (aab + Oao)yielding:4) R2a0 + 2bRao COS (a + 0ah) — Rbo) (4.17)and rearranging terms:0ah(( +R60)— b2 — Ro)— ab (4.18)The scene patch elevation above the reference level is (Figure 13):hYa.RCOSO (4.19)By expanding the slant range values in a binomial series the following approximate solutionresults:0ah C0S1 + COS (Oao + aab))—aab (4.20)52Using the channel B antenna as the arc center, by symmetry the following results:0bh + (0io + 0ab) — aab (4.21)\2-Irb jA more accurate position estimate can be obtained by using the slant range:R— ao+ bo (4222and the off-nadir angle:°ah + 0bh (4.23)measured from the center of the antenna baseline. The sensitivity of the height estimate toerrors in the differential phase can be shown to be:Oh \RsinO (424)2bSfl(O+ab)The single reference track differential phase (Equation 4.15) is independent of the referencetrack positions. Rather, it depends upon the actual slant range to the target from channel B(Rbh) relative to the slant range used for resampling channel B to register with channel A atthe reference level (R&o). If the reference track is segmented along azimuth to better followthe drifts of the aircraft, the differential phase will still be continuous provided the referencelevel is continuous across the segments. This is a very important feature of the single referencetrack approach.Another benefit of the single reference track approach comes from the fact that thedifferential phase wifi vary depending upon the variation of the terrain about the assumedreference level. In other words, the flat earth fringes have been removed. This implies thatthe range spectrum of the interferogram is approximately at baseband. This simplifies thephase smoothing operation because the flat earth fringe rate, inherent in the dual referencetrack approach, does not need to be considered. Typically, the phase smoothing would have53to estimate the local fringe rate and smooth about this slope [19]. But, in either case rangefiltering of the interferogram must be considered carefully so as not to break the requiredsymmetry of the spectrum [34].4.4.2 Dual Reference Track with Unknown TerrainOne thing that is missing in the ideal motion compensation dual reference track interpretation of the differential phase is consideration of the phase errors in each channel that resultfrom motion compensation assuming flat terrain. When these systematic terrain induced phaseerrors are included (or equivalently, motion compensation is backed out after compression)a new double reference track approach is obtained. The approach proceeds as in the singlereference track case by considering the difference in the errors in the applied phase correctionbetween channels. The center of aperture assumption is used here in the same was as in thesingle track case. Therefore, this approach attempts to handle errors due to unknown terraincoupled with slowly varying drifts from the reference tracks.The applied and ideal phase correction in channel A is (Figure 12):gA(applied) = Rtao — Rao)4 (4.25)cA(ideal) = _“f(Rtah— Rah)The channel B (receive-only) applied and ideal phase correction is:cB(applied)— Rao) + (Rtb0 — Rbo))2 (4.26)qB(ideal) f((Rtah — Rah) + (Rtbh— Rbh))The resulting differential phase will be the ideal component (Equation 4.9) plus the differencein the errors (Appendix G):= Rtah—(cA(applied)— qA(Ideal))— (cB(applied)—q5B(ideal)) (4.27)= (Rtao— Rib0 + Rbo Rbh)54The off-nadir angle °ah can be explicitly related to the measured differential phase and otherknown parameters using Equation 4.16:0ah + =(((F + Rtao — Rtbc, + Rbo)2 — Rao“\ (4.28)2bRao )By expanding the slant range values in binomial series the following approximate solutionresults:0ah +i (4.29)cos1 ((i) ( + Rtao — Rtb0) + COS (Oao +Again, the symmetric formulation using the channel B antenna as the arc center leads to:°bh + aab1 ) (4.30)cos1 ((i) ( + Rtao — Rb0) + CoS(Oao + aab))The unbiased position estimates use the same Equations 4.22-4.23 as before. The conventionalapproach measures the off-nadir angle from the reference track but the single and new doubletrack approaches measure from the instantaneous antenna positions. This is consistent with thebacking out motion compensation viewpoint.The differential phase bias error that the new dual reference track approach recovers, thatis inherent to the conventional approach, is shown in Appendix G to be:E(A)-E(B)(27rhbl0N (4.31)— aab)For unsegmented reference tracks, rugged terrain, and typical aircraft motion, this can correspond to height estimation biases of several meters. If the reference tracks are segmented toremain close to the antennas and the antenna and reference track baseline angles are similar,the errors become quite small.The sensitivity of the height estimate to differential phase errors is identical to the singletrack case (Equation 4.24). But, the differential phase depends upon the reference track55position through Rtao and Rtb0. This means that reference track segmentation will resultin discontinuous differential phase across the boundaries which wifi complicate the phaseunwrapping process. One solution to this is to undo the originally applied phase correctionto each channel after compression. The resulting differential phase is then interpreted withrespect to the instantaneous antenna positions, removing the reference track dependence andmaking the differential phase continuous. Another more efficient approach that combines thesingle and dual track approaches to obtain the benefits of both wifi be described later in thischapter (Section 4.5.2).4.4.3 Effects of Unknown TerrainWhen the terrain is assumed to be flat at the reference level and there are higher frequencyflight motions then there will be phase errors that vary along the aperture for each channel.This was analyzed already for the single channel case (Section 3.3.2). For InSAR, of concern isthe difference in these phase errors. The error in each channel has the same form. Specifically,for channel A (from Appendix B):a(t) 8irda(t) (aa(t) + O(t) + t)) sin (t)) (4.32)where Oa(t) is the assumed angle to the target, and zO(t) is the error in that assumed angle.All of these parameters will change along the aperture for a given target due to the changingflight motion. It is shown in Appendix F that for the single reference track case the phase errorfor channel B is the same as channel A plus an additional term:b(t) LqSa(t) + Lqab(t)a(t) + sin (aab(t) + Oao(t) + 9(t)) sin (o(t)) 4.33where b is the baseline length between the two antennas (fixed) and aab(t) is the off-verticalangle of the baseline (aircraft roll angle). The portion of this relative term that is constantalong the aperture will appear as a bias in the differential phase. This is precisely what is56considered in the preceding interpretations of the differential phase that includes the unknownterrain effects, Further approximations, following Appendix F, yields:/ 27rhbJj08(t)/babjt) (4.34)AR sin 0whereI 0ao+h”\b08 bs1n jab + 2 J (4.35)The same analysis can be done for the dual reference track case (Appendix G), yielding:Lq5a(t) 87rda(t) sin (ca(t) + 0a(t) + z0;(t)) (LOa(t))b(t) 8lrdb(t) (t + Obo(t) + LOb(t)) (0(t)) (4.36)= a(t) + /.4ab(t)The difference in the errors is shown in Appendix G to be:ab(t) = (Rtao— Rah + Rtbh—Rtb0 + Rb0 — Rbh) (4.37)Expanding these slant ranges in binomial series yields (Appendix G):4ab(t)I27rbRT’\—A)(cos (aRT + Oao) — (RT + Oah))+() (cos (aab(t) + Oao(t)) — COS (aab(t) + Oah(t))) (4.38)/47rbRT’\. I L0”\ . (L’x0’A ,) Sifl cXRT + Oao +——} sm —/ z0(t)’\ . /L0(t)S1 + 0a(t) + 2 ) Sifl 2The first sinusoidal term involving the reference track is constant for a given target becausethe reference track is fixed over the aperture. The second sinusoidal term is identical to theresidual phase error in the single reference track case. This term varies along the aperturegoverned by aircraft roll and translation. Therefore, the interpretation of the differential phaseis different between the two approaches but the phase errors due to unknown terrain coupledwith higher frequency flight motion are the same (Equation 4.34).57Aircraft RollIt is useful to evaluate how the difference term (Equation 4.34) varies along the aperture,as this will indicate the source of differential phase errors. The effect of a changing roll angledominates over the effect of the changing angles to the target. Assuming the roll angle variationis very small and expanding in a Taylor’s series:ab(t) ab(O) + (a(t) — aab(0)) (4.39)Oaaband9ab 4Kb / /zONabCOS ab(t) + Oao + Sill (4,40)withb108 = bcos (nab +Oao + Oah) (4.41)yields:ab(t) ab(0) +2Khb103(ab(t)— ab(0)))RsinOA I (rl\ 2K iosaa&(t)sin 0where /c(t) is the variation in the roll angle about the center of aperture position. Ofinterested is the time varying term, as the bias will be removed in later processing and doesnot affect compression. This yields a varying differential phase error along the aperture of:21rhboZXca(t) (4.43))Rsin0Angular Velocity For linearly varying roll one might expect to observe a relative mis-positioning between the two channels due to the unknown terrain elevation. Using the resultsof linear phase errors on compression from Appendix C one can define:2Khblosab(t)= (4.44)).RsinO58Differentiating yields:27rhb103 azaab(t)irL (4.45))Rsin8 9tThe resulting differential phase at the peak is:— “rh2b (OLS.aabN2 (446)4K — 2v2,\RsinO\. at )For any reasonable roll this differential phase error is negligible.The resulting differential shift in samples is:L FRF b108h FRF azxaab(t)=— 2K— 2v sin é öt (samples) (4.47)The effect of this differential shift on the variance of the differential phase can be evaluatedby considering the inter-channel correlation coefficient [27]. When the differential shift is lessthan about 1/16 of a resolution width (or about one twelfth of a sample for the CCRS case) thecorrelation coefficient factor is greater than about 0.99. This will be used to limit the alloweddifferential shift. With 1/16 of a sample as a bound on the relative shift one obtains:ô/Q.ab(t) v2sinO (448)at 8b103hPRFNote that the line-of-sight component of the baseline is usually a small fraction of the baselinelength. For typical CCRS parameters (Table 1) and a worse case b108 = 1 m, yields:ôzaab(t) 290 (deg.’\ (4.49)at h \\sJAngular Acceleration If consideration is made for quadratic roll motion (uniform angularacceleration) the phase error is quadratic. When Equation 4.4 is evaluated for a differentialquadratic phase error one obtains the same result as linear FM compression with a quadraticphase error evaluated at the peak of the compressed result (Appendix A), except that it isnonnalized [27]. Therefore, the expected phase bias (which is the angle of the correlation59coefficient) is 1/3 of the value of the quadratic phase error at the band edge. Therefore, thepeak error along the aperture must be less than 1/3 of the maximum allowed differential phaseerror (Appendix A):peak max> (4.50)In terms of uniform angular acceleration:/-ab(t) = (4.51)The resulting differential quadratic phase error, as before, can be transformed into a differentialphase bias error in the compressed image:aroii’rbioshT2-bias (4.52)l2ARsinOWhen processing to the 3 dB width of the antenna beam pattern,ART = T3dB = — (4.53)Dvand conversion to a height estimate bias is made using Equation 4.24, one obtains:Lhbias A2Rbiosaroii 4 54h— 24bvDsin (0+ aab) ( . )For typical CCRS parameters (Table 1), a bias ratio of 1 m in 1000 m, and a worse caseb108 = 1 m, yields:aroll <0.07 —-- (4.55)SIn addition to the differential phase bias there is decorrelation given by the magnitude of thecorrelation coefficient. For quadratic differential phase errors along the aperture the magnitudeof the inter-channel correlation coefficient can be shown to be [27]:/s2(/) +C2(vW)HI = (4.56)60where the Fresnel integrals are defined as:3(x) = fsinu2du0 (457)C(x) = /fcosu2duand ‘I’ is the differential phase error at the band edge. For a typical allowed differential quadraticphase error of 75 mrad at the band edge for the bias result (yields about a 1 m elevation biaserror), this leads to an insignificant decorrelation factor of about 0.995 by numerical integration.Therefore, future decorrelation due to differential quadratic phase errors will be ignored as thebias component is the dominant effect.Sinusoidal Roll The dominant roll motion is often oscillatory. This leads to oscillatorydifferential phase errors along the aperture. It is not possible to analytically evaluate theeffect on the differential phase of such motion, even if it is approximated by a sinusoid. Anapproximate solution can be obtained by considering the instantaneous angular acceleration ofa sinusoid and extrapolating from the angular acceleration results.Let the roll motion be described by:/aab(t) = G (4.58)Troitwhere G and Troij are the amplitude and period of the roll motion, respectively. Thecorresponding instantaneous acceleration by differentiating twice is:f \2,‘27r \ . 2irta,.01j = —G ) Sifl (4.59)\-rollJ 1roUThe effective angular acceleration will depend upon the length of the aperture compared to theperiod of the roll, but will always be less than 1. By using the previous results from Equation4.52, and defining a constant f <= 1, which relates the instantaneous to effective angular61acceleration, the resulting differential bias error is:bias— •fG()sin (4.60)12.\R Sifl 0 Troit TroitTherefore, the bias error will follow the roll motion,The roll motion will also change the expected differential phase. This change can beestimated from differentiating Equation 4.43, yielding:21rblosh/.crab(t)ARsinO‘4612irb108h . 2irtGsrn—SIll 0 TroitWhen the difference between the change in the expected and error in the differential phase arecompared, the following results:-bias — 1exp2irbj03hGsin 2iri (i — f7r2T) (4.62)sin 9 Troii 6TroiiTherefore, the effect of the sinusoidal roll is to smooth out the measured differential phase.When the last bracket in this expression is zero, the roll does not effect the measured differentialphase. The resulting height estimate using the instantaneous antenna positions will thereforehave an error term that follows the sinusoidal roll. When << 1 then f 1, and when>> 1 then f -‘ 0. Thus, the amount of smoothing depends upon the specific parameters ofa given situation. If the smoothing is significant then the conventional (ideal) dual referencetrack interpretation of the differential phase will be less sensitive to sinusoidal roll coupledwith unknown terrain.Translational MotionIt is now helpful to evaluate how translational motion of the aircraft conthbutes to thedifferential phase error. Beginning with the change in the error in the off-nadir angle to thetarget:/cbab(t) sin (crab + Oao + (t)) sin ((1)) (4.63)62and differentiating, yields:________(4 64)-Expanding in an Taylor’s series yields:ab(t) ab(O) + — O(O)) (4.65)By considering displacements M from the center of aperture position one obtains:O(t) — M(O) M108X (0) (4.66)This yields a differential phase error that varies along the aperture governed by the translationaldeviation of the antenna baseline:2irbjj08M108(t)Z.O(O)2r jjosMi08(t)hR2sinOIf this analysis is repeated for changes in the assumed off-nadir angle to the target:/..(I)ab(t, 1) jii (ai + Oao(t) + Sjfl () (4.68)differentiating yields:ö4ab(t, 1) 2Kb108h 4öOa Sifl 0Expanding in an Taylor’s series yields:1ab(t) 1çbab(0) + 0a(t)- 0a(0)) (4.70)By considering displacements M from the center of aperture position one obtains:Oa(t) — 0a(0) Mo3(t) (4.71)63This yields a differential phase error that varies along the aperture governed by the translationaldeviation of the antenna baseline:(i) 2?rbj0Mj6(t)h (4.72)AR smOCombining these two error sources gives:A +b±103M(t)) (4.73)It is easy to show that:bi08iW i(t)+bjj03Mio8(t) = bM(t) Sill(0iqos +0Mlos)(4.74)bM(t)where the angles are measured from the line-of-sight direction. Therefore, the effectivedifferential phase error along the aperture is:( 2rhbM(t)AR smOWhen linear and quadratic motions are analyzed for differential shifts and phase errors,the results are negligible for realistic aircraft motions and rugged terrain. This is presented inAppendix H. Therefore, the differential phase errors along the aperture caused by the unknownterrain coupled with translational motion pose no problem. It should be noted that the effects,such as defocussing, observed for the single channel case will still apply to both channels inthe InSAR case. Therefore, the predominant effect from translational motion is the azimuthshift and defocussing that occurs essentially identically in both channels.4.4.4 FM Rate Induced Differential Phase ErrorsAn analysis of the effect of motion compensation induced FM rate errors for a singlechannel was made in the previous chapter. A quadratic phase error was observed along theaperture that was proportional to the constant offset displacement between the antenna positionand the position used to generate the matched filter (Equation 3.28). In the interferometric case,64the displacement to the track used to generate the matched filters may be different betweenchannels. This leads to a different quadratic phase error in each channel. This differentialphase error along the aperture is:2KVtLRab 476ARwhere /Rab is the difference in the displacements to the track(s) used to generate the filtersor equivalently the inter-channel difference in the slant range matched filter errors. Using theresults from [27] it is known that the differential phase bias error is about 1/3 of the magnitudeof the differential quadratic phase error at the edge of the processed aperture. Therefore, onecan estimate the differential phase error at the peak:1 2Kv2()/Rabbias AR2 (4.77)If processing to the 3 dB azimuth beamwidth is performed,T = T3dB = (4.78)Dvand using the sensitivity in height estimate to differential phase errors result from Equation4.24, one can obtain the bias in the height estimate:A2RLRab Sfl 0/hbias— 2 (4.79)l2bD Slfl(O+Oab)This expression can be used to limit the difference in the displacements for matched filtergeneration. If typical CCRS parameters are used (Table 1), and a 1 m height bias error isallowed, one finds that the difference in the displacements must be:/Rab < 1.3 m (4.80)This gives a good indication of how often each matched filter must be updated in order to keepthe difference in the quadratic phase errors between channels very small. Clearly, the matchedfilters must be updated for every azimuth line across the range swath.65In addition to this, the inter-channel resampling that is performed prior to motion compensation will result in data along azimuth that requires different matched filters. If block azimuthprocessing is performed then the length of the block is limited by the amount of flight motionalong the block. Errors in the matched filter along azimuth that are different between channelswill result in differential phase biases according to the above expressions.4.4.5 RCMC Induced Differential Phase ErrorsIn the previous chapter the effect of motion compensation induced errors when RCMCwas performed for a single channel was described. It was observed that a quadratic phaseerror along the aperture exists that is proportional to the constant offset displacement betweenthe antenna position and the reference track (Equation 3,37). In the interferometric case, thedisplacement to the track used to apply the phase correction can be different between channels.This leads to a different quadratic phase error in each channel. The resulting differential phaseerror along the aperture is:(t) 2Kvt d±i0 (4.81).AR tanO 2where /.djios is the difference in the displacements from the reference track(s). The relativedisplacement is divided by two because of the one-way path difference between the two channels(channel B is receive-only). Using the results from [27] it is known that the differential phasebias error is about 1/3 of the magnitude of the quadratic phase error at the edge of the processedaperture. Therefore, one can estimate the differential phase error resulting from the differencein the displacements to the track(s):‘ 2m2 ,jA L7rv I t_UJj03-biasc” 12,\R tanO 2 (4.82)If processing is to the 3 dB azimuth beamwidth and using Equation 4.24 one can obtainthe corresponding height estimate bias:A2 R cos 0 dji08hbias— 24bD Sin (0+ aab) (4.83)66This expression can be used to limit the difference in the displacements from the referencetrack when RCMC is performed. If typical CCRS parameters are used (Table 1), and a 1 mheight bias is allowed, one finds that the difference in the displacements must be less than:/.dj03 <3.5 m (4.84)For the single reference track case the difference in the displacements to the reference trackare on the order of the baseline length. It would appear that this RCMC induced error maybe significant for the single reference track case.In the dual reference track approach the displacements from the reference tracks will bevery similar in the two cases, greatly reducing the sensitivity to the RCMC errors. Referringto Figure 12 the effective difference in the displacements between the two channels can beestimated from the difference in the applied phase correction:— = = (Rtao — Rib0) + (Rbo — Rao)) (4.85)Expanding the slant range values in binomial series yields:bcos(Oao(t) + aab(t)) — bRTcos (o + aRT) (4.86)The antenna baseline is naturally chosen to be the same as the reference track baseline, thus:zd 1,—a-- 2bsrn(0aRT + aab(t) + °ao + Oao(t)1 (4.87)(—aRT + aab(t) O + O0(t))Defining:Laab(t)= cYab(t) — aRT (4.88)and approximating:11 [ , JJosVaoI) t’ao67which gives:( dj0‘b03 j/.\ab(t) + (4.90)For typical offsets from the reference track and roll of several degrees this effective differencein the displacement from the tracks is negligible. It is clear that if RCMC is to be performed,the dual reference track approach should be used.4.4.6 Effects of No Motion Compensation ResamplingThe idea of neglecting to resample during motion compensation for the single channel casewas introduced in the previous chapter. It is now important to extend this to consider whatdifferential phase errors might occur. More specifically, there may be a slight difference inthe effect on the PSF for each channel that decreases the correlation. If this difference can becharacterized as a differential residual range migration then one can use the results from [27] toevaluate the correlation coefficient. Appendix I considers the difference in the range migrationerrors between channels due to neglecting to resample the data during motion compensation.The resulting differential migration error along the aperture is:zR(t) bjjosMjjos(t) (4.91)where Mjios (t) is the linear deviation over the aperture from the center of aperture position.For a linear flight variation leading to uncompensated differential linear range migration,the inter-channel correlation coefficient can be shown to be [27]:it/3= 2___= Si(!)(4.92)where 3 is the range walk at the band edge in fractions of a range resolution cell. Fortypical CCRS parameters, the effect on the correlation resulting from numerical integration isnegligible (7 > 0.99).68For quadratic flight motion leading to uncompensated differential quadratic range migration,the correlation coefficient can be shown to be [27]:fr = __fsinc(a2)da (4,93)0where e is the uncompensated quadratic migration at the band edge. Again, for typical CCRSparameters the effect on the correlation resulting from numerical integration is negligible(-y > 0.99).404.7 Effects of Inertial Data ErrorsInertial data errors can have a significant impact on the accuracy of InSAR. Errors in thebaseline length or angle translate directly to height estimation errors. In the single channelcase it was shown that errors in the assumed position of the antenna led to phase correctionerrors during motion compensation. In the clutter data case the compression operation doesnot smooth the phase noise like it does in the point target case. The compressed impulseresponse will be broadened and have phase errors. In the InSAR case, the effect of INS errorson the height estimate can be calculated by differential analysis. The differential analysis forthe general case of errors in baseline length and angle has been performed [10]. This workalso applies for the conventional dual track approach with INS errors.The single reference track approach and new dual track method result in the samesensitivities to baseline length and angle as previously reported [10]:RsinO (494)0b btan(O+crab)andRsinO (4.95)öaabIt is therefore straightforward to calculate the effect of INS bias errors on the height estimates.The bias errors will have a small impact at the motion compensation stage, but, can causesignificant errors in the estimated elevation through the above equations.69High frequency INS errors will translate into high frequency differential phase errors. Theeffect on the height estimate can be estimated from Equation 4.24. In addition, there will bedefocussing that will increase smoothly with increased phase noise. Provided the statistics ofthe INS errors are about the same in both channels then the broadening will be about the same.Differential broadening will lead to an increased height estimate variance through Equation 4.8.Given that all three methods, conventional double, new double, and single, all have thesame sensitivity to INS errors there is no preferred approach to handle these errors. In allcases, differential phase smoothing will help reduce the effect of high frequency INS errors.Bias errors are much more difficult to deal with. In mostcases control points will be necessaryin order to obtain height estimates on the order of one meter accuracy.4.4.8 Summary of Errors and Constraints for InSARIn ideal motion compensation all the required geometric parameters are known and thereis no restrictions placed on the processing. This leads to accurate compensation except thatthe range varying phase correction induces a range spectral shift. This can be minimizedby keeping the reference track close to the antenna position by segmenting the track. Idealcompensation is not possible in the InSAR case because the terrain elevation is not known.In realistic InSAR processing situations non-ideal motion compensation must be used. Thefollowing tables summarize the resulting errors and limitations given the assumptions made atthe beginning of the chapter.• The range varying phase correction causes a range spectral shift as in the single channelcase. Subsequent processing, such as inter-channel registration, can lead to decorrelationif baseband data is assumed.• The coupling between the RCMC and the range varying phase correction, which can bedifferent between channels, can lead to differential phase biases.70• Differences in the errors in the FM rates between channels can lead to differential phasebiases.• The effect of unknown terrain that is common to both channels is the azimuth peak shiftand azimuth defocussing seen in the single channel case.• The differential effects of the unknown terrain can lead to differential phase biases whenangular acceleration occurs.• Most other effects will be almost identical in both channels and thus cancel out.• Errors in the inertial data can lead to differential phase errors and biases in the heightestimates. The high frequency INS data errors can be reduced by filtering the interferogrambut INS bias errors will remain.Table 6 outlines the main differences or similarities between the single and dual referencetrack approaches. It is clear that each approach has certain advantages and disadvantages.The issues compared are:• Whether the differential phase is continuous if the reference track(s) are segmented.• Whether there are differential range spectral shifts caused by the range varying phasecorrection (RVPC). The effect of this will depend upon subsequent processing.• Whether there are biases caused by the coupling effect of RCMC and the range varyingphase correction (RVPC).• Whether biases due to slowly varying aircraft drifts coupled with unknown terrain can becompensated for.• Whether there is a range dependent fringe in the interferogram phase even for flat terrainat the reference level.71Table 4 InSAR Motion Compensation Errors and Limitations #1Physical Effect Error or LimitationDisplacement Offset Range Oversampling Factor Or:(range varying phase correction)2cIdjj05l+ARtOBW <ORCMC Induced Bias Errors Differential Phase Error Along Aperture:(range varying phase correction)2rvt /dj05L4(t),\R2tanO 2Differential Phase Bias Error:-biairv2T/djj03S l2ARtanOHeight Estimate Bias (processing to 3 dB):A2 R cos Ozd±105Lhbias 24bD sin (0 + aab)FM Rate Induced Bias Errors Differential Phase Error Along Aperture:27rVt/Rab(t)AR2Differential Phase Bias Error:TV2RabL4bas 6AR2Height Estimate Bias (processing to 3 dB):A2R/Rab sinOhbias— l2bD sin (0 + aab)72Table 5 InSAR Motion Compensation Errors and Limitations #2Effect of Unknown Terrain Error or LimitationEffects Common to Both Channels Phase Error Along Aperture:47rhdjj03(t)lerr(t) AR sinGAzimuth Peak Shift:h Od±102 (sec)V Sifltl at5% Defocussing:XR sinGajj08 < hT2Processing to 3 dB:,\ sin OD2vajj03 <Differential Effects Differential Phase Error Along Aperture:2irhbjj0(t))RsinODifferential Phase Bias due to AngularAcceleration:a0iiirbi8T2bias 12R sinGHeight Estimate Bias due to AngularAcceleration (processing to 3 dB):Lhbias_______________________h 24bvD (0 + ab)73Table 6 Comparison of Single and Dual Reference Track Motion CompensationProperties Single Dual1 Continuity of Differential Phase with Segmented Reference Track(s) Yes No2 Equal Range Spectral Shifts between channels due to RVPC No Yes3 Insensitive to Bias Errors caused by RCMC and RVPC No Yes4 Compensation for Unknown Terrain coupled with low freq. motion Yes Yes5 Removal of flat earth fringes Yes No4.5 Proposed Motion Compensation ApproachesGiven the analysis of the interferometric motion compensation problem presented above,it is now possible to define recommended approaches. This section will begin with a reviewof the system requirements followed by a discussion of the possible processing approaches. Anew approach will be proposed which combines the benefits of the dual and single referencetrack methods.4.5.1 Realization of RequirementsA general purpose InSAR system will require the following in order to generate topographicmaps accurately and efficiently:1. Differential phase fidelity. This means that all bias errors are minimized and randomerrors are reduced by spatial filtering. Optimum smoothing is obtained when the samplesare uncorrelated. This may require full resolution azimuth processing making RCMCnecessary.2. Capable of mapping rugged terrain (+1- 1 kin) under typical flight conditions. This willlikely require segmentation of the reference tracks to follow the flight motion.3. Continuous differential phase over the entire image scene is desirable. This means thatonly one tie point is required to unwrap the phase of the entire image.74The ability to meet the above requirements can be evaluated by considering the summary oferrors and limitations of motion compensation in the previous section. A quantitative analysiscan be carried out for a specific InSAR system. But, a qualitative approach can be formulatedfor the general case.Typical CCRS System ApproachFor the typical InSAR scenario for the CCRS system a fairly simple approach to motioncompensation is adequate. This would include:1. No motion compensation resampling,2. Phase correction to a single reference iTack that may be segmented, if necessary.3. No RCMC.4. Azimuth matched filters determined separately for each channel based on the referencelevel positions.5. Phase unwrapping that assumes continuous differential phase (as far a motion compensationis concerned)6. Differential phase to elevation conversion based on the new fonnulation (Section 4.4.1).Apart from the possible defocussing problems for severe terrain and translational accelerationthis approach should yield a motion compensation approach which contributes on the order ofa meter to the height estimate error budget.Most Accurate ApproachThe most accurate approach would be to implement the following:1. Perform motion compensation resampling.2. Phase correction to dual reference track that may be segmented if necessary.3. Perform RCMC.754. Azimuth matched filters determined separately for each channel based on the referencelevel positions.5. Phase unwrapping that assumes discontinuous differential phase (as far as motion compensation is concerned).6. Differential phase to elevation conversion based on the new dual track formulation (Section4.4.2).The main problem with this approach is the discontinuous differential phase as a result ofsegmentation of the reference tracks. A modification to this approach will now be presentedthat combines the benefits of the single and dual reference track approaches.4.5.2 Combination of Single and Dual Reference TracksA new approach to InSAR motion compensation will now be proposed which is a simple,yet effective, extension to the dual reference track approach. It is clear that performing dualreference track motion compensation is preferred for the compression operation, but, it isdesirable to have the single reference track differential phase for phase smoothing and phaseunwrapping. The proposal is to convert back to the single reference track form after dual trackmotion compensation, optional RCMC, and azimuth compression have been performed.This conversion can be accomplished by a post-compression phase correction that convertsthe channel B reference track to the channel A reference track. Given that the reference trackpositions are parallel over the segment, and the position of each range cell on the groundis fixed if motion compensation resampling is performed, this means simply adding a rangedepending phase constant to each azimuth matched filter. The resulting differential phase ofthis two step approach is equivalent to performing motion compensation in one step to a singletrack. In other words, the differential phase will be zero if each step is performed properly.Again, because of the unknown terrain, errors will arise and the differential phase will be76given by these errors. The combined errors from each step are in fact the same as if onestep of single track compensation were perfonned, except that the RCMC induced phase bias(characteristic of single track motion compensation) will not be present (Equation 4.83). Thiswill now be described in detail,From Equation G.9 of Appendix 0 and Figure 12, the differential phase error resultingfrom the unknown terrain after dual track motion compensation, optional RCMC, and azimuthcompression is:Ed(A) — Ed(B) = Rao— Rtah— Rtb0 + Rbo + Rtbh — Rbh) (4.109)The next step is to apply a phase correction which makes the channel B reference track dataappear as though it were at the channel A reference track. The applied phase correction tochannel B is:qB(applied) = (Rtao—(4.110)Notice that this is constant along the aperture if motion compensation resampling is performedand can therefore be applied as a range varying constant phase term for each azimuth matchedfilter. Otherwise, it must be applied pixel by pixel on the compressed image or to theinterferogram. Due to the unknown terrain the ideal correction should have been:cB(ideal) = Rtah — Rtbh) (4.111)Therefore, there is a differential phase error from this operation which is:E8(B) = Rtao— Rtb0 + Rtbh — Rtah) (4.112)The resulting differential phase after the two operations will be zero if all corrections wereideal. Therefore, the differential phase is equal to the residual errors:= Ed(A)— Ed(B) — E3(B) (4.113)77which results in:=— Rbh) (4.114)which is precisely the same result as if one step single track compensation had been applied(Equation 4.15), except without the RCMC induced errors.Most Accurate and Most Convenient ApproachThe resulting most accurate and convenient approach would therefore be to implement thefollowing:1. Perform motion compensation resampling.2. Phase correction to dual reference tracks that may be segmented if necessary.3. Perform RCMC.4. Azimuth matched filters determined separately for each channel based on the referencelevel positions. Apply the extra phase constant term to each matched filter for ChannelB (Equation 4.110).5. Phase unwrapping that assumes continuous differential phase.6. Differential phase to elevation conversion based on the new dual/single track formulation(Equation 4.114).If motion compensation resampling is not performed the extra phase correction must be appliedafter azimuth compression on the entire image. The benefits versus the extra computation wouldhave to be weighed. The validity of such approaches to motion compensation will be verifiedby simulation and processing of real InSAR data in the next chapter.78Chapter 5 Evaluation of Non-IdealInSAR Motion CompensationIn order to ensure that the modelling and analysis of the previous chapters is correctexperiments were performed using a point target simulator as well as processing real InSAR datafrom the CCRS Convair 580 system [25]. The simulation and real data processing used the zero-Doppler pulse processing and flat earth assumptions of non-ideal InSAR motion compensation.This chapter will begin with a description of the point target simulator and the correspondingresults and conclude with results from real InSAR data processing.5.1 Point Target SimulationA comprehensive simulation of the InSAR data acquisition and processing mechanism wasperformed. Dual channel range compressed raw SAR point target data was generated for anydesired flight motion and target location. An InSAR processor was then used to process thedata into an interferogram and finally into a height estimate for the target. In this section, thesimulator will first be described followed by the modelling of typical flight motion. The resultsobtained using the modelled motions will be described as well as the results from using realmotion data from the CCRS Convair 580 aircraft. When the predicted results are compared tothe simulation results good agreement is found.5.1.1 The Point Target SimulatorThe point target simulator consists of three parts. First, range compressed point target echodata is generated for any aircraft trajectory, including roll effects, and for any target elevation.Next, both channels of the raw data are processed by the InSAR processor. This consistsof motion compensation phase correction, azimuth compression and two-fold oversampling,interferogram creation, phase unwrapping, and differential phase to elevation conversion. The79third component of the simulator is the evaluation routines that monitor the results alongthe way. This consists of performing point target analysis on the two compressed channelsand on the interferogram, as well as analyzing the final height estimate. The simulator wasimplemented in the MATLAB simulation language. The simulator is based on typical C-bandairborne systems, but the specific parameters were taken from the CCRS InSAR system for theConvair 580 aircraft. The simulation algorithm will now be described in more detail.Simulation AlgorithmThe simulation algorithm consists of the following processing steps:1. Generate range compressed point target data for a particular flight and target position fortwo channels (one transmit receive, one receive-only) assuming receive-only channel isregistered with main channel assuming flat terrain. The registration is implicit in thesampling of the raw data. The data is 2048 samples in azimuth and 21 samples in range(variable).2. Apply phase correction to each channel assuming flat terrain and using either one or tworeference tracks. No motion compensation resampling is performed.3. Perform azimuth Fourier Transform.4. Optionally, perform RCMC in the frequency domain by baseband interpolation for eachchannel. An 8 point FIR kernel is used and interpolation is made to 1/16 of a sample.5. Performed azimuth compression using a matched filter for each channel. The slant rangeused to generate the matched filter is correct at the zero Doppler pulse. Terms up to fourthorder in time are used.6. Perform azimuth Inverse Fourier Transform including two-fold zero padding to prepare forinterferogram creation. This is required as interferogram generation doubles the bandwidth.807. Interpolate by a factor of 16 and analyze the compressed peak characteristics in eachchannel.8. Generate the interferogram using the original (not interpolated) compressed data from eachchannel.9. Interpolate by a factor of 16 and analyze the compressed peak characteristics in theinterferogram.10. Unwrap the differential phase using the terrain elevation knowledge.11. Convert the differential phase into a height estimate.5.1.2 Modelling of Typical Aircraft MotionsIn order to use the simulator, aircraft flight motions are required. Antenna displacementdata obtained from CCRS for the Convair 580 was used to come up with classes of motions thatcould be compared with the predicted results (Figures 14 and 15). Data from a number of flightswas analyzed. From this data the following motions were selected for use in the simulator:Table 7 Typical Flight Motion for CCRS Convair 580Motion Type Peak Flight MotionTranslational Velocity 0.5Translational Acceleration 0.01 gBaseline Angular Velocity (roll) 0.2Baseline Angular Acceleration (roll) 0.3Sinusoidal Baseline Variation (roll) 0.15 deg peakDisplacement Offset 10 m5.1.3 Simulation Results of Modelled MotionsThe typical peak flight motion shown in Table 7 were used, one at a time, in thesimulator using the parameters of Table 1 (unless otherwise specified). The motions weredivided into separate components in the line-of-sight direction and perpendicular to this toallow the theoretical error analysis to be evaluated. Both single and dual reference track81Eci)-Jc,)LicciC20zE2ci)Sci)0cci0.C’,00)ci)C0ccicci>0Figure 14 Convair 580 Flight Displacement DataFigure 15 Convair 580 Roll VariationTime (sec)50Time (sec)82motion compensation was performed. Apart from the differential phase errors the results wereessentially identical between the single and dual track approaches. The dual track results willbe presented first, followed by the single track results.Dual Reference Tracks with RCMCTable 8 shows the results for the dual reference track approach for target elevations lessthan or equal to 1000 m from the reference track and for motions that yielded insignificanterrors. In addition, there was no measurable inter-channel shifts. This indicates that theseflight motions are acceptable for the simulator algorithm (including the assumptions containedthere such as the center of aperture assumption, neglecting to resample, and non-ideal motioncompensation). The results for targets below the reference track are similar except for thechange in the angle to the target which makes the above/below results slightly asymmetric.Table 8 Dual Reference Track Motion Compensation with RCMCFlight Motion h <= 1000 mzh Azimuth Azimuth Range Inter-channel(m) Broadening Shift Broadening Broadening(%) (m) (%) (%)llwith <0.05 <1 0 <1 <1d<=lOm4LQs_O5!!! <0,05 <1 0 <1 <1a108 = 0.Olg < 0.05 < 1 0 < 1 < 1<0.05 <1 0 <1 <1The most significant results are shown in the following figures. Figure 16 shows theresulting azimuth shift from translational velocity perpendicular to the line-of-sight direction.The predicted shift of Equation 3.48 agrees well with the simulation result.834.504.003.50E3.00ci)2.500E 2.001.501.000.500.00Velocity Induced Peak Position Shift0.00 0.20 0.40 0.60 0.80 1.00Target Elevation Above Reference Level (km)Figure 16 Unknown Terrain with vjj08 = 0.5 ?. and R 10 kmFigure 17 shows the resulting azimuth defocussing occurring in both channels for a uniformtranslation acceleration perpendicular to the line-of-sight of 0.01 g. By using Equation 3.55there should be about 5 % defocussing for a terrain elevation of 500 m. This agrees quitewell with the simulation results.84Acceleration Induced Azimuth DefocussingTranslational Acceleration = 0.OlgTarget Elevation Broadening(km) (%)0 00.2 00.4 20.6 60.8 151.0 36Figure 17 Unknown Terrain with ajj08 = O.Olg, and R = 10 kmTable 9 shows the results from an angular acceleration of 0.3 for a range of slantrange values. Predicted differential phase errors using Equation 4.52 and height estimatebiases using Equation 4.54 were compared with the simulation results (using no antenna beampattern) yielding good agreement. The larger discrepancy at near range is due to the errors inestimating the effective baseline perpendicular to the line-of-sight as it gets smaller.Table 9 Effects of Angular Acceleration, arc ii = 0.3 , h = 1 km, and no antenna patternSlant Range Processed Mtheory 4)sim htiieory Lhsim(kin) Aperture (mrad) (mrad) (m) (m)(sec)10 3 19 26 -0.5 -0.715 3 14 15 -0.7 -0.715 4.6 32 34 -1.6 -1.620 4.6 30 31 -2.1 -2.120 6 53 55 -3.8 -3.785Figures 18 and 19 show the results from sinusoidal roll of 0.15 deg peak and withperiod = 3 sec. As was suggested by Equation 4.62, the effect of the phase bias error beingabout equal and opposite to the change in the expected phase, and the smoothing performed bythe compression operation, combined to result in a very smooth measured differential phase.The expected phase followed the roll motion thus leading to a sinusoidally varying heightestimation error. This smoothing effect was even more pronounced at far range due to thelonger aperture processed. It should be noted that using a smoothed version of the inertial datato convert the differential phase to elevation will not solve this problem. What will happenis that new height errors will occur even for targets at the reference level. At least underthe current approach it is only targets that deviate appreciably from the reference level whichare affected.CaI-a)Co(aaCla)b1.Sinusoidal Roll and Unknown TerrainSinusodaI Offset (sec)Figure 18 Effect of Sinusoidal Roll on Differential Phase(0.15 deg peak, period = 3 sec, h = 1 km, R = 10 km)866I0.20.10-0.1-0.2-0.3Sinusoidal Roll and Unknown Terrain‘0 0.5 1 1.5 2 2.5Sinusoidal Offset (sec)Figure 19 Effect of Sinusoidal Roll on Height Estimate(0.15 deg peak, period = 3 see, h = 1 km, R = 10 km)Single Reference Track with RCMCThe only difference between the results from the single and dual reference track approacheswas that involving an extra differential phase bias in the single reference track case. This waspredicted by Equation 4.82 to occur when RCMC was performed. Figure 20 shows how thetheory and simulation results agreed very well.5.1.4 Simulation Results of Actual Motion DataActual motion data from Figures 14 and 15 were also used in the simulator. This datacorresponds to pass 6 of the Kananaskis flights in February 1992 of the Convair 580. Pointtarget simulation results were obtained along the flight line every 0.2 seconds and dual referencetrack motion compensation with RCMC and single track without RCMC were performed. Thisis equivalent to simulating a strip of point targets along azimuth separated by about 25 m. TheI Ioo00+01 104 I- -0 0.1-0 0---°-----I I I II I II I- 4-10 I I 010 I I 0o I 001 10 I-t t- -10 0 I8720wC0ce2CoUIa)ISingle Reference Track with RCMCTheorySimulation0.90.80.70.60.50.40.30.20.12.0 3.0 4.0Processed Aperture (sec)Figure 20 Single Track Bias Errors, R = 10 km, and no antenna patternfollowing results indicate the errors that would results from typical Convair 580 flight motioncoupled with unknown terrain and precise inertial data measurements.Dual Reference Tracks with RCMCThe results from dual reference track motion compensation with RCMC and real motiondata were similar to those using the modelled motion data. Most effects were almost identicalin both channels. The significant results of differential phase errors, height estimation errors,azimuth defocussing, and azimuth shifts will now be described.88Differential Phase Errors The dominant differential phase error was due to the smoothingeffect of unknown terrain coupled with aircraft roll (Figure 21). The expected differential phaseitself followed the large scale motion of Rtb0 and Rtao with a small scale ripple due to theaircraft roll. The measured differential phase followed the large scale motion but smoothed outsome of the ripple. This leads to a ripple in the height estimate.-DCuEa)U)aCuCa)a)0Measured and Expected Differentia’ Phase100Time (sec)Figure 21 Measured Versus Expected Differential Phase Dual Track (R = 10 km, h = 1 kin)Height Estimation Errors Given that the differential phase errors follow the roll motion,then the resulting height estimation errors will also follow the roll motion. Figure 22 showsthe resulting height estimation errors using the new double reference track interpretation ofthe differential phase.89Height Estimation Error25 30Time (sec)Figure 22 Height Estimation Errors and Roll (R = 10 km, h = 1 km)Figures 23, 24, and 25 compare the height estimation errors for the conventional (idealassumed) and new dual reference track approaches for different slant range values. Thereference tracks were chosen to minimize the possible biases that can occur in the conventionalapproach. At near range the errors resulting from the aircraft roll are small in the new approachcompared to the residual biases of the conventional approach. At mid and far range the errorsare about the same size but the new approach has higher frequency errors due to the roll. Anexplanation for this was given above for the case of point targets and sinusoidal roll. Theimpulsive height errors come from regions where extreme defocussing occurred.50SCwC0EU)w-0.10)ci)90Conventional Dual Track ApproachE0wC0E-C=E0LiiC0clELii0)ci)=Time (sec)Figure 23 Height Estimation Errors (R = 10 km, Ii = 1 km)Time (see)Figure 24 Height Estimation Errors (R = 15 km, h = 1 km)91New Dual Track ApproachConventional Dual Track ApproachNew Dual Track ApproachConventional Dual Track Approach300Figure 25 Height Estimation Errors (R = 20 km, h = 1 km)Azimuth Defocussing The theory predicted that acceleration perpendicular to the line-of-sight would lead to defocussing in both channels. Figure 26 shows the resulting 3 dB widthazimuth broadening for a target elevation of Ii = 1 km above the reference track. Also plottedis the general variation of the displacement from the reference tracks perpendicular to theline-of-sight direction. The total variation over the processed extent was about 6 m. At afew positions along the flight line significant defocussing occurs. Figure 27 shows the sameexperiment but with h = 500 m. By halving the target elevation the defocussing is reduced toreasonable levels. In either case, the inter-channel differential defocussing is less than a fewpercent provided significant defocussing does not occur. When significant defocussing occursthe interpretation of a 3 dB width loses its meaning as the sidelobes have joined the main lobe.New Dual Track ApproachTime (sec)92CCa)-D0CD4-.CDC).4-.ENUnknown Terrain and Real Motion DataTime (sec)Figure 26 Azimuth Broadening with Unknown Tenain (R = 10 km. h = 1 km)300Azimuth Shift Figure 28 shows the resulting mis-positioning of the peak in azimuth for atarget with Ii = 1 km. The results were the same across the range swath. Using Equation3.48, and by low pass filtering the motion data with a 1 second long box-car filter to estimatethe velocity, the predicted shift was calculated. It is clear that the prediction and measuredresults agree.Single Reference Track without RCMCThe above simulations were repeated but with single reference track motion compensationand no RCMC was performed. Most of the results were essentially the same as the dual0 50 100 150 200 25093Unknown Terrain and Real Motion DataTime (sec)Figure 27 Azimuth Broadening with Unknown Terrain (R 10 km, h = 500 m)reference track with RCMC results. The exceptions to this were the resulting differential phaseand a small bias error. Figure 29 shows how the expected differential phase follows the rollmotion but is otherwise unaffected by the larger drifts. The corresponding height estimationerrors are shown in Figure 30. The small bias at the peak is due to a coupling between therange phase ramp in the interferogram PSF, caused by the difference in the range varying phasecorrections between channels, and a small migration of the compressed peak across range dueto the uncompensated range migration. If the phase at the correct slant range (slant range ofpeak if RCMC had been performed) is evaluated, no bias is found. This bias does not occurfor dual reference track motion compensation without RCMC because there is no range phaseI II I II II I II I II I II II II II I20 -18161412108642C00)Cci)004-DEN-4.I 1Uiios,/I /I’I—I-I‘-FIJ/III/II I——tI50 100 150 200 250 30094IUnknown Terrain and Real Motion Data100Ti me (Sec)Figure 28 Azimuth Shift with Unknown Terrain (h = 1 km)ramp in the interferogram. This error can be minimized for the single reference track case bykeeping the processed bandwidth small.The azimuth broadening and azimuth shift are slightly smaller than the dual track withRCMC case due to the shorter effective aperture processed.Processing with Segmented TracksIn order to investigate the continuity in the differential phase of the new dual/single trackapproach the point target simulator was used. The segmented tracks are handled by overlappingthe azimuth data sufficiently to ensure a constant reference track for a given compressionoperation. The point target simulator was run using the real motion data for two cases of0 10 20 30 40 50 60 70 80 9095Measured and expected differentia’ phase100Figure 29 Measured Versus Expected Differential Phase Single Track (R = 10 km, h = 1 km)reference tracks. The resulting two strips of data were then cut and pasted together at themidway point. This was performed for the dual reference’ track approach and for the newdual/single track approach where a phase correction is applied after compression to convertfrom the dual track to the single track formulation. Figure 31 shows the resulting differentialphase of the interferogram. It is clear that the dual track approach results in discontinuousdifferential phase. The new dual/single track approach has continuous differential phase acrossthe boundary. When the height estimates were compared between the two approaches theresults were identical, as expected (Figure 32).5.1.5 Discussion of Point Target ResultsThere is a convincing correlation between the theoretical error modeffing presented in thisTime (sec)96E011w0E4-Cl)w4-a)=Height Estimation ErrorTime (sec)300Figure 30 Height Estimation Errors Single Track no RCMC (R = 10 km, h = 1 km)work and the simulated results using the point target simulator. As predicted, the followingmain results were observed:1. The majority of errors resulting from non-ideal motion compensation cancel out betweenthe two channels leading to accurate height estimation.2. Unknown terrain coupled with velocity perpendicular to the line-of-sight direction leadsto an azimuth shift of the compressed peak in both channels identically. Provided thisis only a few meters the subsequent differential phase smoothing by spatial filtering willreduce the impact of this error. Alternately, if the error is large the shift can be predictedonce a height estimate is available and corrected by a resampling operation. The predictiveability of the theory is very good. But, this correction should not be necessary for theCCRS InSAR system.0 50 100 150 200 25097ccici)U)Cci0Cuci,0-1.5ECociCuEU)wci,Differential Phase Across Segmentation Boundary1.51 I -2Dual Tracc0.5--Dual! SinsBaseline UBaseline AlBaseline Lngth = 2.819 [nBaseilne Aigle = 40.74 dg._____________Ile Trackngth = 2.805gle=40.24d•11ag.Figure 31Th 50 100 150 200 250 300Time (sec)Differential Phase with Segmented Reference Tracks (R = 10 km, h = 1 km)Dual and Dual/Single Track Height Estimates150 100 150 200 250 300Time (sec)Figure 32 Height Estimates with Segmented Reference Tracks (R = 10 km, h = 1 krn)983. Unknown terrain coupled with acceleration perpendicular to the line-of-sight direction leadsto similar defocussing in both channels. For rugged terrain significant defocussing can occurat a few positions along the flight line for the CCRS InSAR system. Only knowledge ofthe terrain elevation and complete time domain motion compensation and processing canreduce this effect4. The single reference track approach with RCMC has an inherent differential phase biasdue to the difference in the displacements from the reference track coupled with RCMC.This bias did not occur in the double reference track approach, as expected.5. The single reference track approach without RCMC also has a small bias error in the heightestimate at the peak due to the difference in the displacements from the reference trackcoupled with the migration of the compressed peak due to uncompensated RCM. This biasdid not occur in the double reference track approach.6. The dominant source of differential phase errors, and thus height estimation errors, is due toaircraft roll coupled with unknown terrain. The roll motion produces angular accelerationswhich cause bias errors. In addition, the motions are often oscillatory allowing forsmoothing to occur due to the compression operation. This resulted in smoothed measureddifferential phase which did not follow the higher frequency expected differential phase.The conventional approach to interpreting the differential phase is not as sensitive to thiseffect as the new single and new dual track approaches. But, the conventional approach issensitive to the choice of reference tracks for a given segment.5.2 Processing CCRS InSAR DataIn addition to the successful results from the point target simulations, processing with realCCRS InSAR data was carried out to demonstrate the ability to obtain quality interferogramsfor the various proposed approaches. In this section, the InSAR processor will be described99followed by a description and display of the resulting interferograms and a slant rangetopographic map. Some additional experiments with inertial data errors will also be described.A basic InSAR processor was implemented in the MATLAB simulation language followingthe process flow of Figure 2. The main processing steps included:1. Resample channel B to register with channel A at the reference level.2. Perform motion compensation phase correction to each channel using single or dualreference tracks.3. Perform azimuth compression using fast convolution without RCMC.4. Form interferogram.5. Perform phase smoothing and phase unwrapping.6. Convert differential phase to elevation.The processor was designed to operate on one block of range data at a time. The blocks couldthen be pieced together to form a complete image. The processing parameters were those ofTable 1.5.2.1 Three Hills Data ResultsThe data set processed was of the Three Hills area in Alberta. The data was acquired by theCCRS InSAR system in February 1992, pass 10. The scene is of farm land with rolling hills.The raw data was 12000 range lines with a swath width of 2048 samples. The interferogramswere phase smoothed by sample averaging of 10 samples in azimuth and down sampled by 10producing a final image 2048 samples in range by 1024 in azimuth. This yielded about equalsampling intervals of 4 m in each direction.Figure 33 shows the interferogram square root of the magnitude image which is essentiallyequivalent to a single channel SAR image. The image looks distorted because of the slantrange perspective.100Figure 34 shows the resulting interferogram differential phase for the dual reference trackcase. One cycle of the color wheel corresponds to a change of 2w radians. The dominantvariation is due to the flat earth fringes which depend upon the position of the reference tracks.The discontinuity in the reference tracks at range line 512 will clearly cause phase unwrappingproblems. The flat earth fringes have to be removed before the variation due to topographywill be visible in the interferogram.Figure 33 Square Root of Interferogram Magnitude Image101Figure 35 shows the resulting interferogram differential phase for the new dual/single trackcase. One cycle of the color wheel corresponds to a change of r radians. The extra phasecorrection was applied to the phase smoothed and down sampled interferogram to reducecomputation. The new duallsingle approach has the expected continuous differential phaseacross the segmentation boundary.Figure 36 shows the interferogram differential phase for the single reference track case.One cycle of the color wheel corresponds to a change in differential phase of ir radians. ItII U 1111 1Figure 34 Dual Reference Track Interferogram Phase102Figure 35 Dual/Single Reference Track Interferogram Phaseis clear that the differential phase is continuous across the discontinuity at range line number512. This result was compared to the interferogram using the new dual/single reference trackapproach yielding a negligible mean difference of 2 mrad or 6 cm in elevation and 6 mradrms or 3 cm rms in elevation across the image. This result shows that the bias errors thatoccur solely for the single reference track case without RCMC are small when the processedbandwidth is smalLFurther interferometric processing of the single reference track interferogram was per103Figure 36 Single Reference Track Interferogram Phaseformed. The differential phase was phase unwrapped and converted into elevation estimateswithout any control points or phase calibration. Figure 37 shows the resulting slant rangetopography with one cycle of the color wheel corresponding to 50 feet. Figure 38 shows the1:50 000 contour map for the region with contour intervals of 25 feet. A visual comparisonof the two demonstrates the effectiveness of SAR interferometry. The spatial distortion wouldbe removed if a projection to ground range were made.1045.2.2 Inertial Data ErrorsIn order to ensure the equality of all approaches with regard to inertial data errors realdata experiments were carried out. A baseline clutter data set was constructed from theThree Hills data which used the INS data supplied by CCRS for motion compensation andelevation estimation. In the test data set the same raw data was processed but with zero mean,independent, Gaussian errors in the INS data. The single and dual track motion compensationapproaches were used as well as the conventional dual track interpretation of the resulting105Figure 37 Elevation Estimates[cz3,ç_Loj1•1.•‘2975’ j fV‘Sc ‘Figure 38 1:50 000 Contour Map of 3 Hills Area (51°45’ lat. 113°30’ long.),25 ft contours (82 P111/12 Energy, Mines, and Resources Canada, 1990)..-------/V5 4 ‘S S ‘/55/l/- 7 ‘I—‘‘—‘5-(S -- ‘-\S 7 (ii i 5,‘5’55 5’ qj—-5-]--S-, \c)c•— -S . (5,51-5 ZiL_Th ‘EE106differential phase. The height estimates were compared to the height estimates from the baselinedata set and the errors computed. As expected, all methods produced the same variance ofheight errors. As an example of the results for the CCRS processing parameters, a 1 mm rmserror in the position of the receive-only antenna produced about a 0.4 m rms height estimationerror at near range.Another similar experiment was conducted, except that the INS errors were constant overthe aperture. As expected, all methods produced the same bias in the height estimates. Forexample, with a 1 mm bias error in the horizontal and vertical position of the receive-onlyantenna position lead to a height estimation bias of about 0.5 m at near range.It would appear that an important limiting factor for airborne InSAR motion compensationis in the ability to estimate the baseline length and angle. For the single and new dualtrack formulations the errors arise in the differential phase to elevation equations. In theconventional dual track approach the instantaneous baseline must be known in order to motioncompensate accurately to the two reference tracks. The resulting baseline induced errors appearas differential phase biases. In either case, the effect on the height estimates are the same.107Chapter 6 ConclusionsIn this chapter the main conclusions drawn from this work will be listed and brieflydescribed. Following this, the ensuing recommended motion compensation approaches forInSAR will be summarized. The following specific conclusions were made from this work:1. The effect of motion compensation with unknown terrain on focussing can be severe forturbulent flights and moderate terrain or typical flights and rugged terrain. The translationalvelocity causes the compressed peak positions to be shifted (Equation 3.48) and translationalacceleration leads to peak broadening (Equation 3.54).2. Due to the range varying phase correction (RVPC) and the ensuing range spectral shiftinduced, the allowed displacement offset between the reference track and the antenna trackis limited by the range oversampling factor if subsequent baseband processing is used(Equation 3.21). This spectral shift may cause problems for interpolation in RCMC orsubsequent range filtering.3. Differential phase bias, and thus a height estimation bias, can occur that dependent uponthe difference in the displacements between each channel and the reference track due tothe range varying phase correction (RVPC). This difference is much larger in the singlereference track case and can lead to significant differential phase biases. This difference istypically very small in the double reference track case. When RCMC is performed the biasis due to differences in quadratic phase errors (Equations 4.82 and 4.83). When RCMC isnot performed, the bias is due to the peak shift caused by uncompensated RCM coupledwith the range phase ramp of the differential RVPC.4. The effect of motion compensation with unknown terrain on the differential phase is smallin most cases due to the phase errors cancelling out between channels. Differential phasebiases can occur if the offset between the antennas and the reference tracks become large,108combined with rugged terrain (Equation 4.31). These effects can be compensated for inthe differential phase to elevation equations. The single reference track approach handlesthis automatically (Equation 4.15) but the double reference track case must include anextra term in the equation (Equation 4.27). Severe angular (roll) acceleration coupled withrugged terrain can lead to differential phase biases (Equations 4.52, 4.54, and 4.60).5. If block azimuth processing is used then the length of the block is limited by the amount ofrange resampling performed over that block as this determines the slant range required forthe matched filter. The important parameter is the difference in the errors in the matchedfilters between the two channels as this leads to differential phase biases (Equations 4.77and 4.79).6. Neglecting to resample for motion compensation does not have much effect for typicalflights but will cause broadening and a reduction of the useful part of the aperture for veryturbulent flights.7. Errors in the inertial data can lead to differential phase errors and biases in the heightestimates. The high frequency errors can be reduced by filtering the interferogram butbias errors will remain. An important factor for airborne InSAR motion compensationis in the ability to estimate the baseline length and angle. For the single and new dualtrack formulations the errors arise in the differential phase to elevation equations. In theconventional dual track approach the instantaneous baseline must be known in order tomotion compensate accurately to the two reference tracks. The resulting baseline inducederrors appear as differential phase biases. In either case, the effect on the height estimatesare the same, requiring extremely accurate inertial data. In the absence of this correction,control points will be required.8. The single reference track approach yields continuous differential phase even when thereference track is segmented (Equation 4.15). In addition, the flat earth fringes have been109removed.9. The dual reference track approach yields discontinuous differential phase when the referencetracks are segmented (Equation 4.27).10. The benefits of both the single and dual reference track approaches can be realized byusing the dual track for motion compensation followed by a conversion to a single track bya subsequent phase correction. If motion compensation resampling is performed then thisextra correction can be very efficiently applied as a constant to the receive-only channelazimuth matched filter for each azimuth line in the segment (Equation 4.110). If motioncompensation resampling is not performed the correction must be applied pixel by pixel tothe entire receive-only image or to the phase smoothed and down sampled interferogram.6.1 Recommended Approach for CCRS SystemFor the typical C-band InSAR scenario for the CCRS system a fairly simple approach tomotion compensation is adequate. This would include:1, No motion compensation resampling.2, Phase correction to a single reference track that may be segmented, if necessary.3. No RCMC.4. Azimuth matched filters determined separately for each channel based on the referencelevel positions.5. Phase unwrapping that assumes continuous differential phase (as far as motion compensation is concerned)6. Differential phase to elevation conversion based on the new single track formulation.Apart from the possible defocussing problems for severe terrain coupled with translationalacceleration, this approach should yield a motion compensation approach which contributes onthe order of a meter to the height estimate error budget.1106.2 Recommended General ApproachA very accurate approach, which might be required for system parameters other than thatof the CCRS case, would be to implement the following:1. Perform motion compensation resampling.2. Phase correction to dual reference track that may be segmented if necessary.3. Perform RCMC,4, Azimuth matched filters determined separately for each channel based on the referencelevel positions. Apply the extra phase constant term to each matched filter for channel Bto convert the channel B reference track to the channel A reference track.5. Phase unwrapping that assumes continuous differential phase (as far as motion compensation is concerned).6. Differential phase to elevation conversion based on the new dual/single track formulation.111Chapter 7 Suggestions for Future WorkThe following is a list of possible extensions to the current work or new ideas that havecome up:1. Remove the restriction in the current work of a yaw steered antenna and zero Dopplerprocessing, Investigate how this changes the sensitivities and recommended processing formotion compensation for InSAR,2. Investigate the role of Autofocus for InSAR motion compensation. For instance,a. Investigate the potential of absolute phase ambiguity determination by manipulatingthe defocussing sensitivity to unknown terrain.b. Use Autofocus to handle defocussing caused by unknown terrain.3. Investigate the modifications to motion compensation for InSAR that might be required ifthe Split-frequency approach [21] is needed.112Bibliography[1] R. 0. Harger. Synthetic Aperture Radar Systems Theory and Design. Academic Press, NewYork, 1970.[2] C. J. Oliver. Synthetic-aperture radar imaging. Journal of Physics D [Applied Physics],22(7), July 1989,[3] John C. Curlander and Robert N. McDonough. Synthetic Aperture Radar Systems andSignal Processing. John Wiley and Sons, Inc., 1991.[4] I. Cumming, D. Hawkins, and A. L. Gray. All-weather Mapping with Interferometric Radar.In 23rd international Symposium on Remote Sensing of the Environment, Bangkok, April1990.[5] F. W. Leberl, J. Raggam, and M. Kobrick. On Stereo Viewing of SAR Images. IEEETransactions on Geoscience and Remote Sensing, 23(2), March 1985.[6] L. C. Graham. Synthetic Interferometer Radar for Topographic Mapping. Proceeding ofiEEE, 62(6), June 1974.[7] H. Zebker and R. Goldstein. Topographic Mapping From Interferometric SAR Observations.Journal of Geophysical Research, 91 (B5) :4993—4999, 1986.[8] A. K. Gabriel and R. M. Goldstein. Crossed Orbit Interferometry: Theory and ExperimentalResults from SIR-B. International Journal ofRemote Sensing, 9(5):857—872, 1988.[9] M. Seymour. Height Estimation by SAR Interferometry Theory and Error Analysis.Technical report, RX-TN-50—3688 MacDonald Dettwiler and Associates, 1992.[10] E. Rodriguez and J. M. Martin. Theory and Design of Interferometric Synthetic ApertureRadars. lEE Proceedings-F, 139(2), April 1992.[11] A. Moccia and S. Vetrella. A Tethered Interferometric Synthetic Aperture Radar (SAR)for a Topographic Mission. IEEE Transactions on Geoscience and Remote Sensing, 30(1),January 1992.[121 H. A. Zebker, M. J. Sander, T. Farr, and R. P. Salazar. Global Topographic Satellite(TOPSAT) Mission. Jet Propulsion Laboratory, April 1993.[13] C. Elachi. Spaceborne Radar Remote Sensing: Applications and Techniques. IEEE Press,New York, 1988.113[14] J. C. Kirk. Motion Compensation for Synthetic Aperture Radar. IEEE Transactions onAerospace and Electronic Systems, 11(3), May 1975.[15] D. Blacknell, A. Freeman, S. Quegan, I. A. Ward, I. P. Finley, C. J. Oliver, R. G. White,and J. W. Wood. Geometric Accuracy in Airborne SAR Images. IEEE Transactions onAerospace and Electronic Systems, 25(2), March 1989.[16] J. R, Moreira, A New Method of Aircraft Motion Error Extraction From Radar Raw DataFor Real Time Motion Compensation. IEEE Transactions on Geoscience and RemoteSensing, 28(4), July 1990.[17] P. Hoogeboom, P. Snoeij, P. J. Koomen, and H. Pouwels. The PHARUS Project, Resultsof the Definition Study Including the SAR Testbed PHARS. IEEE Transactions onGeoscience and Remote Sensing, 30(4), July 1992.[18] A. L. Gray, P. Vachon, K. Mattar, and P. Farris-Manning. Airborne SAR Interferometry atCCRS: Project Status and Height Bias Error Minimization. In CEOS CAL!VAL Workshop,September 1992, Ottawa, Canada, 1992.[19] R. M. Goldstein, H. A. Zebker, and C. L. Werner. Satellite Radar Interferometry: Two-dimensional Phase Unwrapping. Radio Science, 23(4) :713—720, 1988.[20] C. Prati, M. Giani, and N. Leuratti. SAR Interferometry: a 2-D Phase UnwrappingTechnique Based on Phase and Absolute Value Information. In Proceedings of theInternational Geoscience and Remote Sensing Society, 1990.[21] S. Madsen, H. Zebker, and J. Martin. Topographic Mapping Using Radar Interferometry:Processing Techniques. IEEE Transactions on Geoscience and Remote Sensing, 3 1(1),January 1993.[22] H, Zebker, S. Madsen, and J. Martin. The TOPSAR Interferometric Radar MappingInstrument. IEEE Transactions on Geoscience and Remote Sensing, 30(5), September1992.[23] A. L. Gray and P. J. Farris-Manning. Repeat-Pass Interferometry with Airborne Synthetic• Aperture Radar. IEEE Transactions on Geoscience and Remote Sensing, 3 1(1), January1993.[24] S. A. Hovanessian. Introduction to Synthetic Array and Imaging Radars. Technical report,Hughes Aircraft Company, Radar Systems Group, March 1979.[25] A. L. Gray, K. Mattar, and P. Farris-Manning. Airborne SAR Interferometry for TerrainElevation. In Proceedings of the International Geoscience and Remote Sensing Society,Houston, Texas, 1992.114[26] Leland B. Jackson. Digital Filters and Signal Processing. Kiuwer Academic Publishers,1989.[27] R. Bamler and D. Just. Phase Statistics and Decorrelation in SAR Interferograms. InProceedings of the International Geoscience and Remote Sensing Society, Tokyo, Japan,17-21 August 1993.[28] P. L. Bogler. Motion-Compensated SAR Image ISLR. IEEE Transactions on Geoscienceand Remote Sensing, 25(6), November 1987,[29] I. Cumming. The Effect of Matched Filters on Phase Errors. Technical report, AB-RP50—2412 MacDonald Dettwiler and Associates, 1990.[30] P. Beckmann. Probability in Communication Engineering. Harcourt, Brace & World, Inc.,1967.[31] I. S. Gradshteyn and I. M. Ryzhik. Table of Integrals, Series, and Products (Correctedand Enlarged). Academic Press Inc., 1980.[32] T. Sos. ERIM/ARPA/ARMY TEC IFSARE Program. In Proceedings of the DARPAInterferometric SAR Technology and Applications Symposium, 13-14 April 1993.[33] J. C. Dainty. Topics in Applied Physics: Laser Speckle and Related Phenomena. Springer-Verlag, Berlin, 1984.[34] F. Gatelli, A. Monti Guarnieri, F. Parizzi, P. Pasquali, C. Prati, and F. Rocca.Developments in ERS-1 SAR Interferometry: Use of the Spectral Shift. In InternationalWorkshop on SAl? Interferometry, Naples, Italy, 18-20 May 1993.115Appendix A Effect of Quadratic PhaseErrors on CompressionIt is useful to evaluate the effect of quadratic phase errors on the compression of a linear FMpulse. The ideal signal is:s(t) = eu1t2 (A.1)The signal to be considered is:s(t) = e_iK+t2 (A.2)By requiring a large time-bandwidth product (KT2 > 50, from Oliver, 1989) one can obtain anaccurate approximation to the true compression by correlating with the following matched filter:s(t) = rect (;) ejlrKt2 (A.3)where T is the processed aperture time. The compressed peak is therefore:p(t)=J rect()=e_j((KK)(t2+2tT)) dr (A.4)= ,e_t2J e_j2(KK)trejKTdThe remaining integral can be evaluated resulting in Fresnel integrals. As the phase at the peakis of most importance, one can simplify:Ter2d (A.5)This expression is a complex Fresnel Integral or Cornu Spiral. Of interest is the angle of theCornu Spiral. The following substitution can be used to put this into a more recognizable form:(A,6)116yielding,(0) =____f eUdu/irzK J(A.7)2 1_u=e3 du./7rzK J0Now the Cornu Spiral can be evaluated using an asymptotic expansion for small values of x:S(x) = sinu2du- 7.3!(A.8)C(x) = cosu2du— 5.2!The phase error at the peak is:L(x) = tanl ()_(A.9)Back substituting values gives:rLKT2 (A.l0)In terms of the applied quadratic phase error at the edge of the processed band:= KT2 (A.ll)one obtains:3(A.l2)The exact solution was obtained using numerical integration and compared to this approximation. The approximate relation is an upper bound and is very accurate up to an applied errorof about ir at the band edge. In addition, the phase error across the compressed peak wasobserved from the numerical solutions. It was found that the phase was parabolic about thepeak with the error initially decreasing toward the first null, passes through zero, then increasestoward r at the first null. In other words, the expected zero phase response has zero phaseover the main lobe and a change of r radians across each null. When a quadratic phase error ispresent along the aperture the edges of this step-like phase response are smoothed out. About90 percent of the main lobe has a phase error less than that at the peak and the remaining 10percent is around the null.117The effect of quadratic phase errors on the compressed peak 3dB width is very non-linear. Forapplied phase errors at the edge of the band of less than about ir/2 the broadening is less thanabout 5%. For applied errors of about 2.5 radians the broadening is about 20%. Above thisvalue the peak deformation becomes extreme and the broadening increases very rapidly.The effect of quadratic phase errors on the compressed peak sidelobe levels is fairly linear.Above 2.5 radians the sidelobe level has dropped from —13 dB to —5 dB. It would appear thatto maintain a reasonable peak shape the phase errors should be kept below about 7r/2 at theedges of the processed band.For sinusoidal type motions with periods less than the processed aperture, the phase error mustbe less than ir/2 for a maximum allowable broadening of 10% according to Hovanessian, 1979.118Appendix B Single Channel PhaseErrors due to Unknown TerrainGiven that the off-nadir angle to the imaged patch is initially unknown, this leads to phaseerrors in the applied motion compensation phase correction. Referring to Figure 7 one appliesa correction of:qrnc(app1ied) = —(R0— Rao) (B.l)for the zero Doppler pulse of a particular target. The proper correction is:qmc(ideal)— Rao) (B.2)Therefore, the phase mismatch between the data and the assumed data for matched filtering is:/çberr— Rth) (B.3)One can expand this expression for the phase errors to first order to better understand whatparameters are important by using the cosine law and a truncated binomial expansion:I d dRth Rao/l + — 2—cos(c+ Oah)ao ao (B.4)Rao(l— cos(a+Oah))/ d2 dR0 = Raos/l + —— —2——cos(a+ Oao)ao ao (B.5)Rao(1 cos(a+Oao))Therefore, the phase error isZçbr —(cos (C + 0ah) — COS (c + Oao)) (B.6)definingB700ahoand using.1 1cosa — cosb = —2sin(a+ b)sin(a — b) (B.8)yields:err(t)8d(t)sin ((t) + 0(t)+ 0(t))sin (B.9)119The time dependence has been shown to indicate that the phase error for a particular target willvary across the aperture depending upon the flight motion, the applied motion compensation,and the error in the assumed off-nadir angle to the target. It is clear from this expression thatthe phase error is sensitive to flight motion perpendicular to the line-of-sight of the beam. Thephase error variation over the aperture is therefore,8RClJJOS(t)sinLO(t) (B.lO)where djjo8(t) is the flight motion displacement variation perpendicular to the line-of-sight.One can approximate the error in the off-nadir angle with the upper bound:h (B.ll)Rao Sill 0and another upper bound approximation:.Le ZOsrn—— (B.12)Therefore,err(t) 4 hd103(t) (B,13)ARa0 Sifl 0120Appendix C Effect of Linear PhaseErrors on CompressionIt is useful to evaluate the effect of linear phase errors on the compression of a linear FMpulse. The ideal signal is:s(t) = (C.1)The signal to be considered is:.s(t) = e_j7(Kt2+1t) (C,2)By requiring a large time-bandwidth product (KT2 > 50, from Oliver, 1989) one can obtain anaccurate approximation to the true compression by correlating with the following matched filter:s(t) = red (;) ejt2 (C.3)where T is the processed aperture time. The compressed peak is therefore:p(t)= J rect ()ej(K(t)2+L(t+T) ejKr2dr= fej(Kt2+2KtT(t+T))dr (C.4)= e_i(I2t)fe_j2(tTe_jTdr[e_j2tTe_jTdr = — T sin K(2tK + L) (C.5)J ‘,‘r(2tK+L)Therefore, the final solution is:p(t)= _Tsin7r(2tK + L)e_j1r(Kt2+Lt) (C.6)Tr(21K + L)The position of the peak is no longer at time zero but has been shifted to:Ltpeak = —- (sec) (C.7)In addition, there is a phase error at the peak:KL2Oerr =-jj- (rad) (C.8)121Appendix D Directional Random Walk ModelThe complex pixel result from azimuth compression with random phase noise can be describedby a directional random walk model. The statistics of the complex summation result for thecase of communication theory problems has been studied by Beckmann, 1967. The summationcan be characterized by:Re0 = X + jY=(D.l)where each contribution is independent and the expectation of Y is zero (expect zero phase).The A,’ s can be regarded as the modulated scattering strength along azimuth including anycompression windowing. According to Beckinann, the probability density function (PDF) forthe resultant phase is:Ke[_5B2(1+J)1 2= 2 [i + GeG (1 + erf(G))] (D.2)27r(K cos2 0 + sin 0)where___________________G=BKcosO4/ 1+K 2 (D.3)V 2(Kcos2O+siri 0)B=x (D.4)+ <p2>K= I <>V < (x—)2> (D.5)erf(G) = 2fe_Z2dz (D.6)Of interest is the variance of the resultant phase:cr=f62p(O)dO (D.7)A closed form solution to this integral does not exist. But, in most cases of interest anapproximate result is sufficient.When the variance of X and Y are the same (K = 1) one obtains:(0) = e [i + Bcos0(1 + erf(Bcos0))eB2c0s9] (D8)122Furthermore, when the mean of the real part of the resultant is much larger than its variance(B >> 1):B cos 8erf(BcosO) / ez2d 1 (D.9)andcos2O1—O (D.lO)yielding:,,,BcosO—B28e ( . )Due to the very narrow peak of this function one can ignore the cos 0 term. Therefore, thePDF of the resultant phase is approximately Gaussian:B—B28eThe variance of the resultant phase is clearly:02B2 (D.l3)Of interest is to evaluate the phase variance for the case of Gaussian phase noise applied tothe signal prior to compression as a result of high frequency inertial data errors. It is assumedthat the phase noise is zero mean and independent from element to element. The result ofcompression including compression weighting is:N N NRe8 = X + jY=Ae = cos cb + j A sin qj (D.14)The variance to be solved is:2 = <(x X)2 = <x2> = (D.15)X2 X2Furthermore,2 — N (i — N> = e d) A= 2A (D.16)= C1 OSe_d) A = e# A2 (D.17)123The exponential integrals are evaluated using Gradshteyn, 1980, p480. The resulting phasevariance is:N(1_e_20” ZA2‘ i=1 (D.18)2e° (N 2i=1where is the variance of the Gaussian phase noise. When there is uniform weighting thisreduces to:(1 — e_2)2 (DJ9)2N e°124Appendix E Inter-channel Azimuth BroadeningIt is possible to calculate the inter-channel correlation coefficient for distributed homogenoustargets given the point spread function of each channel. The point spread functions will beassumed to be weighted sinc functions that are scaled relative to each other. To evaluate theeffect of inter-channel azimuth broadening one needs to evaluate the following expression:Pai(x)Pa2(x)dx7= (E.1)Fai(x)Pai(x)d) (1Pa2(x)Pa)dx)wheresin 2irWxPai(x) ®w(x) (E.2)27rWxand2irWxPa2(x)=•27rW Øi(x) (E.3)are the azimuth point spread functions for channel 1 and 2, respectively, and where r > 1(® means convolution) and (x) is the inverse Fourier Transform of the frequency domainapplied window. These integrals are evaluated by invoking the following property of FourierTransforms:.T{g(t)}= (f)g(t)dt=(O) (E.4)-Therefore,00 00 2Wf Pai(x)Fa2(x)dx = j {(s c)(s (x)) }dx—00—00Isin2rWx 12irWzØw(x)j = WrectW) x w(f)27rWx (E.6){s1nØü(x)} = rect() x w(f)f Fai(x)Pa2(x)dx = [(rect()w(f)) ® (rect()w(f))]125then,wf Fai(x)Pa2(x)dx = 4W2f w2(f)df00 wf Pai(x)Pai(x)dx = 42 f w2(f)df-00-ww00 77fPa2(x)Padx = 4W2f w2(f)dfResulting in:wJ w(f)dfw-f w2(f)df f w2(f)df-w ji77For example, if there was 10 % inter-channel azimuth broadening this would reduce the inter-channel correlation coefficient by a factor of 0.95.wf w2(f)dff w2(f)df-wIf w(f) = 1 on [-W,W], then1(E.9)(E.10)126Appendix F Single Reference Track Phase ErrorsIn the single reference track case both channels are motion compensated to the same referencetrack (Figure 13). The channel A (transmit/receive) applied phase correction is:qA(applied)—(R0— Rao) (F.1)What should have been applied if the target elevation had been known is:A(ideal) = —(Rth — Rao) (F.2)Therefore, the channel A phase correction error, due to the target height being above thereference level, for this particular range cell and this particular azimuth pulse is:E(A) = kA(applied) — qA(idea1)47r= T(Rto — Rao — (fith — Rao)) (F.3)= -f(Rto--Rth)The channel B (receive-only) applied phase correction is:qB(applied) =—((R0— Rao) + (R0— Rbo)) (F.4)What should have been applied if the target elevation had been known is:cB(ideal) Rth — Rao) + (Rth— Rbh)) (F.5)Therefore, the channel B phase correction error for this range cell and this azimuth pulse is:E(B) = cbB(applied) — qB(ideal)— Bao — Rb0 + Bao + Rbh — 2Rth) (F.6)The difference in the two errors, which is the differential phase, is:= E(A)— E(B) = —(Rbo — Rbh) (F.7)By expanding the slant range values into binomial series one obtains:27rb= —,—(cos (cxab + 0ah) — COS (“ab + Oao)) (F.8)defining00aho (F.9)127and usingcosa—cosb= —2 sin (a+b) sin (a—b) (RiO)yields4Kb (TSifl ab+Oao+jsrn (F.11)4irbjj0A 2But,h (R12)2 2RaosjnOtherefore,2Kb103h (F.i3)—ARa0 Sjfl 0128Appendix G Dual Reference Track Phase ErrorsIn the dual reference track approach each channel is motion compensated to its own track(Figure 12). The applied phase correction in channel A is:A(applied) Rtao— Rao) (G.1)What should have been applied if the target elevation had been known is:qA(ideal) = (Rtah— Bao) (G.2)Therefore, the channel A phase correction error for this range cell and this azimuth pulse is:E(A) = cA(applied)— bA(ideal)47r (G.3)= (Rao— Rtah)The channel B (receive-only) applied phase correction is:B(applied) = Rtao — Rao) + (Rtb0— Rbo)) (GA)What should have been applied if the target elevation had been known is:cbB(zdeal) = — Rao) + (Rtbh— Rbh)) (G.5)Therefore, the channel B phase correction error for this range cell and this azimuth pulse is:E(B) kB(applzed) — qB(ideal) =— ((Rtao — Rao) + (Rtbo — Rbo) — (Rtah — Rao) — (Rtbh — Rbh)) (G.6)In the dual reference track case, one would expect that the differential phase will be the pathlength difference phase, using the reference track positions, plus a component due to the errorsin the applied phase correction because of the unknown target height.The component of differential phase due to pure triangulation is:= Rtah — Rtbh) (G.7)Combining the conventional term with the new term results in:c1 f(Rtah— Rtbh) + E(A) — E(B) (G.8)whereE(A)— E(B) = (Rtao — Rao—Rtah + Baa— Rtb0 + Rb0 + Rtbh — Rbh) (G.9)129which results in:DP = f(Rtao Rtb0 + Rbo Rbh) (G.lO)It is helpful to express the difference in the errors between the channels in a more understandableform. By pairing together terms that can be converted into known parameters and invokingthe parallel ray approximation this yields:Rtao— Rtb0 —bRT (Oo + aRT)Rtbh— Rtah bRT COS (Oh + aRT) (G.l 1)Rb0— Rbh b(cos (Oao + aab) — COS (Oak + cxab))Therefore,E(A)— E(B) (2bRT) (co (o + aRT) — cos (oh + aRT))+(G.12)() (cos (Oah + aab) — COS (Oao +Extending this by defining:(G.13)= °ah — Oaoand usingcos a — cos b = —2 sin (a + b) sin (a — b) (G.14)yields:E(A)-E(B)____I zO\ . (ze’) srn RT + Oao + 2) Slfl 2) — (G.15)/4’7rb’\ / /O(t)’\ /M(t)Sjfl ab(t) + Oao(t) + 2 ) Slfl 2When the displacement from the reference track is small and the antenna baseline is the sameas the reference track baseline (which it is):E(A)-E(B)/LO’\ / / LO’\ / /O’\N (G.16)sin Sifl aRT + °ao + — S1 ab(t) + Oao +using(0.17)R0 sin 0130one obtains:E(A)-E(B)(_4rbh_\ ( (aRT + aab (a — cxabN “s (G.18)RsinO} COS 2 + Oao + 2) SIll 2 ))which gives:(2irhbj03” (0.19)Rsin8)RT — aab)This represents the error recovered by this new approach that is not considered in the conventional dual reference track approach.131Appendix H Unknown Terrain andTranslational MotionIt is helpful to evaluate how translational motion of the aircraft contributes to the differentialphase error. Using the result of Equation 4.75 for the differential phase error along the aperturethe possible effects on the interferogram can be analyzed:2irhbM(t)H 1)R2sin9For linearly varying deviations one would expect to observe a relative mis-positioning betweenthe two channels (Appendix C):irLt 27rhbM(t) (H.2))R smODifferentiating:27rhb OM(t)7TL.\R2 sin o a (H.3)The resulting differential phase at the peak is:irL2 Trh2b 7öM’ 2i—I (H.4)4K 2vAR3sin2 0 \. at IFor any reasonable parameters this differential phase error is negligible.The resulting differential shift in samples is:LFRF hbPRFOM(t)LF =— 2K 2vRsinO a(samples) (H.5)Using 1/16 of a sample as an upper bound we obtain:OM(t) v2RsinOat <8hb FRF (H.6)For typical CCRS parameters (Table 1) this yields:aM(t) 18100 mh (H.7)Clearly, this is not a very limiting factor.If one consider quadratic roll motion (uniform angular acceleration) the peak error along theaperture must be less than 1/3 of our maximum allowed differential phase error (Appendix A):max2?R sinOIn terms of uniform acceleration:12)R S1flO/bjasa< 7rhbT2 (H.9)For typical CCRS parameters (Table 1) this gives:a1742g (H.10)were g is the acceleration due to gravity. Again, it is clear that the differential phase errorsdue to translational motion coupled with unknown terrain elevation is negligible.132Appendix I Motion Compensation ResamplingWhen resampling for motion compensation is neglected there is a residual range migrationerror in each channel that varies along the aperture depending upon the geometry and flightmotion. It is important to explore the nature of the difference in the errors between the twochannels (Figure 39).RTFigure 39 Differential Resampling GeometryAssuming that the point P is the correct sample to follow (flat terrain) the energy of a targetfor flight motion parallel to the reference tracks, the residual migration in channel A is:zR(A)— Rao— BasThe residual migration error in channel B is:The difference between them is:zR(A) = Rb0- Rb8(1.1)(1.2)iR(A) - LR(B) = Rao— Rbo + Rb8 - Ras (1.3)BbARboRts133This can be approximated from:Rbo = RaoJ1 + + 2L COS (aab + 0ao)ao ao (14)Rao (i + COS (aab + Oao))and____________________________Rbs = Ras/1 + + COS (aab + Oas)as as (1.5)Ras (i + COS (nab + Qas))Yielding:/R(A) — zR(B) —(cos (aab + Oao) — COS (aab + Oas)) (1.6)definingOOas17LiO ao °asand usingcos a— cos b = —2 sin (a + b) sin — b) (1.8)yields:LR(A)— zR(B) 2b sin (aab + + Sin (1.9)With further approximations from Figure 11:x(t)cos(t) (110)R(A) — LR(B) bsin(aab + O)x(t)cos(t) (1.11)Provided the roll angle Cab(t) does not vary significantly over the aperture the resulting erroris proportional to the product of the flight deviation and antenna separation perpendicular tothe line-of-sight direction:zR(t) b105Mjjo5(t) (1.12)134Appendix J MATLAB Program ListingsThe following is a sample of the macros used for the MATLAB simulator and InSAR processor.They are written in the MATLAB simulation language.135% simulatorl.mDual channel point target simulation with aircraft motion, motioncompensation to one or two reference tracks, RCMC, azimuth compression,interferogram generation, and point target analysisInitial release, modified from rcmcl.m to add RCMCadded interpolation by 2 of ch A and B before i’gram generationshifted azimuth profile by one sampleadded sinusoidal roll% Initialize Parameters% A/C velocity rn/s% Hz% wavelength (m)% antenna azimuth length Cm)% Range sampling rate Hz% speed of light (m/s)% slant range pixel spacing Cm)% slant range resolution (in)% reference altitude (m)% slant range at nearest approach (m)% from ‘nominal” antenna position (deihlc,delvlc)Ka = 2*V’2/(wave*Rao);% Azimuth FM rate Hz/setac = 0; % azimuth time of nearest approach (sec)Na = 2048; % number of azimuth samples generatedBWproc = 3; % processed bearnwidth (deg.)%Nmf = BWproc*(pi/180)*Rao*PRF/V; % length of azimuth matched filter (saxnpi.)Nmf = 1024; % number of azimuth samples processedK = 50; % azimuth beamwidth constant (3.1 deg one-way)M = pi/(l.125*sires); % range point spread function constantNi = 16; % number of azimuth elements about peak to interpolateNr = 21; % number of range cells (must be odd!)Ir = 16 % interpolation factor for peak analysis% interpolation length% beta used for matched filter kaiser window% avoids divides by zero in sinc(x)% wave number% interpolate in RCMC to i/interp_size of a sample% length of RCMC kernel% Cutoff for RCMC kernel =1 Fr/2% RCMC kernel kaiser beta%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% generate motion compensation datah = 1000; % elevation of target above reference level (m)offset = 0; % offset for extracting azimuth peakBc =1; % distance from antenna a to reference trackthetaBc = 50*pi/l80; % off horizontal angle of antenna a to ref. trackdelhc = Bc*cos(thetaBc); % horizontal offset of channel a from referencedeivc Bc*sin(thetaBc); % vertical offset of channel a from referencedabh = 1.812; % offset between antenna a and b horizontaldabv = 2.141; % offset between antenna a and b verticalb = sqrt(dabh”2÷dabv”2); % antenna baseline distancealphart = atan2(dabh,dabv); % antenna baseline anglebrt = b; % reference track baseline distanceMh = 5; % flight motion horizontal amplitudeMv = 5; % flight motion vertical amplitudeMr = 0.l5*(pi/l80); % roll variation (rad)eta = _(Na_l)/(2*PRF) l/PRF (Na..l)/(2*PRF); % azimuth slow timeprd = 3; % period of sinusoid in seconds%% Change control% June 8/93% June 10/93% June 17/93% Dec 16/93V = 131.0;PRF = 337;wave = 0.0565646;D = 1.0;Fr = 37.5E6;c = 2.9979E8;sipix = c/(2*Fr);sires = 5.6;H = 6000;Rao = 10000;NNi = Ni*Ir;beta = 0;• div_zero = 1E-li;k = 2*pi/wave;interp_size = 16;Nk = 8;Wn = l/interp_size;betai = 3.5;136pphase = 0.25delva = zeros(size(eta));delvb = zeros(size(eta));deiha = zeros(size(eta));delhb = zeros(size(eta));% phase offset of motion in seconds% linear flight motion (+1- M/2 peak)%delva = (Mv*PRF*(Na/Nmf)/(Na_l)).*eta ÷%delha = (Mh*PRF*(Na/Nmf)/(Na_l)).*eta +%delvb = (Mv*PRF*(Na/Nmf)/(Na_l)).*eta +%delhb = (*PRF*(Na/Nmf)/(Na_l)) *eta +% quadratic flight motion ( 0 to +M peak)delvc;delhc;delvc;delhc;%Av = Mv*4*PRF2*(Na/Nmf)2/= 4J.*4*pRF”2*(Na/Nmf)’2/%delva = Av*(eta_pphase)/2%deiha = j*(ta_pphase)A2%delvb = Av*(eta_pphase)/’2%delhb = *(eta_pphase)2(Na—i) ‘2;(Na—i) 2;+ delvc;+ delhc;+ delvc;+ delhc;%delva%delha%delvb%delhb% roll motion% linear roil (+7- Mr/2)+ delvc;+ delhc;+ delvc;+ delhc;%roll = (Mr*PRF*(Na/Nmf)/(Na_l)).*eta + alphart;% quadratic roll (0 to +Mr peak)%Ar = Mr*4*PRF2*(Na/Nmf)’2/(Na_l)2;%roil = Ar*(eta).2 + alphart;% sinusoidal roll (+7- Mr)roll = Mr*sinU2*pi/prd) .*(eta_pphase))% roll+ alphart;delva = delvc*ones(size(eta));deiha = delhc*ones(size(eta));delvb = b*cos(roil) - dabv + delvc;delhb = b*sin(roll) - dabh + delhc;%% mocomp displacements at center of aperture%delvac = delva(Na/2);delhac = delha(Na/2);delvbc = deivb(Na/2);deihbc = delhb(Na/2);% Calculate registered channel b range to h=0 target for middle of aperturedelv = dabv + delvb - delva;delh = dabh + delhb - deiha;alphaab = atan2(delh,delv);alphaabc = alphaab(Na/2);thetaaco = acos((H+delvac)/Rao);Rbo = sqrt(Rao”2 + b’2 ÷2*Rao*b*cos(alphaabC+thetaaco));% Calculate antenna to reference track geometry along aperture% sinusoidal flight motion (+7- M peak)%= Mv*sin( (2*pi/prd) *(eta_pphase)= *sin((2*pi/prd)*(etapphase)= Mv*sin( (2*pi/prd) .*(eta_pphase)= Mj*sin( (2*pi/prd) .*(etapphase)alphala = atan2(delha,deiva)137alphalb = atan2(delhb,delvb);ha = sqrt(delha/’2+delva.”2);has = lla.*lla;llb = sqrt(delhb.2÷delvb.”2);llbs = llb.*llb;%% mocomp displacements%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Generate range varying matched filters for each channelrg = 0;etairf = _(Nmf_l)/(2*PRF) : h/PRF : (Nmf_l)/(2*PRFhfor Ras = Rao_(Nr_l)*slpix/2÷div_zero : sipix : Rao÷(Nr)*slpix/2+div_zerorg = rg + 1;thetal&s = acos((H+delvac)/Ras);Rbs = sqrt(Ras2 + bA2 + 2*Ras*b*cos(thetalas÷alphaabc));Ka = 2*V2/(wave*Ras);mfa = exp(j*pi*Ka*(etac_etamf) .A2_(j*pi*V4*(etac_etaInf) .‘4)/ (wave*2*Ras3))MFA(rg,:) fft(mfa,Na);Kb = 2*V/2/(wave*(Ras+Rbs)/2); Vmfb = exp(j*pi*Kb*(etac_etamf)/.2_(j*pi*VA4*(etac_etaxnf) /‘4)/ (wave*2* C (Ras+Rbs) /2) 3H;MFB(rg,:) = fft(mfb,Na);endV%% Assume that delva,delha,delvb,delhb are measured from individual ref. tks%% calculate deviation of antenna from center of aperture position13a = sqrt((delva-delvac) /‘2+(delha-delhac) /‘2);l3as = 13a.*13a;alpha3a = atan2(delha-delhac,delva-delvac);13b = sqrt((delvb-delvbc) .‘2-i-(delhb-delhbc) /‘2);l3bs = 13b.*l3b;alpha3b = atan2 (delhb-delhbc,delvb-delvbc);% calculate actual range from center of aperture to target for channel bthetaah acos((H-,-delvac-h)/Rao);Rbh = sqrt(RaoS2÷b2+2*Rao*b*cos(alphaabc+thetaah));thetabh = acos( (H+brt*cos(alphart)+delvbc_h) /Rbh);% generate signal data then mocomp it using flat earth for each channelfprintf(’raw data generation in progress\n’);az=0; Vfor eta = _(Na_1)/(2*PRF)+div_zero : h/PRF : (Na_l)/(2*PRF)+div_zeroaz = az + 1;if (rem(az,64) == 0)fprintf(’hine %d \n’,az);end% calculate Ra in range-elevation plane from antenna a to target: assuming% that range to target from ‘nominal” track is always Ras independent% of h (this is for simplicity of equations for changing h)Rla = sqrt(Rao”2-i-l3as(az) +2*Rao*13a(az)*cos(thetaah+alpha3a(az)));% calculate Rb in range-elevation plane from antenna b to target VRlb = sqrt(Rbh”2-t-l3bs(az) +2*Rbh*13b(az)*cos(thetabh+alpha3b(az)));% calculate instantaneous range from antenna a to target(RCM included)138Ra = sqrt(RlaA2+V2*(eta_etac)s2);Aa = (sin(K*V*eta/Ra)/(K*V*eta/Ra))/2; % antenna a beam pattern% calculate instantaneous range from antenna b to target(RCM included)Rb = sqrt(Rlb2+V’s2*(eta_etac)2);Ab = (sin(K*V*eta/Rb)/(K*V*eta/Rb))2; % antenna b beam pattern% Generate range compressed data in rangerg = 0;for Ras = Rao_(Nr_l)*slpix/2+div_zero slpix : Rao÷(Nr)*slpix/2+div_zerorg = rg + 1;thetalas = acos((H+delva(az))/Ras);Rbs sqrt(Ras’2 + b’2 + 2*Ras*b*cos(thetalas+alphaab(az)));thetaibs = acos((H+brt*cos(alphart)+delvb(az))/Rbs);dataa(rg,az) = Aa*exp(_j*2*k*Ra) * sin(M*(Ras_Ra))!(M*(Ras_Ra));datab(rg,az) = Ab*exp(_j*k*(Ra÷Rb))*sin(M*((Ra÷Rb)/2_(Ras÷Rbs)/2))/(M*((Ra÷Rb)/2_(Ras+Rbs)/2));Calculate required motion phase correction, dual trackRtao = sqrt(Ras”2 + llas(az)- 2*Ras*lla(az)*cos(thetalastalphala(az)));deira Rtao - Ras;Rtbo = sqrt(Rbs2 + llbs(az)- 2*Rbs*llb(az)*cos(thetalbs+alphalb(az)));delra = Rtao - Ras;delrb = Rtbo - Rbs;now mocomp it to reference paths using flat earth assuniptiondata_a(rg,az) = data_a(rg,az)*exp(_j*k*2*delra);data_b(rg,az) = data_b(rg,az)*exp(_j*k*(delra+delrb));endendfprintf(’raw data generation complete\n’);%% Calculate expected differential phase from center of aperture positionsthetalas = acos((H+delvac)/Rao);thetaibs = acos( (H+brt*cos(alphart)+delvbc) /Rbo);Rtao = sqrt(Rao’2 + llas(Na/2)- 2*Rao*lla(Na/2)*cos(thetalas+alphala(Na/2H);Rtbo = sqrt(Rbo’2 + llbs(Na/2)- 2*Rbo*llb(Na/2)*cos(thetalbs÷alphalb(Na!2)));DPP = k*(Rtbo - Rtao + Rbh - Rbo);DPe = atan2(sin(DPP),cos(DPP));% generate interpolation kernel for RCMCkernel = interp_size*firl(interp_size*Nk,Wn,kaiser(interp_size*Nk+l,betaiH;% extract every 8th value to generate 16 different filters% - one for interpolating to the ith/l6 sample intervalfor 1=1: interp_size;for m=l:Nk;filt(l,m) = kernel((m_l)*interp_size÷l+l);endendtt = l:Na;cdata = zeros(size(data_a));fftdata = zeros(size(tt));cfftdata = zeros(size(tt));Fa = 0 PRF/Na : (Na_l)*PRF/Na;phaser = exp(j*2*pi*(Nmf/(2*PRF)).*Fa);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% perform RCMC on channel A%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%139for rg = l:Nr;fftdata(rg, :) = fft (data_a(rg, :));endfor az = l:Na;if (rem(az,64) == 0)fprintf(’line %d \n’,az);endif (az <= Na/2+l)shift = (2*Fr/c)*0.5*V2*((az_l)*PRF/Na)S2/(KaS2*Rao); % in sampleselseshift = (2*Fr/c)*O.5*V/S2*((Na÷1_az)*PRF/Na)S2/(Ka2*Rao); % in samplesend% determine integer shift part%mt shift = floor(shift + l/(2*interp_size));% determine fractional shiftfrs = round(interp_size*(shift — mt_shift));% perform integer shift on datafftrcmcdata(l:Nr-jnt_shift,az) = fftdata(int_shift+l:Nr,az);% perform fractional interpolation and discard Nk/2 points at start% leaving them zero and leave last Nk/2+l points zero, note also% that the last mt_shift non-zero values are throwaway%filt_no = interp_size - frs;for rg l:Nr-Nk-i-l;fftrcmcdata(rg+Nk/2-l,az) = sum(fftrcmcdata(rg:rg+Nk-l,az) .filt(filt_no, :) •‘);endend% compress channel Afor rg = l:Nr;MF MFA(rg,:);cfftdata = fftrcrncdata(rg, :) *phaser;cdata_a(rg,:) ifft(cfftdata);end% perform RCMC on channel Bfor rg = l:Nr; Vfftdata(rg, :) = f ft (data_b(rg, :));endfor az = l:Na;if (rem(az,64) == 0)fprintf(’line %d \n’,az);endif (az <= Na/2+l)shift = (2*Fr/c)*0.5*V2*((az_l)*PRF/Na)’2/(Ka2*Rao); % in sampleselseshift = (2*Fr/c)*0.5*V/2*((Na+l_az)*PRF/Na)F2/(Ka2*Rao); % in samplesend% determine integer shift partmt_shift = floor(shift + l/(2*interp_size));% determine fractional shift140ifs = round(interp_size*(shift- mt_shift));perform integer shift on datafftrcmcdata(l:Nr-int_shift,az) = fftdata(int_shift÷l:Nr,az);% perform fractional interpolation and discard Nk/2 points at start% leaving them zero and leave last Nk12+l points zero, note also% that the last mt_shift non-zero values are throwawayfilt_no = interp_size - frs;for rg = l:Nr-Nk+l;fftrcmcdata(rg+Nk/2-l,az) = suxn(fftrcmcdata(rg:rg-i-Nk-1,az) .filt(filt_no, :) •‘);endend% compress channel Bfor rg l:Nr;MF = MFB(rg,:);cfftdata = fftrcmcdata(rg, :) .*MF.*phaser;cdata_b(rg,:) = ifft(cfftdata);end.fprintf ( ‘data compressed\n’)% interpolate channels A and B by 2 to prepare for i’gram generationfor rg = l:Nr;cdata_ar(rg, :) =interp(cdata_a(rg,offset+Na!2-(Ni-l)/2:Na/2+(Ni-l)/2+offset),2);cdata_br(rg, :)interp(cdata_b(rg,offset+Na/2—(Ni-l)/2:Na/2+(Ni-l)/2+offset),2);endfor az l:2*Ni;cdata_arr(:,az) = interp(cdata_ar(:,az),2);cdata_brr(:,az) = interp(cdata_br(:,az),2);endinter = cdata_arr . conj(cdata_brr);cd = inter(Nr/2:3*Nr/2_l,Ni/2:3*Ni/21);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% perform point target analysis of interferogram%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% interpolate in rangecdf = fft(cd); % fft columns (each of length Nr)Nri = Ir*Nr;cdfi = zeros(Nri,Ni);cdfi(l:(Nr—l)/2÷l,:) = cdf(l:(Nr—l)12+l,:);cdfi(Nri—(Nr—l)/2÷l:Nri,:) = cdf((Nr—l)/2+2:Nr,:);cdi = ifft(cdfi);% interpolate in azimuth%cdff = fft(cdi.’); % fft columns (each of length Ni)Nai = Ir*Ni;cdfii zeros(Nai,Nri);cdfii(l:Ni/2+l, :) = cdff(l:Ni/2+l, :);cdfii(Nai—Ni/2+2:Nai, :) = cdff(Ni/2+2:Ni, :)cdii = ifft(cdfii).’;% find peak and plot results141[Val,Row] = max(abs(cdii));[val,midxc] = max(Val);midxr = Row(midxc);mid = Nai/2;[w3a, sia) = widiobe(cdii (midxr, :));[w3r, sir) = widiobe(cdii ( ,miclxc));peakphase = angle (cdii (midxr,midxc));% calculate estimated elevationabs_phase = peakphase + ((DPP_rem(DPP,2*pi))/(2*pi)_l)*2*pi;h_est = H+delvac_Rao*cos (-alphaabc÷acos ( ( ((wave*abs_.phase/ (2*pi) )÷Rtao-Rtbo+Rbo)2 - b2_Rao’2)/(2*b*Rao)));n = i:Nai;axis([1 Nai —val val]);plot(n,abs(cdii(midxr,:)),n,angle(cdii(midxr,:))*.3183*val);title(’Azimuth Cross—section through point target’);xlabel C ‘Azimuth samples’);ylabel(’Magnitude and phase’);text (mid, .9*val, sprintf(’Error in estimated elevation %2.3f m’ ,h-h_est));text (mid, .8*val,sprintf(IMeasured Peak phase %2.3f rad.’, peakphase))text(mid, .6*val,sprintf(IExpected Peak phase %2.3f rad.’, DPe))text(mid, .4*val,sprintf(lPeak phase error %2.Of mrad.’,1000* (peakphase-DPe)))text(mid, .2*val,sprintf(fAz. 3dB width %2.3f m’, V*w3a/(PRF*2*Ir)))text(mid,0,sprintf(’Az. First side lobe %2.2f dB’, sla))text(mid,_.2*va1,sprintf(Rg. 3dB width %2.2f m’, slpix*w3r/(2*Ir)))text(mid,_.4*val,sprirltf(IRg. First side lobe %2.2f dB’, sir))text(mid,_.E*val,sprintf(fPeak at %2.0f/%2.Of range’,midxr,Nri));text(mid,_.8*val,sprintf(Peak at %2.0f/%2.Of azimuth’,midxc,Nai));grid%meta plotpausen = l:Nri;mid = Nri/2;axis([l Nri -val val]);plot(n,abs(cdii(:,midxc)),n,angle(cdii(:,midxc))*.3l83*val);title( ‘Range Cross-section through point target’);xlabel(’Range samples’);ylabel(’Magnitude and phase’);text (mid, .8*val, sprintf(’Measured Peak phase %2.3f rad.’, peakphase))text(mid, .6*val,sprjntf(IExpected Peak phase %2.3f rad.’, DPe))text (mid, .4*val,sprir1tf(Peak phase error %2.Of mrad.’,1000* (peakphase-DPeH)text(mid,.2*val,sprintf(Az. 3dB width %2.3f m’, V*w3a/(PRF*2*Ir)))text(mid,0,sprintf(’Az. First side lobe %2.2f dB’, sla))text(mid,_.2*val,sprintf(IRg. 3dB width %2.2f m’, slpix*w3r/(2*IrH)text(mid,_.4*val,sprintf(Rg. First side lobe %2.2f dB’, sir))text(mid,_.6*val,sprintf(Peak at %2.Of/%2.Of range’,midxr,Nri));text(mid, _.8*val, sprintf(’Peak at %2.Of/%2 .Of azimuth’ ,midxc,Nai));gridpause%% extract peak area and do 2D interpolation of channel A%cd = cdata_a(l:Nr,offset+Na/2—(Ni-1)/2-l:Na/2+(Ni—l)/2÷offset—l);% interpolate in rangecdf = fft(cd); % fft columns (each of length Nr)Nri = Ir*Nr;cdfi = zeros(Nri,Ni);cdfi(1: (Nr—l)/2+l, :) = cdf(l: (Nr—l)12÷1, :)cdfi(Nri—(Nr—l)/2-t-l:Nri,:) = cdf((Nr—l)/2+2:Nr,:);cdi = ifft(cdfj);142% interpolate in azimuth%cdff = fft(cdi.’); % fft columns (each of length Ni)Nai = Ir*Ni;cdfii = zeros(Nai,Nri);cdfii(l:Ni/2+l, :) = cdff(l:Ni/2+l, :)cdfii(Nai—Ni/2+2:Nai, :) = cdff(Ni/2+2:Ni, :);cdii = ifft(cdfii).’;% find peak and plot results[Val,Row] = xnax(abs(cdii));[val,midxc] = max(Val);midxr = Row(midxc);mid = Nai/2;[w3a, sia] = widiobe (cdii (midxr, :));[w3r, slr) = widiobe(cdii( : ,midxc));peakphase = angle (cdii (midxr,midxc));Pae = angle (exp(_j*4*pi*Rtao/wave));n = l:Nai;axis([i Nai -val vail);plot(n,abs(cdii(midxr, :)),n,angle(cdii(midxr, :))*.3183*val);title(’Azimuth Cross-section through point target Ch A’);xlabel ( ‘Azimuth samples’);ylabel(’Magnitude and phase’);text (mid, .8*val, sprintf(’Measured Peak phase %2.3f rad.’, peakphase))text (mid, .6*va1,sprintf(Expected Peak phase %2.3f rad.’, Pae))text (mid, .4*va1,sprintf(Peak phase error %2.0f mrad.’,1000* (peakphase-PaeH)text(mid,.2*val,sprintf(IAz. 3dB width %2.3f m’, V*w3a/(PRF*Ir)))text(mid,0,sprintf(’Az. First side lobe %2.2f dB’, sla))text(mid,_.2*val,sprintf(fRg. 3dB width %2.2f m’, slpix*w3r/Ir))text(mid,_.4*val,sprintf(?Rg. First side lobe %2.2f dB’, sir))text(mid,_.6*val,sprintf(Peak at %2.0f/%2.0f range’,midxr,Nri));text(mid,_.8*val,sprintf(Peak at %2.Of/%2.Of azimuth’ ,midxc,Nai));grid%meta plotpausen = l:Nri;mid = Nri/2;axis([i Nri —val val]);plot(n,abs(cdii(: ,midxc)) ,n,angie(cdii(: ,midxc))*.3183*val)title(’Range Cross-section through point target Ch A’);xlabei ( ‘Range samples’);ylabel(’Magnitude and phase’);text(mid, .8*val, sprintf(’Measured Peak phase %2.3f rad.’, peakphase))text (mid, .6*val,sprintf(Expected Peak phase %2.3f rad.’, Pae))text (mid, .4*val, sprintf(’Peak phase error %2.Of mrad.’,1000* (peakphase-PaeH)text(mid,.2*val,sprintf(Az. 3dB width %2.3f m’, V*w3a/(PRF*Ir)))text(mid,0,sprintf(’Az. First side lobe %2.2f dB’, sla))text(mid,_.2*val,sprintf(Rg. 3dB width %2.2f m’, slpix*w3r/Ir))text(mid,_.4*val,sprintf(sRg. First side lobe %2.2f dB’, slr))text(mid,_.6*val,sprintf(Peak at %2.Of/%2.Of range’,midxr,Nri));text(mid,_.8*val,sprintf(Peak at %2.Of/%2.Of azimuth’ ,midxc,Nai));gridpause% extract peak area and do 2D interpolation of channel B%cd = cdata_b(l:Nr,offset-i-Na/2-(Ni-l)/2-l:Na/2+(Ni-1)/2+offset—l);% interpolate in rangecdf = fft(cd); % f ft columns (each of length Nr)Nri = Ir*Nr;cdfi = zeros(Nri,Ni); 143cdfi(l: (Nr—l)/2+l, :) = cdf(l: (Nr—l)/2÷l, :)cdfi(Nri—(Nr--l)/2+l:Nri, :) = cdf((Nr—l)/2+2:Nr, :);cdi = ifft(cdfi);% interpolate in azimuth%cdff = fft(cdi.’); % fft columns (each of length Ni)Nai = Ir*Ni;cdfii = zeros(Nai,Nri);cdfii(l:Ni/2+l,:) = cdff(l:Ni/2+l,:);cdfii(Nai—Ni/2+2:Nai, :) = cdff(Ni/2+2:Ni, :);cdii = ifft(cdfii).’;% find peak and plot results[Val,Ràw) = rnax(abs(cdii));[val,midxc] = max(Val);midxr = Row(midxc);mid = Nai!2;[w3a, sla) = widlobe(cdii (midxr, :));[w3r,slr) = widlobe(cdii(:,midxc));peakphase = angle(cdii(midxr,midxc));Pbe = angle(exp(_j*2*pi*(Rtao+Rtbo)/wave));n = l:Nai;axis([l Nai -val val]);%psfl000 = abs (cdii (midxr, :)%save psfl000lO24.mat psfl000plot(n,abs(cdii(midxr, :)),n,angle(cdii(midxr, :))*.3183*val);title(’Azimuth Cross-section through point target Ch B’);xlabel (‘Azimuth samples’);ylabel(’Magnitude and phase’);text(mid, .8*val,sprintf(Measured Peak phase %2.3f rad.’, peakphase))text(mid, .6*val,sprintf(fExpected Peak phase %2.3f rad.’, Pbe))text (mid, .4*val,sprintf(lPeak phase error %2.Of mrad.’,1000* (peakphase-Pbe)))text (mid, .2*val,sprintf(Az. 3dB width %2.3.f m’, V*w3a/(PRF*Ir)))text(mid,0,sprintf(’Az. First side lobe %2.2f dB’, sla))text(mid,_.2*val,sprintf(IRg. 3dB width %2.2f m’, slpix*w3r!Ir))text(mid,_.4*val,sprintf(IRg. First side lobe %2.2f dB’, sir))text(mid,_.6*val,sprintf(Peak at %2.Of/%2.0f range’,midxr,Nri));text(mid,_.8*vaisprintf(Peak at %2.Of/%2.Of azimuth’ ,midxc,Nai));grid%rneta plotpausen = l:Nri;mid = Nri/2;axis([l Nri -val val]);plot (n,abs (cdii(: ,midxc)) ,n,angle(cdii(: ,midxc))*.3l83*val)title(’Range Cross-section through point target Ch B’);xlabel ( ‘Range samples’);ylabel(’Magnitude and phase’);text (mid, .8*val,sprintf(Measured Peak phase %2.3f rad.’, peakphase))text(mid, .6*val,sprintf(Expected Peak phase %2.3f rad.’, Pbe))text(mid, .4*va1,sprintf(Peak phase error %2.Of mrad.’,1000* (peakphase-Pbe)))text(mid, .2*val,sprintf(Az. 3dB width %2.3f m’, V*w3a/(PRF*Ir)))text(mid,0,sprintf(’Az. First side lobe %2.2f dB’, sla))text(mid,_.2*va1,sprintf(Rg. 3dB width %2.2f m’, slpix*w3r/Ir))text(mid,_.4*val,sprintf(Rg. First side lobe %2.2f dB’, sir))text(xnid,_.6*val,sprintf(Peak at %2.Of/%2.Of range’ ,midxr,Nri));text(mid,_.8*vai,sprintf(Peak at %2.Of/%2.Of azimuth’ ,midxc,Nai));grid144% widlobe.m% 1 D point target analysis[width,slobe] = WIDLOBE (input_array)- input_array must be amplitude units- Calculate 3dB width, first side lobe levelNote: Peak and sidelobes must not be at an edge% change control% feb 92 - initial release% nov 12/92 - fixed up index checkingfunction [widt.h,slobe] = widlobe (x)% Convert to power unitspow = 20*loglo(abs(x));[maxval,maxidx) = max(pow);pow = pow - rnaxval;% locate 3dB width (to sample spacing accuracy only)%k = maxidx-l;threshold = -3.0;1 = length(x);while ((k>l) & (pow(k) > threshold))k = k - 1;endif (pow(k) > threshold)sprintf(’could not find left 3dB point’);endleft = k + (-3.0 - pow(k))/(pow(k÷1)-pow(k));k = maxidx-t-l;while ((k<1) & (pow(k) > threshold))k = k + 1;endif (pow(k) > threshold)sprintf(’could not find right 3dB point’);endright = k - (-3.0 - pow(k))/(pow(k-l)-pow(k));width = right - left;% Locate first sidelobek = maxidx—l;while ((k > 2) & (pow(k-l) < pow(k)))k = k - 1;endif (pow(k-l) < pow(k))sprintf(’could not find left sidelobe’);endleftlobe = max(pow(l:k-l));k = maxidx+l;while ( (k < (1-1)) & (pow(k+l) < pow(k)))k=k+l;endif (pow(k÷l) < pow(k))sprintf(’could not find right sidelobe’);endrightiobe = max(pow(k÷l:length(pow)));slobe = max(leftlobe,rightlobe);145% regl.m% Register Channel B to Channel A at reference level% Change control% Jan 18/94 Initial release% Jan 20/94 modified for multi-segments% Initialize Parametersatodfreq=37 .5;delay=[13.289, 13.2737];wave=0 . 05 65 646c=299.79;delsr=c/ (2*atodfreq);rgd=5l.99;rbase=(rgd-delay(1) )*c/2;H=6059;K=2 *pj /wave;D = 1;% generate interpolation kernelkernel = interp_size*firl(interp_size*Nk,Wn,kaiser(interp_size*Nk÷1,betai));% extract every 8th value to generate 16 different filters- one for interpolating to the ith/16 sample intervalfor 1=1: interp_size;for m=l:Nk;filt(l,m) = kernel( (m—l) *interpsize+l+l);endendload kpl0prf % load in motion dataddelh = delh(:,2)—delh(:,l);ddelv delv(:,2)-delv(:,l);b = sqrt(ddelh.”2 + ddelv.”2);alphaab = atan2(ddelh,ddelv);n_range= 12000; # of range lines no. available to processsrange= 1; % first range line no. to processerange=srarige÷n_rarige-l; % last range line no. to processn_azimuth= 180; % # of azimuth lines to processpstart= 1; % starting azimuth linepend = pstart+n_azimuth-l;% last azimuth lineuniend=(segmerit-1) *(nazimuth_lO)+n_azimuth; % unimat last pixel #unistart= uniend—ri_azimuth + 1; % unimat file start range pixel no.% Load in the dataload plObsk=pstart:perid;Ra=(delsr* (k-t-uriistart-2) )+rbase;% Perform the resampling% Range sampling rate (MHz)% SAW delay (usec) A,B% wavelengthspeed of light in m/microsec% slant range pixel spacing (4 m)% range gate delay in microsec.% base range, only ch A used% elevation of reference track% wavenumber% antenna length (in)segment = 6; % segment numberinterp_size = 16; % interpolate to l/interp_size of a sampleNk = 8; % length of kernelWn = l/interp_size;% Cutoff for kernel =1 - Fr/2betai = 3.5; % kernel kaiser beta146rdata = zeros(size(data));for l=srange : erangetemp=(H-f-delv(l, 1)) . IRa;temp2=find(temp>l);temp(temp2)=ones(size(temp2));theta_a=acos (temp);Rb = sqrt(Ra/’2 + b(lY’2 + 2*Ra*b(l) •* (theta a+alphaab(l)));chB pixel # - chA pixel # = dein for coincident scene patchesdeln = _(Ra_Rb_c*(delay(2)_delay(l)))./(2*delsr);intshift = floor(deln + lI(2*interp_size)); % integer sample shiftfrs = round(interp_size*(deln - mt_shift)); % fractional sample shiftfprintf(’ Processing range line %4.Of\n’, 1);% process integer and fractional shiftsindex = 1;for k=pstart:pendistart = k + int_shift(index) - Nk/2 + 1;iend = istart + Nk - 1;if ((istart >= pstart) & (iend <= pend))filt_no = iriterp_size - frs(index);rdata(k,l) = sum(data(istart:iend,l).*filt(filt_no, :).‘);elserdata(k,l) = 0;endendenddata = rdata;save plobrs data147% igraml.m% Single channel motion compensation and azimuth compressión% Change control% Jan 17/94 Initial Release% Jan 18/94 modified for dual reference tracks% Jan 20/94 corrected transpose problem, modified for segment processing% Jan 21/94 modified for dual track% Initialize Parametersatodfreq=37 .5;delay=[13 .289, 13.2737];wave=0 . 05 65 646c=299.79;delsr=c/ (2*atodfreq);rgd=51.99;rbaset(rgd-delay(1) )*c/2;H=6059;K=2 *pj /wave;D = 1;%segment 6;channelsw=2;% Range sampling rate (MHz)% SAW delay (usec)% wavelengthspeed of light in m/microsec% slant range pixel spacing (4 m)% base range, only ch A used% elevation of reference track% wavenuxnber% antenna length (m)% image segment number across range% processing channel A (1) or B (2)n_range= 12000; % # of range lines no. available to processsrange= 1; % first range line no. to processerange=srange+n_range-l; % last range line no. to processn_azirnuth= 180; % # of azimuth lines to processpstart= 1; % starting azimuth lineuniend=(segment_l)*(rl_azimuth_lO)÷n_azimuth; % unimat last pixel#unistart= uniend-n_azimuth + 1; % unimat file start range pixel no.swathend = 2048; % last range pixel no. in swath to process% load in the motion compensation file and extract relevant dataload kpl0prf % mocomp file name for single reference trackdelha=delh(srange:erange, 1);delva=delv(srange:erange, 1);brt=2.8013; % reference tk baseline length (0 for single track)alphart=40.22*(pi/180); % reference track baseline angle% transform insar channel to its own reference trackdelhb=delh(srange:erange, 2)_brt*sin(alphart);delvb=delv(srange:erange, 2) _brt*cos (alphart);delv = brt*cos(alphart) + delvb - delva;delh = brt*sin(alphart) + delhb - delha;b = sqrt(delv.”2+delh.’2); % antennaalphaab = atan2(delh,delv); % antennaalphala = atan2(delha,delva); % channelalphalb = atan2(delhb,delvb); % channelha = sqrt(delha.”2+delva.”2);ilas = lla.*lla;llb = sqrt(delhb/’2+delvb/’2);llbs = llb.*llb;% calculate azimuth processing parametersbaseline lengthbaseline angleA angle to reference trackB angle to reference track% Obtain processing instructions148mspacing=mean (spacing (srange erange));maxr=rbase÷delsr*(swathend_l); % maximum range, channel Aproc_phi = 2.23*(pi/180); % azimuth beamwidth to processnf ft = 13824; % forward FFT lengthmaxrefsize = proc_phi*maxr/mspacing; % maximum matched filter length% load in registered raw dataif ( channelsw == 1 );load plOas;elseload plObrs;end% phase correct the raw data based on the navigation dataindex = 1;for k=pstart :pstart+n_azimuth-1;Ra=(delsr*(k+unistart_2))+rbase;temp= (H÷delva) IRa;temp2=find(temp>l);temp(temp2)=ones(size(temp2H;theta_a=acos (ternp);Rtao = sqrt(Ra”2 + has - 2*Ra*lla.*cos(theta_a+alphala));delra = Rtao - Ra;if channelsw==l;delr=2*delra;elseRb = sqrt(Ra”2 + b.2 + 2*Ra.*b.*cos(theta_a+alphaab));temp= (H+delvb+brt*cos (alphart)) /Rb;temp2=find(temp>l);temp(temp2)=ones(size(temp2));theta_b=acos (temp);Rtbo = sqrt(Rb/’2 + llbs - 2*Rb.*llb.*cos(theta_b+alphalb));delrb = Rtbo - Rb;deir = delra+delrb;Rbmf(index) = Rb(n_rangel2);end % of if channelsw ==lfprintf(’ Phase correcting azimuth line %4.Of, range=%8.lf\n’, k,Ra);phc=exp(_sqrt(_l)*K*delr);data(k,l:n_range)=data(k,l:n_rarige).*phc(1:n_range).;index = index + 1;end; % end of for k= loop% Loop to azimuth process the dataouti = zeros (n_azimuth, n_range-maxrefsize);index=1;for k=pstart :pstart+n_azimuth-l;if channelsw==1;r=(delsr* (k+unistart-2) )+rbase;elser=Rbmf(index);endfprintf(’ Processing azimuth line %4.Of, range=%8.lf\n’, k,r);temp=fft(data(k,1:n_range), nfft);ref=refgen(r, mspacing, wave, proc_phi);nref=length(ref);temp=ifft(fft(ref, nfft) .*temp);startindex= (nref÷maxrefsize) /2;endindex=n_range- (maxrefsize-nref) /2;outl(index, 1 :endindex-startindex+1)=temp(startindex:endindex);inc1ex=index--1;149endif channelsw==1;save sega.xnat outlelsesave segb.mat outiend%%%%% Generate azimuth matched filter %%%%%%%%%%%%%%%%%%%%%function [ref]=refgen(r, spacing, wave, Proc_phi)Nref=l+2*round(r*proc_phi/2/spacing);for j=l:Nrefref(j)=exp(i*2*pi*(spacing)A2*(1_Nref/2)A2/r/wave);end150% elevl.m% Estimate terrain height from interferogram% Change control% Feb 1/94 Initial release% Initialize Parametersatodfreq=37.5; % Range sampling rate (MHz)delay[l3.289, 13.2737]; % SAW delay (usec) A,Bwave=0.0565646; % wavelengthc=299.79; % speed of light in m/microsecdelsr=c/(2*atodfreq); % slant range pixel spacing (4 m)rgd=5l.99; % range gate delay in microsec.rbase=(rgd_delay(l))*c/2; % base range, only ch A usedH=6059; % elevation of reference trackK=2*pi/wave; % wavenurnberD = 1; % antenna length (m)% load in igramload cigrams%% load in the motion compensation file and extract relevant data%load kploprf % mocomp file name for single reference trackn_range= 1024; % # of range lines no. available to processn_azimuth= 2045; % # of azimuth lines to processdelha=delh(l:10:10*n_range,l);delva=delv(l :10: 10*n_range, 1);%brt=2.8013; % reference track baseline lengthbrt=0; % reference track baseline lengthn = 1:512;alphart = zeros(size(delha));alphart(l:5l2)=40.22*(pi/l80)*ones(size(n)); % reference track baseline anglealphart(5l3:l024)=40.8*(pi/180)*ones(size(n)); % reference tk baseline angle% transform insar channel to its own reference trackdelhb=delh(l:lO:10*n_range,2)_brt.*sin(alphart);delvb=delv(l:l0:l0*n_range,2)brt.*cos(alphart);delv = brt*cos.(alphart) + delvb - delva;delh = brt*sin(alphart) + delhb - delha;b = sqrt(delv/’2÷delh.2); % antenna baseline lengthalphaab = atan2(delh,delv); % antenna baseline anglealphala = atan2(delha,delva); % channel A angle to reference trackalphalb = atan2(delhb,delvb); % channel B angle to reference tracklla = sqrt(delha.”2+delva/’2);has = lla.*lla;llb = sqrt(delhb.2+delvb.2);llbs = llb.*llb;% perform further phase smoothingfiltnr = 6; % filter size rangefiltna = 6; % filter size azimuthfor i=l:filtnrfor j=l:filtnafilt(i,j) = l/(filtnr*filtna);end 151endcigramsf = conv2(cigrams,filt, ‘same’);cigrarnsf = rot90(cigramsf);a = angle(cigrarnsf);for i=1:80a(:,i) = zeros(size(lla));end%% phase unwrap assuming no residues!aa(l,:) =unwrap(a(1,:));for i=1:n_azimuthfor j=2:n_rangediff = a(j,i) — a(j—l,i);if (abs(diff) < pi)aa(j,i) = aa(j—l,i) + diff;elseif (diff < 0)aa(j,i) = aa(j—l,i) + 2*pi + diff;elseaa(j,i) = aa(j—l,i) — 2*pi + diff;endendendend.% Calculate motion compensation data and height estimate%h=zeros (size (aa));for k=l:n_azimuthRa=(delsr*(k_l))+rbase;temp= (H+delva) IRa;temp2=find(temp>1);temp(temp2)=ones(size(temp2));theta_a=acos (temp);Rtao = sqrt(Ra2 + has - 2*Ra*hla.*cos(theta_a+alphala));Rbo = sqrt(RaA2 + b.2 + 2*Ra.*b.*cos(theta_a+alphaab));temp=(H+delvb+brt*cos(alphart)). /Rbo;temp2=find(temp>l);temp(temp2)=ones(size(temp2));theta_b=acos(temp);Rtbo = sqrt(Rbo/’2 + libs - 2*Rbo.*llb.*cos(theta_b÷alphalb));calculate height estimatetempO = (wave* (aa(:,k)_2*pi)/(2*pi));temp00 (tempO +Rbo)/’2;templ = acosUtemp0O - b/2_RaA2)../(2*b.*RaH;h(: ,k) = (H÷delva_Ra.*cos(_alphaab+ tempi));fprintf(’processed azimuth line %d\n’,k);endsave h.mat hh_disp=rem( (64*3.28/100) *(h+4o) ,64);index = 0;for i=1:l6:l024temp(1:16,1:80) = index*ones(16,80);h_disp(i:i+15,l:80) = temp(1:16,l:80)index = index + 1;endimage (h_disp)title(’3 Hills, Alberta: Elevation’)xlabel(’Slant Range Samples’)ylabel ( ‘Azimuth Samples’)152

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065119/manifest

Comment

Related Items