Geomatica-insar-user-guide.pdf

InSAR User Manual Geomatica 2017

April 2017

Revision History Version

Date

Description

Author

1.0

February 2017

Initial Document

John Wessels

1.1

April 2017

Minor Edits to DEMADJUST and INSCALDEFO)

John Wessels

Executive Summary: This InSAR user manual is intended to provide the reader with an overview of the processing steps required to generate interferometric products using PCI’s Geomatica™ package. The current InSAR package can be used to generate topographic products to characterize digital surface models (DSMs) or deformation products which identify subtle changes in surface elevation due to land subsidence. Although it is possible to generate very high quality DSMs, our emphasis is placed upon the generation of multiple layers of deformation products which can be easily co-registered into temporally ordered stacks for long term analysis of subsidence patterns. Our objective is to provide a commercial, end-to-end, interferometric processing, temporal analysis, and visualization capability designed to significantly reduce the processing requirements and operational constraints encountered by non-SAR experts. In general, the user can accept the default values to automate the processing steps. Each module is designed to be flexible and intuitive. The developed modules can be executed either independently or sequentially via an EASI™ or Python™ script. The extracted motion and/or displacement estimates can be reformatted to virtually any map projection supported by Geomatica™. The DSM and subsidence output can be combined with optical, vector or cadastral data to enhance the spatial context of the extracted information. The targeted application areas include monitoring surface changes due to mineral and water extraction, oil and gas extraction, glacial rebound, volcanic activity and landslides. The output products can also be used to monitor the integrity of urban infrastructure (e.g. bridges), transportation networks (e.g. road and rail) and pipelines. The current InSAR package supports single-look, complex valued, slant range data acquired in strip map mode from single or multi-polarized configurations. End-to-end support is available for the RADARSAT-2 and TerraSAR-X sensors. The ability to generate interferograms is currently available for TANDEM-X, Cosmo-SkyMed, PALSAR-2, and Kompsat-5 and end-to-end support for these sensors is anticipated in future releases. The workflow sequence can be summarized as: 

Creation of Interferograms o SAR data ingest, extraction of metadata and geocoding information o Automated image to image co-registration o Generation of raw or filtered interferogram



Generation of Interferometric Products o Removal of flat earth and topographic effects o Removal of residual orbital errors o Unwrapping the interferometric phase o Adjustment of Interferometric Output Products o Convert to Orthorectified Final Product



Temporal Analysis of Surface Subsidence o Temporal parameter extraction from unwrapped interferometric stacks o Enhanced visualization of temporal information.

This document provides a summary of each InSAR module; identifies and provides a brief description of the required input and output parameters; and gives an example for each processing step. Each module is presented in the sequence in which it would normally be executed. Input parameters may be entered using the EASI™ or Python interface.

Acknowledgments: This project was partially funded by the Canadian Space Agency under the Earth Observation Application Development Program (EOADP) contract number 9F043-130644/006/MTB entitled “Application Development for Environmental Monitoring and Remediation.” The author wishes to acknowledge the significant contributions toward the successful completion of this initiative by a number of individuals. In particular, Pierre Nadeau of the Canadian Space Agency; Naomi Short and Sergey Samsonov of the Canada Centre for Mapping and Earth Observation; as well as Guylain Greer, Juliusz Ostrowski, Gabriel Gosselin, and Michael Horlings of PCI Geomatics.

Table of Contents Revision History ........................................................................................................... i Executive Summary: ................................................................................................... ii Acknowledgments: .................................................................................................... iv Overview .................................................................................................................... 1 Quick Start for Deformation Products ......................................................................... 3 Quick Start for Topographic Products .......................................................................... 5 Automation of Processing Modules ............................................................................. 7 Supported Sensors and Data Types ............................................................................. 8 Helpful Hints ............................................................................................................. 10 Selecting Reference and Dependent Data Sets ......................................................................10 Selecting a DEM ...........................................................................................................................10 Assessing the Quality of Interferograms ...................................................................................11

Interferometric Products........................................................................................... 13 Digital Terrain Models ..................................................................................................................13 Subsidence Maps .........................................................................................................................14 Ancillary Interferometric Products ..............................................................................................16 Temporal Analysis Products .......................................................................................................17

Troubleshooting ....................................................................................................... 18 Creation of an Interferogram .................................................................................... 20 INSINFO (Analysis of Interferometric Datasets) ......................................................................21 SARINGEST (Data Ingest, Extraction of Metadata and Geocoding)....................................21 INSCOREG (Automated Image to Image Coregistration) ......................................................22 INSRAW (Generation of Raw or Adaptively Filtered Interferogram) ....................................23 Example of How to Create an Interferogram ............................................................................23 Example of How to Display an Interferogram using FOCUS .................................................29

Topographic and Deformation Products .................................................................... 31 Overview ........................................................................................................................................31 What You Need to Know .............................................................................................................32 INSTOPO (Removal of Flat Earth and Topographic Phase) .................................................33 INSADJUST (Removal of Residual Orbital Errors) .................................................................35 INSUNWRAP (Unwrapping the Interferometric Phase) ..........................................................38 How to Unwrap a Topographic Interferogram to Generate a DSM .......................................40 How to Unwrap a Deformation Interferogram (Subsidence Map) .........................................42

DEMADJUST (Adjustment of Topographic Products) ............................................................44 Steps to Calibrate Unwrapped Topographic Products ............................................................45 GEOCODEDEM (Convert DSM to Geocoded Product) .........................................................47 INSCALDEFO (Adjustment of Deformation Products) ............................................................48

Temporal Analysis of Deformation Products.............................................................. 52 ORTHO (Orthorectification of Deformation Products) ............................................................52 INSSTACK (Temporal Stacking of Deformation Products) ....................................................54 An example illustrating how the INSSTACK output is calculated .........................................55 Temporal Analysis of Interferometric Stacks ............................................................................59

Appendix A ............................................................................................................... 61 A.1 InSAR Data Analysis.............................................................................................................61 A.2 InSAR Data Ingest.................................................................................................................63 A.3 Coregistration of Reference and Dependent Images ......................................................65 A.4 Generation of Raw (or Filtered) Interferogram ..................................................................69

Appendix B ............................................................................................................... 71 B.1 Flat Earth and Topographic Correction ..............................................................................71 B.2 Orbital Correction and Residual Fringe Removal .............................................................74 B.3 Phase Unwrapping ................................................................................................................77 B.4 Adjustment of Topographic Products .................................................................................82 B.5 Conversion from Slant Range to Orthorectified Product .................................................84 B.6 Adjustment of Deformation Products ..................................................................................86

Appendix C ............................................................................................................... 88 C.1 Create a Stack of Displacements or Velocities .................................................................88 C.2 Temporal Analysis Tool ........................................................................................................91

Complex Value Display in FOCUS ............................................................................... 92

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Overview Radar interferometry is a technique which combines the coherent phase from a pair of synthetic aperture radar (SAR) signals to provide three-dimensional information about the Earth’s surface. The technique requires the coregistration and processing of two SAR signals from very similar viewing geometries which have been acquired over the region of interest. By coherently combining the signals from the two acquisitions, the interferometric phase between the received signals can be determined for each imaged point. The interferometric phase represents the difference in the path lengths to the imaged point. The interferometric phase is affected by topography, atmospheric effects, surface motion and system noise. The effects of topography can be estimated and removed through the use of a (high quality) digital elevation model (DEM) and knowledge of the radar characteristics (wavelength, range to target) and relative orbital positions (baseline length and drift). Once the effects of topographic, atmospheric and residual noise have been removed, the remaining phase differences are due to subsidence in the direction of the radar’s line of sight (LOS). If multiple data sets are available, a long term temporal assessment of motion can be created for each target. Our objective is to provide an end-to-end processing system designed for non-SAR experts capable of fully exploiting the rich information content that can be extracted from interferometric SAR data sets in an user friendly and intuitive way. The software modules can process both topographic and deformation information derived from multiple SAR sensors using data sets ranging from a single interferometric pair to multiple stacks of single, dual, compact and fully polarimetric data to produce reliable and valuable map products in a timely manner. The intuitive and flexible interface significantly reduces the complexity and redundancy of the processing steps. The processing sequence can be summarized as 

Creation of Interferograms o SAR data ingest and extraction of metadata and geocoding information o Automated image to image co-registration o Generation of raw or filtered interferogram



Generation of Interferometric Products o Removal of flat earth and topographic effects o Removal of residual orbital errors o Unwrapping the interferometric phase o Adjustment of Interferometric Output Products o Convert to Orthorectified Final Product



Temporal Analysis of Surface Subsidence o Temporal parameter extraction from unwrapped interferometric stacks o Enhanced visualization of temporal information.

Version 1.1

April 2017

1

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

The modules described in this document can be executed either separately or sequentially and currently support single-look, slant range data acquired in strip map mode from single and multipol configurations. A brief description of the steps required for interferometric processing is presented in the nominal order of InSAR processing. For additional information, the user is referred to the appendices where each module is identified by name, followed by an indication of whether the module is considered to be optional or mandatory and a brief description of its purpose. The required input and output parameters are described and an example of the expected output is provided where applicable. Required parameters and default values (if applicable) are indicated. The configuration showing the nominal InSAR workflow sequence and the interconnectivity of the modules is shown in Figure 1. Mandatory Process Vendor SAR SLC Optional Process Optional Ancillary Information

Data Ingest (SARINGEST)

Deformation Products Topographic Products

Report (INSARINFO)

Image Coregistration (INSCOREG)

Raw/Filtered Interferogram (INSRAW)

Flat Earth and Topo Removal (INSTOPO)

DEM

Orbital Refinement (INSADJUST) Phase Unwrapping (INSUNWRAP)

Adjustment Points

Vector Mask

Calibrated DSM (DEMADJUST)

Calibrated Subsidence Map (INSCALDEFO)

Orthorectified DSM (GEOCODEDEM)

Orthorectification (ORTHO)

Adjustment Points

Line of Sight Measurement (INSSTACK)

Figure 1: Overview of Interferometric Workflow Temporal Analysis

Version 1.1

April 2017

2

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Quick Start for Deformation Products What you Need Before You Start Vendor Supplied Reference File Vendor Supplied Dependent File Digital Elevation Model

Step 1:

X Y Z

Create an Interferogram Process

Ingest Reference File Ingest Dependent File Coregister Reference and Dependent File Create a Raw Interferogram

PPF

Input

Output

SARINGEST SARINGEST INSCOREG INSRAW

X Y A,B A,C

A B C D

Recommended Parameter Settings for Interferogram Generation CALIBTYP = “Sigma” FLSZ = 7

! Recommended calibration type ! Recommended filter size

Example of How to Create an Interferogram

FILI = Reference Dataset FILO = Reference.pix Run SARINGEST

! This is the original reference file “X” ! This will be the reference file in PCIDSK format “A” ! Read in the reference file

FILI = Dependent Dataset FILO = Dependent.pix Run SARINGEST

! This is the original dependent file “Y” ! This will be the dependent file in PCIDSK format “B” ! Read in the dependent file

FILREF = Reference.pix FILI = Dependent.pix FILO = Resample.pix Run INSCOREG

! This is the PCIDSK reference file “A” ! This is the PCIDSK dependent file “B” ! This will be the name of the coregistered dependent file “C” ! Co-register the dependent file to the reference file

FILREF = Reference.pix FILI = Resample.pix FILO = Raw_Interferogram.pix Run INSRAW

! This is the PCIDSK reference file “A” ! This is the name of the coregistered dependent file “C” ! This is the file name to be given to the raw interferogram “D” ! Create the interferogram

Step 2:

Create a Calibrated Deformation Product Process

PPF

Create A Wrapped Deformation Product Fix Orbit Effects for Wrapped Deformation Product Unwrap Phase for Deformation Map Adjust for Random Deformation Offset

Version 1.1

April 2017

INSTOPO INSADJUST INSUNWRAP INSCALDEFO

Input

Output

D,Z E F G

E F G G

3

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Example of How to Create a (Single) Channel Deformation Product

FILI = Raw_Interferogram.pix FILEDEM = DEM.tif FILO = Raw_deformation.pix Run INSTOPO

! This is the raw interferogram “D” ! This is the digital elevation model “Z” (used to remove topography) ! This will be the wrapped deformation file “E” ! Correct the raw interferogram for flat earth and topographic effects

FILI = Raw_Deformation.pix DBIC_OCS =1 FILO =Orbit_Adjusted_Deformation.pix Run INSADJUST

! This is the wrapped deformation file “E” ! This is the channel number to use to calculate the orbital correction ! This will be the deformation file adjusted for orbital effects “F” ! Automatically remove any residual phase due to orbital error

FILI = Orbit_Adjusted_Deformation.pix DBIC = 1 FILO = Unwrapped_Deformation.pix Run INUNWRAP

! This is the phase wrapped orbit adjusted file “F” ! This is the channel to unwrap (can only do 1 at a time) ! This will be the name of the unwrapped deformation “G” ! Unwrap the subsidence layer

FILE = Unwrapped Deformation.pix DBIC = 2 DBOC = 3 TFILE = Adjustment_Points.txt Run INSCALDEFO

! This is the unwrapped deformation file “G” ! The deformation (in meters) with random offset is on channel 2 ! We want to write the calibrated output to channel 3 ! This is text file containing the image locations of adjustment points ! Calibrate the deformation product

Step 3:

Perform Temporal Analysis of Surface Deformation Process

PPF

Temporal Merging of Deformation Products Temporal Analysis and Visualization

INSSTACK FOCUS (Charts)

Input

Output

G1…Gn H

H 

Example of How to Create a Temporally Ordered Stack for Analysis

MFILE = List_Of_Defo_Products.txt FILO = Deformation_Stack.pix OQUANT = “DISP” UNWRAP = “NO” DUNITS = “cm” TUNITS = “year” Run INSSTACK

! this is a text file of all the deformation products (G1......Gn) ! this is the name to be given to the deformation stack “H” ! select displacement as the output parameter ! select NO because we already calibrated the output ! we want to measure cumulative displacement in centimeters ! we want to set the time axis (for analysis by FOCUS ) in years

To perform temporal analysis  Load Deformation_Stack.pix “H” into FOCUS  AnalysisInSAR Temporal Chart will show temporal displacement in (centimeters) with time axis in years

Version 1.1

April 2017

4

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Quick Start for Topographic Products What you Need Before You Start Vendor Supplied Reference File Vendor Supplied Dependent File Coarse DEM (for optional adjustment)

Step 1:

X Y Z

Create an Interferogram Process

Ingest Reference File Ingest Dependent File Coregister Reference and Dependent File Create a Raw Interferogram

PPF

Input

Output

SARINGEST SARINGEST INSCOREG INSRAW

X Y A,B A,C

A B C D

Recommended Parameter Settings for Interferogram Generation CALIBTYP = “Sigma” FLSZ = 7

! recommended calibration type ! recommended filter size

Example of How to Create an Interferogram

FILI = Reference Dataset FILO = Reference.pix Run SARINGEST

! This is the original reference file “X” ! This will be the reference file in PCIDSK format “A” ! Read in the reference data set

FILI = Dependent Dataset FILO = Dependent.pix Run SARINGEST

! This is the original dependent file “Y” ! This will be the dependent file in PCIDSK format “B” ! Read in the dependent data set

FILREF = Reference.pix FILI = Dependent.pix FILO = Resample.pix Run INSCOREG

! This is the PCIDSK reference file “A” ! This is the PCIDSK dependent file “B” ! This will be the name of the coregistered dependent file “C”

FILREF = Reference.pix FILI = Resample.pix FILO = Raw_Interferogram.pix Run INSRAW

! This is the PCIDSK reference file “A” ! This is the name of the coregistered dependent file “C” ! This is the file name to be given to the raw interferogram “D” ! Generate the raw interferogram

Version 1.1

April 2017

5

PCI Geomatics Enterprises Inc.

Step 2:

InSAR User Guide Geomatica 2017

Create a Calibrated Topographic Product Process

PPF

Create A Wrapped Topographic Product Fix Orbit Effects for Wrapped Topographic Product Unwrap Phase for Topographic Map Adjust for Random Topographic Offset Convert DSM from Slant Range to Orthorectified

INSTOPO INSADJUST INSUNWRAP DEMADJUST GEOCODEDEM

Input

Output

D J K L L

J K L L M

Example of How to Create a Topographic Product FILI = Raw_Interferogram.pix FILO = Raw_Topographic.pix Run INSTOPO

! This is the raw interferogram “D” ! This will be the wrapped topographic file “J” ! Generate an interferogram of topographic fringes

FILI = Raw_Topographic.pix DBIC_OCS =1 FILO =Orbit_Adjusted_Topographic.pix Run INSADJUST

! This is the wrapped topographic file “J” ! This is the channel number to use to calculate the orbital correction ! This will be the topographic file “K” adjusted for orbital effects ! Remove residual phase due to orbital errors

FILI = Orbit_Adjusted_Topographic.pix DBIC = 1 FILO = Unwrapped_Topographic.pix Run INUNWRAP

! This is the phase wrapped orbit adjusted file “K” ! This is the channel to unwrap (can only do 1 at a time) ! This will be the name of the unwrapped DEM “L” ! Unwrap DEM (output contains relative elevations)

LASC = 0 NUMFORM = “GEOCOORD,Suppress” FILE = Coarse_DEM.tif DBIC = 1 SAMPLING = 20, 20 TFILE = Elevation_Points.txt Run NUMWRIT

! Clear the last segment created value ! Suppress the NO DATA Areas ! Coarse DEM to be used for adjustment ! Channel containing elevations ! Example of DEM sampling interval (in image coordinates) ! This will be the text file containing location of the adjustment points ! Generate the adjustment Text File

FILV = Elevation_Points.txt VECUNIT = “Long/Lat” FILE = Unwrapped_Topographic.pix Run VREAD

! This is the text file with the adjustment points ! Specifies the projection of adjustment points from text file ! Adjustment points will be added to unwrapped topographic product ! Copy the vector control points to the DEM

FILE = Unwrapped_Topographic.pix PCIOP = “Add” PCIVAL = 0,0,0,1 Run PCIMOD

! This is the name of the file to add another channel ! We are going to add a adjustment channel ! The new adjustment channel is real valued ! Add an adjusted channel

FILV = Unwrapped_Topographic.pix FILEDEM = Unwrapped_Topographic.pix DBVS = LASC DBIC = 2 DBOC = 3 TFILE = Adjustment_Points.txt Run DEMADJUST

! The control points are in the unwrapped file “L” ! This is the unwrapped deformation file “L” ! The vector segment is the last one created ! The topography (in meters) with random offset is on channel 2 ! We want to write the calibrated output to channel 3 ! The text file containing the image locations of adjustment points ! Adjust the DEM to absolute value

MFILE = Unwrapped_Topographic.pix DBEC = 3 MAPUNITS = “UTM” FILO = Topographic_Final.pix Run GEOCODEDEM

! File to be orthorectify “L” ! Calibrated elevation is on channel 3 ! output to be written in UTM ! Name of final orthorectified DEM Produce an orthorectified DSM

” Version 1.1

April 2017

6

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Automation of Processing Modules Many of the classic processing algorithms developed for interferometry are based upon trigonometric relationships. In theory, trigonometric techniques work well. However; for satellite based interferometry where the satellite ranges to the target are always orders of magnitude larger than the baseline distance between the satellites; even very small angular errors (i.e. hundredths of a degree) can lead to a severe degradation of the interferometric product. Until recently, the use of the classical model was sufficient because the accuracy of the orbital information provided by commercial vendors was also limited. The lack of precision of the orbital data often made it necessary to process interferometric datasets multiple times by iteratively adjusting the satellite positions and velocities until a satisfactory result was achieved. Subsequently there was no real advantage to improving the classic model since interferometric processing was essentially a trial and error process anyway. Newer satellites incorporate GPS technologies, have much better orbital control and provide a much higher level of precision in the orbital metadata. Since interferometric processing is essentially a three dimensional problem requiring precise knowledge of the satellite positions, velocities, and target range, we reformulated the classic model to incorporate vector geometry which is a better representation of the four dimensional (t , x, y, z ) aspect of interferometric processing. We developed a system capable of automatically extracting the improved metadata and reprojecting the pertinent information into a standardized, four dimensional coordinate system. The vector based mathematics is used to accurately determine all of the required interferometric parameters such as the interferometric gradient derived from the 3D baseline and target location calculated as the look vector intersection with the earth model. Using vector based calculations for each pixel location minimizes the accumulation of errors because each calculation is independent of the calculation for any other pixel. The accuracy of the final result is still constrained by the accuracy of the metadata. A major advantage of vector based analysis is the ability to determine the “direction” of the error. We exploit this knowledge to automatically adjust the relative satellite positions until we converge to an optimum solution. This ability to automatically correct the satellite’s position significantly reduces the complexity of the processing since manual adjustments and subsequent reprocessing are no longer required. The benefits derived from the use of vector based algorithms (i.e. improved precision and ability to self-correct) make it much easier for non-SAR experts to process interferometric data.

Version 1.1

April 2017

7

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Supported Sensors and Data Types The 2017 InSAR package supports single look complex (SLC) data sets acquired in strip map mode. Support for single and multi-pol data is provided. Full support is available for:  

RADARSAT-2 TerraSAR-X

The 2017 InSAR package also provides support for the generation of interferograms for:    

PALSAR-2, Cosmo-SkyMed, TanDEM-X Kompsat-5

Full support for these sensors will be available in subsequent releases.

Strip map is acquired with the radar antenna pointing in a fixed direction (usually perpendicular to the flight direction). Each radar pulse generates a scan line in the range direction. The beam mode (i.e. swath width and incident angle range configuration) is selectable There are no spatial gaps in the azimuth direction. The strip map acquisition concept is shown in Figure 2.

Figure 2: Strip Map Mode

Version 1.1

April 2017

8

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Support for spotlight mode is not currently available. Spotlight mode achieves high resolution by steering the radar beam over areas of interest but some areas are not imaged. The spotlight acquisition concept is shown in Figure 3.

Figure 3: Spotlight Mode

Support for ScanSAR mode or Sentinel-1 TOPS mode is not currently available. ScanSAR and TOPS cover a much larger swath area by illuminating multiple overlapping swaths at different ranges. The improved coverage comes at the expense of resolution. The ScanSAR acquisition concept is shown in Figure 4.

Figure 4: ScanSAR Mode

Version 1.1

April 2017

9

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Helpful Hints Selecting Reference and Dependent Data Sets Prior to selecting your data sets, it is important to know what type of interferometric product you wish to produce. In general, it is a good idea to use data sets with the smallest time difference between acquisitions. This will minimize the effects of temporal decorrelation. If possible, you should avoid acquisitions that occur during or shortly after significant weather events such as rain or snow falls. If you plan to generate a topographic product such as a DEM, you should select data pairs with a perpendicular baseline in the range of hundreds of meters. If the perpendicular baseline is too small (i.e. tens of meters) the DEM becomes overly sensitive to phase noise and the results are not reliable. Conversely if the perpendicular baseline is too large (i.e. more than a third of the critical distance), it becomes difficult to unwrap the phase due to geometric decorrelation. If you plan to generate subsidence (deformation) products you should select data sets with the smallest possible perpendicular baseline. Small baselines improve the overall coherence and reduce the effect of positional or elevation errors in the DEM used to remove topographic phase. In general, perpendicular baselines in the tens of meters should be used. If you plan to generate a series of layers showing cumulative displacements over a long period of time, you should try to keep the temporal spacing between the acquisitions fairly evenly spaced. This is because the cumulative displacements are derived by integrating the estimated velocities over the time interval represented by the interferogram. Any error in the velocity estimate integrated over a long time interval will bias the estimate of the cumulative displacement.

Selecting a DEM Digital elevation models are used to determine the phase adjustment necessary to remove the interferometric phase due to elevation. If the ELEVATION_DATUM metadata tag is set, the appropriate elevation model is used. If the elevation model is not specified, the DEM is assumed to be in meters above mean sea level (MSL). The accuracy of the topographic phase removal step is directly related to the accuracy of the elevations extracted from the DEM. For typical radars with a slant range of 1000 km. and an incident angle of 35 degrees, the error in surface displacement is about 1 mm. per 40 meter error in elevation. Since we would like to measure displacement to millimeter accuracy, the vertical accuracy of the DEM should be better than  20 m. To avoid elevation errors due to interpolation, the spacing of the DEM should be within a factor of 10 spacing of the SAR data. Digital elevation models from the SRTM mission can be

Version 1.1

April 2017

10

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

freely downloaded at http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1 or at https://earthexplorer.usgs.gov/

Assessing the Quality of Interferograms To generate an interferogram it is essential for both SAR signals to have been acquired using identical wavelengths and polarizations. In addition, the radar viewing geometry must be identical using datasets with identical pass direction, look direction, and beam mode. However; confirming the radar viewing geometry and acquisition parameters are valid does not necessarily ensure a good result since the quality of the interferogram will be degraded if the radar targets have changed in any way or moved more than a quarter of a wavelength between acquisitions. In general, some degradation in quality is expected if the time between acquisitions is very large (i.e. years). A reduction in interferometric quality can also occur in the short term due to the effects of severe weather (e.g. windstorms, flooding, snowfall etc.) or changes in the target scattering (e.g. vegetation growth, harvested fields, random movement etc.). For areas with a large amount of vegetation, the quality of the interferogram may be improved by using a satellite with a longer wavelength. Poor interferometric quality can also occur in areas with low signal levels due to low backscattering (e.g. water) or in regions of radar shadow. The quality of the final interferometric product (DEM or subsidence map) is directly related to the quality of the raw interferogram. Prior to generating the final product, it is prudent to visually examine the raw interferogram. The Geomatica™ visualization tool FOCUS, provides the ability to readily assess the quality of an interferogram by examining the coherence (i.e. amplitude). High coherence areas show up as “bright” areas in the display. Figure 5 shows the coherence for a given reference file with three different dependent data sets. The difference in acquisition times is 48 days, 96 days and 456 days respectively. The 96 day file should not be used because of the low coherence. There was a significant snowstorm which occurred between the first and second acquisition which explains the reduction in coherence.

48 days

96 days

456 days

Figure 5: Comparison of Various Coherence Levels

Version 1.1

April 2017

11

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

The quality of the interferogram can also be assessed by examining the phase pattern. To show the both the coherence and interferometric phase, select one the “Amplitude and Phase” viewing option. Areas with random phase (i.e. no discernable fringe pattern) denote regions of low reliability while areas with clear fringes indicate regions where high quality results can be expected. Figure 6 shows an interferogram displayed in FOCUS as “Amplitude and Phase” which clearly shows high coherence regions (bare surface) and low coherence regions (due to rapidly moving glacier)

Figure 6: Interferogram showing areas of high and low coherence

Version 1.1

April 2017

12

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Interferometric Products Digital Terrain Models If the orbital position and velocity of the SAR sensors is precisely known, the interferometric phase between the two signals can be used to estimate topography. For the generation of digital terrain models, the optimum distance between satellite positions (baseline distance) is in the range of hundreds of meters (typically 100 to 600 m.). A relative terrain model can be generated from a coregistered pair of coherent SAR data sets by removing the interferometric fringes due to the known orbital geometry (automatically derived from the metadata). The product can be further enhanced by providing a few ground control points to derive a low order polynomial correction which can be applied to remove any residual atmospheric effects and correct the random offset to produce the absolute elevation. Figure 7 shows an example of a digital terrain model generated by the package

0

Elevation in Meters

700

Figure 7: Example of Digital Terrain Model Product

Version 1.1

April 2017

13

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Subsidence Maps If the local topography is known, the effects of elevation on the interferometric phase can be estimated and removed through the use of a high quality DEM. The quality of the DEM used for correction should be as high as possible. In general, the spacing of the elevations should not exceed ten times the spacing of the interferogram. The vertical accuracy should be within ten meters. The remaining interferometric phase represents subtle changes in the earth’s elevation (i.e. subsidence or deformation). Since the time between the acquisitions is also known, these changes in surface elevation can be represented as either a user defined distance measurement representing subsidence as shown in Figure 8 or as a velocity measurement as shown in Figure 9.

.

-2.0

Cumulative Subsidence in Centimeters

2.0

Figure 8: Example of Cumulative Subsidence

-25

Subsidence Velocity in Centimeters per Year

Figure 9: Example of Subsidence Velocity

Subsidence products can be used in a number of applications to precisely measure any changes in the surface elevation due to water/oil/gas/mineral extraction, volcanic activity, melting permafrost or glacial rebound. These products can also be used to monitor subtle shifts in urban infrastructure (e.g. buildings, bridges etc.), road and rail networks, and pipelines. Spatial registration and temporal ordering of multiple subsidence products provides an accurate

Version 1.1

April 2017

14

25

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

and reliable mechanism for the long term monitoring of surface displacement and/or surface velocity over the area of interest. An illustration of a temporal stack is shown in Figure 10.

Figure 10: Multiple Deformation Layers Organized as a Temporal Stack

An example illustrating the output from the temporal analysis tool is shown in Figure 11 where the X axis is time in years and the Y axis shows cumulative displacement in centimeters.

Figure 11: Temporal Analysis Tool output showing Cumulative Displacement at Cursor Position Version 1.1

April 2017

15

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Ancillary Interferometric Products The generation of a number of primary and ancillary interferometric products is also supported. These include the following: 

Calibrated SAR Data (to  0 ,  0 ,  0 )

       

Coherence (including support for multiple polarizations) Raw Interferogram (wrapped) Topographic and/or Deformation Phase (wrapped) Topographic and/or Deformation Phase (unwrapped) Orbital Phase Adjustment (removal or orbital fringes) Calibrated Topographic Product(s) Calibrated Subsidence Map(s) Atmospheric Phase Screen

Examples of ancillary interferometric products are shown in Figure 12

Calibrated Intensity

Coherence

Raw Interferogram

Topographic Phase

Calibrated Digital Elevation Model

Orbital Phase Adjustments

Atmospheric Phase Screen

Calibrated Subsidence Map

Figure 12: Examples of Interferometric Products

Version 1.1

April 2017

16

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Temporal Analysis Products The InSAR package also supports the generation of temporal products shown in Figure 13. These include:  Layers of Cumulative Displacement (in user defined units)  Velocity Estimates at temporal midpoints (in user defined units)  Average Linear Velocity  Average Coherence

Cumulative Displacement

Best Linear Velocity

Average Coherence

Figure 13: Examples of Temporal Analysis Products

Version 1.1

April 2017

17

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Troubleshooting The quality of the interferometric process is highly dependent upon the quality of the two signals. Areas that have undergone any changes between acquisitions will have low coherence. Areas with very weak radar returns such as water bodies will also result in low coherence. As much as possible, we have tried to anticipate potential processing issues and mitigate the problem by passing the required metadata and/or providing default values. The best way to increase the likelihood of getting a good interferometric result is:        

Minimize the time between acquisitions for the reference and dependent data sets Use data sets with precise orbital information (usually available 2 days after acquisition) Use high resolution data sets Use co-polarized data sets (i.e. HH, VV) which have a stronger signal rather than cross polarized (i.e. HV) data sets which are generally weaker Avoid very steep (less than 25°) or shallow (more than 50°) incident angles For the generation of DEMs use data pairs with perpendicular baselines in the range of 100 – 800 meters For the generation of subsidence maps use data pairs with perpendicular baselines in the range of 0 - 100 meters Use the highest quality DEM available for topographic phase removal and/or adjustment

A list of potential problems and recommendations follows: 1. Cannot coregister the data sets (INSCOREG)  Increase the search radius of the image matcher  Increase the number of points to use  Decrease the required minimum score 2. Not getting a high quality interferogram (INSRAW)  Use data sets with reduced time between acquisitions  To minimize effects of uncompensated orbital drift, reduce the size of the area to be processed 3. Interferometric artefacts remain after removal or topographic phase (INSTOPO)  Ensure the quality of the DEM is a high as possible  Ensure the quality of the orbital information is as high as possible 4. Phase unwrapping takes a long time (INSUNWRAP)  Mask out the low coherence areas  Process smaller areas independently  Increase the size of available computer memory

Version 1.1

April 2017

18

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

5. Tiling Artefacts appear in the unwrapped data (INSUNWRAP)  Process smaller areas  Increase the size of available computer memory 6. Calibrated Output has systematic error (DEMADJUST)  Use a higher density of adjustment points by changing the sampling interval

Version 1.1

April 2017

19

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Creation of an Interferogram The first requirement is to create an interferogram. The generation of an interferogram requires the coregistration of multiple SAR data sets; acquired with similar radar characteristics and viewing geometries over a non-changing and relatively stationary area of interest. These restrictions are required to ensure that the coherence which is an estimate of the quality of the interferogram remains high. Coherence is reduced when the target changes significantly between satellite acquisitions. For example this might happen if one of the data acquisitions occurs after a significant snowfall or flooding. Coherence is also reduced if the target is not stationary. As a rule of thumb, targets which move more than a quarter of a wavelength will have a low coherence. Therefore the first task is to ensure that the SAR data sets can be used to generate a high quality interferogram. To address the issue of whether the data sets were acquired using similar radar characteristics, it is necessary to know the radar wavelength (or frequency) and the polarization (which describes the orientation of the transmitting and receiving antennae. The radar polarization is specified as the transmit polarization followed by the receive polarization. The nomenclature is defined as “H”, “V”, “R” and “L” for horizontal, vertical, right circular, and left circular respectively. Therefore “HH” indicates horizontal transmit with horizontal receive while “RV” would imply right circular transmit and vertical receive. To address the issue of similar viewing geometries, it is necessary to ensure the beam mode; look direction and pass direction are identical. The beam mode defines the range of incident angles and data resolution. The look direction is defined as the direction of data acquisition (either right or left) relative to the satellite flight direction and the pass direction is defined as ascending when the satellite is moving from south to north or descending when moving from north to south. Ancillary information such as the processing elevation which provides information about the elevation used for the earth model as well as the distance from the radar to the first and last pixel are required for calculations in the downstream processing. The critical baseline represents the maximum allowable distance between the satellites (baseline distance). The closer the baseline distance is to the critical baseline, the lower the coherence (quality) of the interferogram. To address the issue of target stationarity, it is important to know the acquisition time of each data set and the time difference between the two SAR acquisitions. If the time difference is too great, it may result in a poor interferogram due to moving or changing targets. All of the required information mentioned above is automatically extracted and stored in the metadata.

Version 1.1

April 2017

20

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

INSINFO (Analysis of Interferometric Datasets) The INSINFO module is an optional module and can be used to report the interferometric parameters for SAR data sets. To use this module it is important to identify the “reference” data set. The reference data set is the data set from which all the other measurements are made. Although you can only identify a single reference file, it is possible to have many “dependent” data sets. If a single dependent data set is to be used, it can be specified directly by its file name. If many dependent data sets are to be analyzed, (e.g. as part of an interferometric stack) they can be placed in a common directory or alternatively listed in a text file. To execute the INSINFO module, it is necessary to identify the reference file, (FILREF), the dependent file(s) and the name of the text file which will contain the output containing the list of interferometric parameters. If multiple dependent files are specified (in MFILE), the output text file (TFILE) contains information for each reference/dependent combination. The INSINFO module can read data either directly from the vendor supplied format or from the PCIDSK format after data ingest (see SARINGEST). If raw vendor data is specified as input; the viewing geometry parameters such as slant ranges, incident angles and baseline information (first; middle and last) refer to the entire image dataset since no image subset can be specified. If the PCIDSK format (i.e. after reading in the data with SARINGEST) is used as input, the slant range, incident angle range and baseline information correspond to the image subset as specified by the database input window (DBIW) parameter. Additional information and some examples on how to run the INSINFO module can be found in Appendix A.1

SARINGEST (Data Ingest, Extraction of Metadata and Geocoding) The SARINGEST module is optionally used to convert the vendor supplied data sets into the PCIDSK format. The vendor supplied input data set is specified by the FILI parameter. The name to be given to the PCIDSK formatted output file is given by the FILO parameter. Although it is not mandatory to run this module, it does provide some additional flexibility. In particular, the interferometric processing can be limited to a small area of interest. The small area of interest is specified using the database input window (DBIW). The DBIW parameter requires the user to specify the image offsets and the size of the image output in raster coordinates. The XOFFSET value represents the first image column of the subset. Similarly YOFFSET is the first image row of the image subset. The size (in raster units) of the output is given by XSIZE which is the number of columns and YSIZE which represents the number of rows to extract. SARINGEST also provides the user with the ability to calibrate the radar data. Calibration is controlled by the CALIBTYP parameter and is generally used for change detection purposes where it is essential to monitor changes in backscattering. Setting the CALIBTYP to SIGMA converts the backscattering to ground range units. Setting the calibration type to BETA is used to calibrate in slant range units and GAMMA is used to calibrate to the plane perpendicular to slant range. Although calibration is not required for interferometric processing, it is a good idea to calibrate

Version 1.1

April 2017

21

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

all of your interferometric data in a consistent manner. (e.g. all sigma, beta or gamma) because the coherence estimates used in interferometric processing are based upon the amount of backscattering. Using different calibration techniques will result in a statistical bias when certain interferometric parameters such as coherence are calculated. The internal organization of the PCIDSK data can be specified by the DBLAYOUT parameter. The DBLAYOUT parameter is only important for very large data sets because. This parameter determines how much disk space has to be searched to retrieve the data. The recommended setting is TILED256 which organizes the data into local tiles of 256 pixels by 256 lines. The final parameter is the POPTION parameter. The recommended POPTION parameter is AVERAGE (which creates overview pyramids). Selecting AVERAGE will reduce the image loading time and provide smoother looking overviews when using the Geomatica™ display module FOCUS. However; the AVERAGE option does require additional storage space on disk. The SARINGEST module has to be executed for both the reference file and the dependent file. SARINGEST automatically imports and organizes the required metadata for all of the available polarizations. The file’s metadata, geocoding information, and organizational structure can be seen using the FILE PROPERTIES METADATA tab in FOCUS. Additional information and some examples on how to run the SARINGEST module can be found in Appendix A.2

INSCOREG (Automated Image to Image Coregistration) The INSCOREG module is used to automatically co-register the dependent data set (as specified by FILI) to the reference data set (specified by FILREF). Prior to generating an interferogram, it is necessary to align the dependent file with the reference file on a pixel by pixel basis. INSCOREG uses a single reference channel (specified by DBIC_REF) and finds the best image match on the dependent channel (as specified by DBIC_DEP) to generate a RESAMPLED output file (specified by FILO) which contains the identical number of pixels (columns) and lines (rows) as the input reference file and is aligned one-to-one with the reference file. The INSCOREG module selects a number of points from the reference file (as specified by NUMPTS) and searches the corresponding area in the dependent file. The size of the area to search (in raster units) is specified by the search radius (SEARCHR). Within the search area the point (to sub-pixel accuracy) with the highest the probability of being a match is found. If the estimated probability exceeds the MINSCORE value, the point is selected as a candidate point. Once all of the candidate points have been identified, a low order warping function is calculated. Candidate points with absolute residual errors exceeding 0.1 pixel are rejected. The low order polynomial is used to resample the data base input channels (specified by the DBIC parameter) of the dependent file to match the spacing of the reference file on a one-to-one basis. The match statistics are written to the text file specified by the REPORT parameter. Additional information and some examples on how to run the INSCOREG module can be found in Appendix A.3

Version 1.1

April 2017

22

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

INSRAW (Generation of Raw or Adaptively Filtered Interferogram) The final step in the creation of an interferogram is to convolve the signal of the reference dataset (specified by FILREF) with those of the resampled (i.e. co-registered) data set (FILI). The reference input channels specified by DBIC_REF will be convolved with the channels specified by the resampled (i.e. coregistered dependent file) channels (DBIC) and written to the output file (FILO). The user has the option to apply an adaptive filter to the output. The adaptive filter size (FLSZ) must be odd and of size 5 or bigger. If no filtering is to be applied, the FLSZ parameter must be left blank. However; if no filtering is applied, the intensity of each pixel in the interferogram will be one. See “How to display an interferogram” . Any pixels in the reference file which do not overlap with the dependent file are set to NO DATA VALUE. Additional information and some examples on how to run the INSRAW module can be found in Appendix A.4

Example of How to Create an Interferogram We have described all the steps needed to generate an interferogram. The first requirement is to establish whether a pair of SAR data sets constitutes a “good” interferometric pair. To establish this we want to know if they have the same radar characteristics and viewing geometry. The reference and dependent datasets must have identical wavelengths, look direction, pass direction, and beam mode. In addition, the percentage of image overlap should be high and the baseline distances should be low compared to the critical baseline. In our example, we have placed two vendor supplied data sets organized by acquisition date in separate directories in “C:\MyProject”. Within each directory, the vendor files are named “product.xml’. We will designate the dataset from 2014-07-17 as the reference file and the one from 2014-09-03 as the dependent file. To determine if they meet the requirements to form an interferometric pair, we will run INSINFO and write the text output to the file 2014_07_17_2014_09_03_INSAR.rpt. The EASI script is shown below: EASI>FILREF = “C:\MyProject\2014_07_17\product.xml EASI>MFILE =“C:\MyProject\2014_09_03\product.xml EASI>TFILE= “C:\MyProject\2014_07_17_2014_09_03_InSAR.rpt EASI>Run INSINFO

! Identify the reference file ! Identify a single dependent file ! Name of output text file

The output written to the text file (specified by TFILE) is given below. The critical parameters are highlighted. C:\MyProject\2014_07_17\product.xml Acquisition Time: 2014-07-17T12:34:32.276546Z

Version 1.1

April 2017

23

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Wavelength (m): 0.0554658 Beam Mode: U5 Polarizations: HH Look Direction: Right Pass Direction: Descending Near Incidence Angle (degrees): 33.248592 Far Incidence Angle (degrees): 34.551216 First Pixel Slant Range (m): 950606.296000 Last Pixel Slant Range (m): 938357.103000 Pulse Repetition Frequency (Hz): 1620.110596 Processing Elevation (m): 186.092712 Critical Baseline (m): 11745.982920 Nominal ECEF Satellite Position (m): 715559.123599 -2189273.657145 6781850.446254 Reference: C:\MyProject\2014_07_17\product.xml Dependent: C:\MyProject\2014_09_03\product.xml Acquisition Time: 2014-09-03T12:34:32.027653Z Wavelength (m): 0.0554658 Beam Mode: U5 Polarizations: HH Look Direction: Right Pass Direction: Descending Near Incidence Angle (degrees): 33.248783 Far Incidence Angle (degrees): 34.551407 First Pixel Slant Range (m): 950606.296000 Last Pixel Slant Range (m): 938357.103000 Pulse Repetition Frequency (Hz): 1620.110596 Processing Elevation (m): 187.197952 Critical Baseline (m): 11746.067392 Nominal ECEF Satellite Position (m): 715875.375205 -2189150.064812 6781856.525945 Time Difference: 47:23:59:59.751107 Percent Overlap: 98.065612% Ambiguity Height (m): 62.758471 First Line Baseline Distance (m): 287.381451 First Line Parallel Baseline Distance (m): 144.274378 First Line Perpendicular Baseline Distance (m): -248.541752 Middle Line Baseline Distance (m): 287.504478 Middle Line Parallel Baseline Distance (m): 144.339800 Middle Line Perpendicular Baseline Distance (m): -248.646027 Last Line Baseline Distance (m): 287.626860 Last Line Parallel Baseline Distance (m): 144.404925 Last Line Perpendicular Baseline Distance (m): -248.749730

Since the wavelengths, beam modes, look directions, and pass directions for the reference and dependent files are identical and have significant image overlap with baseline distances considerably less than the critical baseline, this pair meets the radar and viewing requirements needed to generate an interferogram. The final quality of the interferogram will depend upon whether the area of interest remained coherent during the 48 days between acquisitions The next step is to select and calibrate the image subset representing the area of interest from the reference and dependent datasets. The easiest way to identify the area of interest is to load the original vendor data set into FOCUS and identify the raster image positions of the upper left and lower right corners for the area of interest. The DBIW offsets (XOFFSET and YOFFSET) will correspond to the upper left corner. The extent of the region of interest is given

Version 1.1

April 2017

24

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

by the XRANGE and YRANGE parameters which are simply the differences in pixel and line coordinates between the lower right and upper left corner. (i.e. XRANGE = XLowerRight - XUpperLeft and YRANGE = YLowerRight - YUpperLeft ). In our example, assume we wish to use the image subset defined in the vendor formatted reference file with upper left position (4900,3700) and lower right position (6900,5700). Similarly we wish to use the area defined by the positions (3050, 1700) for the upper left and (5050, 3700) for the lower right in the vendor supplied dependent file. The SARINGEST module automatically modifies the metadata to align with the region of interest. The acquisition times, slant ranges, and incident angles are automatically adjusted to correspond to the image subset. In addition, assume we want to calibrate the reference and dependent data sets to ground area (i.e. CALIBTYP = sigma). This can all be accomplished using the module SARINGEST. We’ll write the output in tiled format and use the averaging option to create the overviews. The EASI script to extract the area of interest and create the reference file in PCIDSK format is given by the following: EASI>FILI= “C:\MyProject\2014_07_17\product.xml EASI>DBIW = 4900, 3700, 2000, 2000 EASI>FILO = “C:\MyProject\2014_07_17_reference.pix EASI>POPTION =”AVER EASI> DBLAYOUT =”TILED256 EASI>CALIBTYP = “sigma” EASI>Run SARINGEST

! This is the vendor formatted reference file ! Identify the area of interest (offsets and ranges) ! Name to be given to output reference file ! Use averaging to create overviews ! Use 256 x256 tiling for output file ! Calibrate to units of ground area ! Run the module to create the reference file

Since only the first three parameters need to be changed to create the dependent file we can simply use the following to import the dependent file: EASI>FILI= “C:\MyProject\2014_09_03\product.xml EASI>DBIW = 3050, 1700, 2000, 2000 EASI>FILO = “C:\MyProject\2014_09_03_dependent.pix EASI>Run SARINGEST

! This is the vendor formatted dependent file ! Identify the area of interest (offsets and ranges) ! Name to be given to output dependent file ! Run the module to create the dependent file

Now that we have extracted the metadata and defined our area of interest for the reference and dependent files, we need to co-register the dependent file to the reference file. To accomplish this, we use the INSCOREG module. In our example, we will use the first channel of the reference and dependent files to compute the image warping coefficients. We will use 500 evenly distributed points with a search radius of 30 pixels and only accept points with a minimum correlation of 0.85. The image matching statistics will be written to the specified report file. An example of an EASI script to coregister the dependent file to the reference file and create the resampled output is given by: EASI>FILREF= “C:\MyProject\2014_09_03_reference.pix EASI>DBIC_REF = 1

Version 1.1

! This is the calibrated subset of the reference file ! We want to use the first reference channel for image matching

April 2017

25

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

EASI>FILI= “C:\MyProject\2014_09_03_dependent.pix ! This is the calibrated subset of the dependent file EASI>DBIC_DEP = 1 ! We want to use the first dependent channel for image matching EASI>DBIC =1, 2, 3, 4 ! We want to resample all four channels of the dependent file EASI>NUMPTS = 500 ! Use 500 candidate points from the reference file EASI>MINSCORE = 0.85 ! Reject candidate points with scores below 0.85 EASI>SEARCHR = 30 ! Use a search radius of 30 pixels to find a match EASI>FILO = “C:\MyProject\2014_07_17_2014_09_03_resample.pix ! Name to be given to the resampled dependent file EASI>REPORT = “C:\MyProject\2014_07_17_2014_09_03_resample.rpt ! Name to be given to the coregistration report file EASI>Run INSCOREG ! Run the module to create the resampled file

The first section of the report file will contain a list of the candidate points from the reference file with scores exceeding the minimum threshold. The first section lists the reference file image coordinates, followed by their geocoded position and the correlation score. The second section of the report contains the list of candidate points (from the first section) that are within a tenth of a pixel in both the X and Y directions. The third section of the report gives the residual errors of the candidate points exceeding the tenth of a pixel limit. A minimum of 10 points is required to compute the warping coefficients. An example of the report file is given below. The first section of the report file shows all candidates with scores exceeding the MINSCORE Automatic GCP Point Collection FILI GCP Point ID

FILREF PIXEL LINE

X

Y

Z

Score

57784515_4314_001 211.0 10.0 574829.768652 8108823.254799 0.000000 0.844461 57784515_4314_002 211.0 183.0 574948.190127 8108493.119410 0.000000 0.878764 57784515_4314_003 211.0 356.0 575066.160867 8108163.501780 0.000000 0.889292 57784515_4314_006 211.0 875.0 575420.721751 8107174.438704 0.000000 0.886918 57784515_4314_007 211.0 1048.0 575538.796001 8106844.683343 0.000000 0.880030 57784515_4314_008 211.0 1221.0 575657.171892 8106514.865637 0.000000 0.894924 57784515_4314_009 211.0 1394.0 575775.407453 8106184.989056 0.000000 0.888751 57784515_4314_010 211.0 1567.0 575893.209064 8105855.525933 0.000000 0.872189 ………………………………………………………………………………………………………………………………………………… 57784515_4314_420 57784515_4314_421 57784515_4314_422 57784515_4314_423

3823.0 3823.0 3823.0 3823.0

3124.0 585252.715567 8105841.007861 3297.0 585370.840401 8105511.127014 3470.0 585489.036466 8105181.446235 3643.0 585607.552764 8104851.966024

0.000000 0.000000 0.000000 0.000000

0.914336 0.899319 0.894103 0.895230

The second section of the report file lists candidates with scores exceeding the MINSCORE and with absolute residual error less than 0.1 in both X and Y direction. The RMS errors (in pixels) are also provided. 146 GCPs are saved into GCP segment 7, 277 points have been deleted Residual Report for 146 GCPs after refinement: GCP ID X Y ---------------- ----------- ----------------------57784515_4314_023 0.01 -0.00 57784515_4314_024 0.03 -0.04 57784515_4314_025 -0.07 0.03 57784515_4314_028 -0.07 -0.02 57784515_4314_033 -0.03 -0.03 57784515_4314_045 0.02 0.06 57784515_4314_046 0.07 -0.01 57784515_4314_050 0.01 0.03

Version 1.1

April 2017

26

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017 57784515_4314_051 57784515_4314_052 57784515_4314_053 57784515_4314_054 57784515_4314_055 57784515_4314_056 57784515_4314_057 57784515_4314_068 57784515_4314_069 57784515_4314_071 57784515_4314_074

0.00 0.01 -0.02 0.03 0.01 -0.02 0.02 -0.01 0.09 0.06 -0.04

0.04 -0.02 0.03 -0.03 -0.01 0.03 -0.06 0.09 -0.07 -0.07 -0.06

RMS: 0.05 0.04 The third section of the report file lists candidates with scores exceeding the MINSCORE but with absolute residual error more than 0.1 in either the X and/or Y direction Error Report for 277 deleted Points after refinement GCP ID X Y ---------------- ---------------- ----------57784515_4314_001 0.45 -0.25 57784515_4314_002 0.48 0.02 57784515_4314_003 0.42 -0.04 57784515_4314_004 0.29 0.04 57784515_4314_005 0.59 -0.04 57784515_4314_006 0.45 -0.04 57784515_4314_007 0.41 -0.05 57784515_4314_008 0.49 0.02 57784515_4314_009 0.50 0.09 57784515_4314_010 0.40 -0.12 57784515_4314_011 0.36 -0.06

Note: It is important to note that the resampled (dependent) file as specified by the FILO parameter in INSCOREG now represents the dependent data. The final step in the creation of an interferogram requires two co-registered data sets; namely the original reference file and the resampled dependent file. The reference and resampled datasets are used as input to create the raw interferogram. The raw interferogram is the interference pattern created by multiplying the amplitude and phase of the two coherent signals. The raw interferogram still contains the interferometric fringes due to the viewing geometry and the topographic fringes due to topography. The module INSRAW is used to generate the raw interferogram. This module requires the user to specify the name of the reference file (FILREF) and the reference channels to be processed (DBIC_REF), the name of the resampled dependent file (FILI) and the corresponding dependent channels (DBIC). The user has the option to apply an adaptive Lee filter to the interferogram. The size of the filter (FLSZ) must be odd and at least five. The Lee adaptive filter preserves sharp edges and generally provides a “smoother looking” interferogram. If no filtering is to be applied, the FLSZ parameter must be left blank. However; for unfiltered data, all of the pixels in the output interferogram(s) will all have an intensity of one. If the FOCUS view option is set to “grayscale” or “pseudocolor” an unfiltered interferogram will be displayed as a constant

Version 1.1

April 2017

27

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

(black) value since the default for complex data is intensity which will be one everywhere. To properly display an interferogram in FOCUS see the section “How to Display an Interferogram”. An example of an EASI script to generate a raw interferogram using an adaptive filter size of seven is given by: EASI>FILREF= “C:\MyProject\2014_09_03_reference.pix EASI>DBIC_REF = 1,2,3,4 EASI>FILI= “C:\MyProject\2014_07_17_2014_09_03_resample.pix EASI>DBIC =1, 2, 3, 4 EASI>FLSZ =7 EASI>FILO = “C:\MyProject\2014_07_17_2014_09_03_raw.pix EASI>Run INSRAW

!this is the calibrated subset of the reference file !we want to create raw interferograms for all four channels !name of resampled dependent file !use all four channels of the (resampled) dependent file !apply an adaptive filter of size 7x7 !name to be given to the raw interferogram !run the module to create the raw interferogram

Figure 14 shows the interferometric phase in pseudocolor for an unfiltered (left) and 7x7 adaptively filtered dataset (right)

Figure 14a: Unfiltered Raw Interferogram

Version 1.1

April 2017

Figure 14b: Adaptive Lee (7 x 7) Filtering

28

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Example of How to Display an Interferogram using FOCUS The final step is to display the interferogram(s). For an interferometric file containing complex values (designated by FOCUS as C32R for complex 32 bit reals, or as C16U for complex 16 bit unsigned integers) the default is to display the layer as intensity. The data interpretation of the complex layer can be changed by right clicking on the file in the MAPS tab and selecting Data Interpretation and selecting one of the following options Intensity (default), Amplitude (which represents coherence), Phase in Radians, Phase in Degrees, Decibel, Real Part, Imaginary Part. The display will now interpret the complex layer using the selected option. For non-filtered data, the image will initially appear black. This is because the intensity for unfiltered data will be one everywhere. In this case, one of the following Data Interpretations (Phase in Radians, Phase in Degrees, Real Part, Imaginary Part) should be used. Figure 15 shows the identical interferogram interpreted as intensity on the left and as phase in degrees on the right.

Figure 15: Various Data Interpretations (Intensity and Phase in Degrees) for Complex Layers

FOCUS can also display an interferogram as a combination of amplitude and phase. The interferometric amplitude is displayed in black and white while the phase is displayed as a color ramp. To do this, right click the complex valued raster layer, select the view option and select Amplitude and Phase. Alternatively, you can select the grayscale or pseudocolor option for the layer of interest. Figure 16a shows the interferogram displayed as amplitude and phase. The Figure 16b shows the same interferogram viewed as pseudocolor and with the data interpretation (as described previously) set to phase in degrees.

Version 1.1

April 2017

29

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Figure 16a: Complex Layer Displayed by FOCUS as Amplitude and Phase

Version 1.1

Figure 16b: Complex Layer Displayed by FOCUS as Pseudocolor and Interpreted as Phase in Degrees

April 2017

30

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Topographic and Deformation Products Overview The raw interferogram produced during the previous step contains interferometric fringes caused by a combination of satellite viewing geometry, local topography, local subsidence, atmospheric effects and signal noise. If we want to generate a topographic product such as a DEM or DSM, we need to begin by removing the phase due to viewing geometry. If we wish to generate a deformation product such as a subsidence map, we need to remove the phase due to viewing geometry and local topography. The ability to remove interferometric fringes due to viewing geometry depends upon very precise modelling of the satellite orbits. The information required to model the orbit of each satellite is automatically extracted from the metadata. If we wish to create a subsidence product we will also have to remove the interferometric fringes caused by local topography. The ability to do this accurately is highly dependent upon the precision of the user supplied digital elevation model. The ability to remove the effects of viewing geometry and optionally the topographic effects is supplied by the INSTOPO module. This module assumes the orbital information is correct and if a DEM is provided, the module estimates the slant range elevation for each pixel and removes the topographic phase. The metadata of the output is automatically updated to indicate whether the topographic phase has been removed. The module INSADJUST should be executed after INSTOPO. This module removes any residual fringes due to small errors in the orbital model. The level of orbital adjustment is written to the report file. The module INSUNWRAP is used to convert from wrapped phase (which is discontinuous and ranges between   and  ) to continuous “unwrapped” phase. An example showing the difference between wrapped and unwrapped phase is shown in Figure 17.

Figure 17: Wrapped and Unwrapped Phase

Version 1.1

April 2017

31

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

After converting the unwrapped product from cycles of wavelength to distance, the INSUNWRAP module produces either a digital surface model (DSM) or a subsidence map depending upon the interferometric product metadata tag set by INSTOPO. The phase unwrapping introduces a random offset to the DSM or subsidence product. The DSM can be calibrated by providing a few ground control points. The adjustment is applied using the module DEMADJUST. Once the DSM has been generated the final step is to convert the output from slant range to ground range. This is accomplished by the module GEOCODEDEM. The unwrapping process also introduces a random offset to the subsidence product. To eliminate this offset, the user is required to identify the location of a few stable (i.e. non-moving) points. A low order polynomial is fitted through the adjustment points to produce the final result. The adjustment for deformation products is applied by the module INSCALDEFO.

What You Need to Know The next step after the creation of the interferogram is to remove all of the fringes caused by a combination of the radar viewing geometry and local topography. The removal of the fringes caused by the viewing geometry of the interferometric configuration is generally known as “Flat Earth Phase Removal”. For flat earth phase removal it is necessary to calculate the fringe rate (i.e. change in phase) for a specified earth model (ellipsoid) given the baseline (the relative distance between the satellites) as defined by the orbital state vectors. During the image ingest step (see SARINGEST), the orbital state vectors and all of the other parameters required for interferometric processing were stored in the reference and dependent files’ metadata. The state vectors are used to compute both of the reference and dependent orbits as well as the baseline distance as a function of time. By knowing the baseline distance at the time the reference point is acquired and the distance from each satellite to the ellipsoidal model (adjusted for processing elevation), it is possible to calculate the interferometric phase required to correct the interferogram for the “flat earth” effects. If this phase is applied to the raw interferogram, the remaining fringes represent the local topography. As a rule of thumb, each interferometric cycle represents a change in elevation given by the Ambiguity Height parameter. (see INSINFO). It is also possible to remove the interferometric fringes caused by the local topography. To do this requires a digital elevation model (DEM). The original earth model is adjusted to the known elevations extracted from the DEM and adjusted for the processing elevation. The adjusted elevations are then projected from ground range into slant range (distance from the reference satellite). The slant range elevation is used to compute an additional phase correction which can be used to remove interferometric fringes due to topographic effects. If the topographic correction is applied, the remaining fringes represent subtle changes in the earth’s

Version 1.1

April 2017

32

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

surface. These small changes can be used to map local deformation (subsidence) caused by mineral extraction, urban infrastructure, volcanic activity, or glacial rebound. If you are interested in generating digital surface models (DSMs), the removal of the flat earth phase is sufficient and no DEM is required. If you are interested in deformation mapping, a DEM is required as input. In general, the quality of the deformation product is highly dependent upon the quality of the DEM. The final step required to create a topographic or deformation product is to unwrap the interferometric phase. The raw interferometric phase values range (in radians) from   to  . The process needs to determine the “best” integer number of cycles (N ) to add to the phase associated with the two-dimensional unwrapping. This is given mathematically by  unwrap   wrap  2 N . For noisy images, this is a difficult problem and many techniques have been developed. The unwrapping technique we have implemented is built upon the open source package known as SNAPHU™ which is a statistical cost and network flow algorithm. SNAPHU requires a significant amount of computer memory and this should be a consideration when processing interferometric data. Our approach is to extract the pertinent metadata and perform all of the pre-processing steps required to launch SNAPHU. It is worth noting that the phase unwrapping step can take a very long time for areas with low coherence such as water bodies or areas in radar shadow. We have provided the option to mask out areas of low coherence prior to the phase unwrapping step. Masking may significantly reduce the processing time required for unwrapping because the masked areas are simply interpolated. The processing modules to convert raw interferograms into topographic or deformation products are described in the following sections.

INSTOPO (Removal of Flat Earth and Topographic Phase) The INSTOPO module is used to generate either a topographic or deformation interferogram. Topographic interferograms have phase adjustments to remove only the flat earth effects and can subsequently be used to generate digital surface models (DSMs). Deformation interferograms have both the flat earth and topographic phase effects removed and can be used to identify subtle changes in surface elevations (e.g. land subsidence). If a DEM is supplied, the INSTOPO module assumes the output is a deformation product and automatically writes all of the appropriate metadata required for the downstream processing (e.g. the phase unwrapping step). If no DEM is supplied, the assumption is that a topographic interferogram is to be generated and the output metadata is modified accordingly. The INSTOPO module requires the name of the raw interferogram (FILI) which was generated by INSRAW, the list of database input channels (DBIC) to process, (optionally) the file name of the DEM (FILEDEM) and the name of the output file (FILO). Any portion of the interferogram which is outside the boundaries of the DEM will be marked as NO DATA VALUE.

Version 1.1

April 2017

33

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

An example of an EASI script to generate a topographic interferogram from the raw interferogram is given by the following: EASI>FILI= “C:\MyProject\2014_07_17_2014_09_03_raw.pix EASI>DBIC_REF = 1,2,3,4 EASI>FILEDEM = EASI>FILO = “C:\MyProject\2014_07_17_2014_09_03_topo.pix EASI>Run INSTOPO

!this is the file name of the raw interferogram (from INSRAW) !create topographic interferograms for all four channels !leave blank since no topographic correction is to be applied !name to be given to the topographic interferograms !run the module to create the topographic interferogram

An example of the topographic output displayed in FOCUS as amplitude and phase is shown in Figure 18. Each cycle of the interferometric fringe (e.g. from red to red) corresponds to the ambiguity height. An example of an EASI script to generate deformation interferograms from the same raw interferogram is given by the following script: EASI>FILI= “C:\MyProject\2014_07_17_2014_09_03_raw.pix EASI>DBIC_REF = 1,2,3,4 EASI>FILEDEM =“C:\MyProject\Project_DEM.tif EASI>FILO = “C:\MyProject\2014_07_17_2014_09_03_defo.pix EASI>Run INSTOPO

!this is the file name of the raw interferogram !create deformation interferograms for all four channels !use DEM to remove topography leaving only deformation !name to be given to the deformation interferograms !run the module to create the deformation interferogram

Figure 18: Topographic Product with interferometric fringes due to elevation

Version 1.1

April 2017

34

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

An example of the deformation output displayed as amplitude and phase is shown in Figure 19. For deformation products, each cycle of the interferometric fringe (e.g. from red to red) corresponds to a change in elevation of half a wavelength. Note: Since the lower left corner was not in the area covered by the DEM it is automatically flagged with NO DATA VALUE and appears as part of the white background.

Figure 19: Deformation Product with interferometric fringes due to small changes in elevation

Additional information and some examples on how to run the INSTOPO module can be found in Appendix B.1 The topographic and deformation corrections in INSTOPO rely exclusively on the accuracy of the orbital state vectors. Minor errors in the satellite’s velocity vector can result in residual fringes remaining in the interferogram. These residual fringes are due to the positional and drift errors in the calculated orbits of the reference and/or dependent satellite.

INSADJUST (Removal of Residual Orbital Errors) The module INSADJUST is used to remove any residual fringes in the topographic or deformation interferograms. This module uses a Fourier transform to convert the interferogram to spectral components. The module iteratively estimates the horizontal and vertical frequency components and applies a phase correction based upon the adjusted viewing geometry.

Version 1.1

April 2017

35

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

To run the INSADJUST module, the user has to identify the input file (FILI) and the channel to be used for the orbital correction (DBIC_OCS). Even though only one channel is specified to calculate the orbital correction, the orbital phase correction is applied to all channels specified by the DBIC parameter since orbital errors affect the phase of all interferometric channels in an identical manner. The MAXITER parameter specifies the maximum number of Fourier iterations to apply. The default number of iterations is 10 up to a maximum of 20. The program may stop before the specified number of iterations is reached. The program terminates when the applied correction is less than 0.001 cycles across the image; if the calculated orbital adjustment is greater than the previously corrected adjustment (oscillating) or if the specified number of iterations has been reached. If an oscillation is detected, the final (oscillating) correction is not applied. The phase of the specified output file (FILO) is modified to correct for the baseline changes due to the errors in satellite position and drift. The corrections are summarized in the text file specified by the REPORT parameter. The orbital adjustment can be applied to topographic and/or deformation interferograms. An example of the report file is given below. Iteration 1 of 20 Detected residual fringes in cycles: horizontal 1.1334, vertical 1.2635 Equivalent orbital adjustments: Midpoint correction for orbital offset: -3.766145 meters Correction for orbital drift 0.028988 meters per second Iteration 2 of 20 Detected residual fringes in cycles: horizontal 0.1536, vertical 0.2130 Equivalent orbital adjustments: Midpoint correction for orbital offset: -0.510372 meters Correction for orbital drift 0.004887 meters per second Iteration 3 of 20 Detected residual fringes in cycles: horizontal 0.1438, vertical 0.2503 Equivalent orbital adjustments: Midpoint correction for orbital offset: -0.477993 meters Correction for orbital drift 0.005743 meters per second Oscillating correction detected in iteration 3. (Note: Correction (#3) is not applied because vertical frequency is oscillating) Mode adjustment of interferogram phases: -0.50 degrees, -0.00873 radians. Cumulative adjustments in cycles: horizontal 1.2869, vertical 1.4765. Equivalent orbital adjustments: Midpoint correction for orbital offset: -4.276517 meters Correction for orbital drift: 0.033875 meters per second

An example of an EASI script to adjust the orbit for a topographic interferogram is given by: EASI>FILI= “C:\MyProject\2014_07_17_2014_09_03_topo.pix !file name of the topographic interferogram EASI>DBIC_OCS = 1 !use the first channel to determine the adjustment EASI>DBIC=1, 2, 3, 4 !apply the correction to all four channels EASI>FILO = “C:\MyProject\2014_07_17_2014_09_03_topo_adjusted.pix !name of adjusted topographic interferogram EASI>MAXITER = 20 !try up to 20 iterations EASI>REPORT = “C:\MyProject\2014_07_17_2014_09_03_topo_adjusted.rpt !name of report file EASI>Run INSADJUST !create adjusted topographic interferogram

Version 1.1

April 2017

36

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

A comparison of the topographic output derived from the original orbit (left) and adjusted orbit (right) is illustrated in Figure 20.

Figure 20a: Topographic Interferogram with Original Orbit

Figure 20b: Topographic Interferogram with Modified Orbit

An example of an EASI script to adjust the orbit for a deformation interferogram is given by: EASI>FILI= “C:\MyProject\2014_07_17_2014_09_03_defo.pix EASI>DBIC_OCS = 1 EASI>DBIC=1, 2, 3, 4 EASI>FILO = “C:\MyProject\2014_07_17_2014_09_03_defo_adjusted.pix EASI>MAXITER = 20 EASI>REPORT = C:\MyProject\2014_07_17_2014_09_03_topo_adjusted.rpt EASI>Run INSADJUST

!file name of the deformation interferogram !use the first channel to determine the adjustment !apply the correction to all four channels !name of adjusted deformation interferogram !try up to 20 iterations !name of report file !create adjusted topographic interferogram

A comparison of the deformation output derived from the original orbit (left) and adjusted orbit (right) is illustrated in Figure 21.

Version 1.1

April 2017

37

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Figure 21a: Deformation Interferogram Using Original Orbit

Figure 21b: Deformation Interferogram Using a Modified Orbit

Additional information and some examples on how to run the INADJUST module can be found in Appendix B.2

INSUNWRAP (Unwrapping the Interferometric Phase) Once the residual orbital phase errors have been removed from the topographic or deformation interferogram, we are ready to unwrap the interferometric phase. The interferometric phase is cyclic and limited to the range   to  . The INSUNWRAP module is used to convert the wrapped phase to unwrapped phase (which can be directly scaled to relative topographic or relative deformation). INSUNWRAP is a wrapper module that prepares the data, determines the processing parameters, invokes SNAPHU to perform the actual unwrapping, and ingests the results back into the standard PCIDSK format. SNAPHU is a freely available, two-dimensional, statistical cost network flow algorithm. The source code, manuals and references can be found at http://web.stanford.edu/group/radar/softwareandlinks/sw/snaphu The INSUNWRAP module extracts the required metadata from the input data file (FILI) and performs all of the preprocessing steps required by SNAPHU. The metadata determines the output product (FILO). If the topographic phase has not been removed, the output is a digital surface model (DSM). If the topographic phase has been removed the output is a deformation

Version 1.1

April 2017

38

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

map. (see INSTOPO). The type of input (TOPO or DEFO) is stored in the interferometric layer metadata as InSARContent. SNAPHU requires a significant amount of computer memory. The default is to use half of the available computer memory. The user can alter the amount of memory allocated for the phase unwrapping by specifying the number of megabytes to be assigned to the process with the MEMSIZE parameter. If the memory requirements of SNAPHU are too large to fit into the allocated buffer size, the input is automatically tiled and each output tile is unwrapped separately. This may result in obvious artefacts in the unwrapped phase. Due to the large memory requirements, INSUNWRAP can only handle one channel at a time. The channel number to be unwrapped is given by the DBIC parameter. SNAPHU tries to optimally resolve all of the residuals (i.e. incompatible 2D phases) found in the interferogram. Since low coherence areas such as water can have many residuals, unwrapping over water bodies or other low coherence areas can take a considerable amount of time and memory. The unwrapped results in low coherence areas are generally unreliable and can lead to unwrapping errors. To avoid this problem, the user can specify a MASKFILE. The mask file contains geocoded polygons which encompass areas of low coherence. To create a mask file for an area of interest use FOCUS. The task is accomplished by loading the image which covers the area of interest, creating a vector layer (FILE tab  New Vector Layer), delineating the polygon area to be masked out (using the vector editor in FOCUS) and saving the output as a separate file. The geocoded vectors will have a segment number which is given by the MASK parameter. If the mask is contained in the file (FILI) to be unwrapped, you only need to identify the segment number using the MASK parameter without the MASKFILE. If the MASK and MASKFILE parameters are left blank, all of the pixels will be unwrapped. SNAPHU offers two choices for the cost initialization. The first option is the Minimum Spanning Tree (MST) and the second is the Minimum Cost Flow (MCF). In INSUNWRAP the selection is set by the INITALGO parameter. For more information on each option the user should refer to the SNAPHU manual. The PROCMODE parameter specifies which processing mode should be used to launch SNAPHU. The RUN option automatically creates the SNAPHU configuration file based upon the available metadata and image size, preprocesses the data, launches SNAPHU, accepts the output from SNAPHU and reformats it into a PCIDSK file, and updates the PCIDSK metadata. The PREP option simply creates the SNAPHU configuration file in the output directory but does not launch SNAPHU. The user can modify the configuration file (e.g. change the size or number of tiles, modify the overlap between tiles etc.) and then manually execute the command written in the REPORT file. The command will launch SNAPHU. If the PREP option is used, the user can rerun the module with the FINISH option which will copy the output from SNAPHU and reformat to PCIDSK. The INSUNWRAP output will contain two layers. The first layer is the unwrapped phase (stored as polar complex values). The output of the second layer (which is real valued), is relative

Version 1.1

April 2017

39

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

elevation (in meters) for topographic input or relative displacement (in meters) for deformation input.

How to Unwrap a Topographic Interferogram to Generate a DSM An example of an EASI script to unwrap the topographic interferometric phase using a mask file to identify low coherence areas, using the minimum spanning tree algorithm, and 24 Gigabytes of memory is given by the following: EASI>FILI= “C:\MyProject\2014_07_17_2014_09_03_topo_adjusted.pix !file name of the adjusted topo interferogram EASI>DBIC = 1 !unwrap the first channel EASI>MASK =3 !do not unwrap geocoded polygon in this layer EASI>MASKFILE = C:\MyProject\WaterMask.pix” !name of mask file EASI>FILO = “C:\MyProject\2014_07_17_2014_09_03_topo_unwrapped.pix !name to be given to unwrapped topo file EASI>INITALGO= “MST !use the minimum spanning tree algorithm EASI>PROCMODE= “RUN !preprocess, execute and reformat to PCIDSK EASI>MEMSIZE=24000 !use up to 24 Gbytes (24000 Mbytes) of memory EASI>REPORT = “C:\MyProject\2014_07_17_2014_09_03_topo_unwrapped.rpt !name of report file EASI>Run INUNWRAP !create unwrapped topographic interferogram

The configuration file created by INSUNWRAP can be manually edited if the PROCMODE parameter is set to “PREP”. The user should refer to the SNAPHU manual for an explanation of the parameters. An example of the configuration file is given below. Generally it is not advised to change anything except for the Tile Control Parameters since the number of tiles and amount of overlap between tiles is based upon our conservative estimate. These can be modified to reduce (or increase) the number of tiles used in the unwrapping. INFILE LINELENGTH OUTFILE CORRFILE LOGFILE

C:\MyProject/2014_07_17_2014_09_03_topo_adjusted.img 4000 C:\MyProject/2014_07_17_2014_09_03_topo_unwrapped_OutputPhase.img C:\MyProject/2014_07_17_2014_09_03_topo_unwrapped_Correlation.img C:\MyProject/2014_07_17_2014_09_03_topo_unwrapped.log

################### # Runtime options # ################### STATCOSTMODE TOPO INITMETHOD MST VERBOSE FALSE ################ # File formats # ################ INFILEFORMAT FLOAT_DATA CORRFILEFORMAT FLOAT_DATA OUTFILEFORMAT FLOAT_DATA ############################### # SAR and geometry parameters #

Version 1.1

April 2017

40

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

############################### ORBITRADIUS 7162288.611757 EARTHRADIUS 6371000.000000 BPERP 249.735613 TRANSMITMODE REPEATPASS NEARRANGE 942753.396000 DR 1.331 DA 2.040 RANGERES 2.69 AZRES 3.47 LAMBDA 0.0554658 NCORRLOOKS 14.2535 ################################## # Decorrelation model parameters # ################################## RHOSCONST1 0.0001 RHOSCONST2 0.3 RHOMINFACTOR 1 ################ # Tile control # ################ NTILEROW 1 NTILECOL 1 ROWOVRLP 30 COLOVRLP 30 NPROC 1 RMTMPTILE TRUE

Figure 22 shows the unwrapped topographic interferogram displayed as amplitude (coherence) and phase (unwrapped). The unwrapped phase is shown with maximum elevation in red and minimum elevation in purple. The unwrapped phase is directly proportional to the relative elevation.

Figure 22: Unwrapped Topographic Interferogram Version 1.1

April 2017

41

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

The second channel of the output is the computed relative elevations. Figure 23 shows the DEM output using FOCUS as grayscale (left), pseudocolor (middle) and DEM Edit (right)

Figure 23: Unwrapped Relative Elevations shown in FOCUS Grayscale (left), Pseudocolor (middle) and DEM Edit (right)

How to Unwrap a Deformation Interferogram (Subsidence Map) An example of an EASI script to unwrap the deformation interferometric phase using a vector water mask file, the minimum cost flow algorithm, and half of the available memory is given by the following:

EASI>FILI= “C:\MyProject\2014_07_17_2014_09_03_defo_adjusted.pix !file name of the adjusted defo interferogram EASI>DBIC = 1 !unwrap the first channel EASI>MASK = 3 !use segment 3 (waterbodies) from the mask file EASI>MASKFILE =“C:\MyProject\Water_Mask.pix !name of the file containing waterbodies EASI>FILO = “C:\MyProject\2014_07_17_2014_09_03_defo_unwrapped.pix !name to be given to unwrapped defo file EASI>INITALGO= “MCF !use the minimum cost function EASI>PROCMODE= “RUN !preprocess, execute and reformat to PCIDSK EASI>MEMSIZE= !use up to half of available memory EASI>REPORT = “C:\MyProject\2014_07_17_2014_09_03_defo_unwrapped.rpt !name of report file EASI>Run INUNWRAP !create unwrapped deformation interferogram

Figure 24 shows the unwrapped deformation interferogram displayed as amplitude (coherence) and phase (unwrapped). The unwrapped phase is shown with maximum displacement in red and the minimum displacement in purple. The unwrapped phase is directly proportional to the relative subsidence.

Version 1.1

April 2017

42

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Figure 24: Unwrapped Deformation Interferogram

The second channel of the deformation output is the computed relative subsidence estimates. Figure 25 shows the subsidence as grayscale (left), pseudocolor (middle) and DEM Edit (right)

Figure 25: Unwrapped Relative Deformations shown in FOCUS Grayscale (left), Pseudocolor (middle) and DEM Edit (right) Additional information and some examples on how to run the INUNWRAP module can be found in Appendix B.3

Version 1.1

April 2017

43

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

DEMADJUST (Adjustment of Topographic Products) The output from the INSUNWRAP module for topographic input is a DSM product with a random offset and possibly a “tilt”. The random offset is an artefact of the phase unwrapping process since “zero” phase is random and can represent any elevation. The “tilt” may come from an error in the unwrapping process or from atmospheric effects. Each missed cycle in the unwrapped topographic product will add or subtract an elevation corresponding to the AMBIGUITY HEIGHT (see INSINFO). The correction of the elevation values to absolute values requires some “adjustment” points. If the true elevations at specific points are known, it is possible to project the points to the slant range image projection, calculate the difference between the unwrapped elevation and the true elevation and apply an interpolated correction for each pixel. The first step is to create a simple text file which contains the geocoded position and elevation of the known adjustment points. The easiest way to do this is to use the NUMWRIT module. The NUMWRIT module writes the positional information in the format specified by NUMFORM from the input file (FILE) and channel (DBIC) at the user specified spacing (SAMPLING). The output is written to the text file specified by the TFILE parameter. The MASK parameter is used to specify a subarea of the DEM to process and can be left blank. Since the DBGEO parameter has to be specified when generating geocoded output, it should be set to one. The second step is to project the geocoded points generated by NUMWRIT into the projection of the unwrapped elevation model. To accomplish this we use the module VREAD. The text output from NUMWRIT is used as the database vector file (FILV). The points are extracted to a vector segment in the user specified projection (VECUNIT) and stored in the unwrapped topographic file (FILE). Once the vector segment has been added to the unwrapped topographic product (by VREAD), it is important to note the vector segment number. The LASC (Last Segment Created) parameter can be used in a script or you can load the image into FOCUS and note last vector segment number. This value will be required as one of the inputs to the DEMADJUST module as the DBVS (Data Base Vector Segment) value. The DEMADJUST module requires the input file (FILV) containing the vector segment (i.e. the FILE parameter from VREAD). The module also requires the “new” vector segment number (DBVS) and the fieldname (FLDNME). Since we are using the elevations (i.e. Z coordinates) for the correction, the FLDNME parameter should be set to “ZCOORD. The file containing the DEM (FILEDEM) is the original unwrapped topographic product which in our case is identical to the FILV parameter. The input channel (DBEC) is given by the original elevation channel. The adjustment method to be applied is specified by ADJMETHO and the corrected output channel is written to the channel specified by the DBOC parameter. The available options for the ADJMETHO parameter are described in the Help files. Since the offset can be quite large and there may be a “systematic tilt” to the output, it is recommended that all of the input points (the default) be used.

Version 1.1

April 2017

44

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

If the DBEC and DBOC parameters are identical, the original elevations will be overwritten. If you do not want to overwrite the original unwrapped data, it is necessary to add additional channel(s) to the unwrapped topographic product using PCIMOD or FOCUS. To get an estimate of the applied correction which is a combination of unwrapping errors, atmospheric effects, and noise simply subtract the “calibrated elevations” from the original unwrapped elevation. Note: If the datum of the adjustment points (e.g. Mean Sea Level) is different from the datum of the original data set, (e.g. WGS84), the elevation of the input points are automatically modified, but the final output is written using the original datum. Additional information and some examples on how to run the DEMADJUST module can be found in Appendix B.4

Steps to Calibrate Unwrapped Topographic Products In our example, we will use a coarse DEM to generate the control points. The DEM contains elevation values on the first channel. Since we may want to use the output many times (for other data sets) we will generate a text file written in geocoded coordinates. We also wish to suppress the areas of NO DATA. We are not applying a MASK and setting the DBGEO segment to one. An example of how to run NUMWRIT to generate geocoded output, suppressing the NO DATA areas and writing the output to a text file, is given below: EASI>LASC = 0 EASI>FILE= “C:\MyProject\Project_DEM.tif EASI>DBGEO = 1 EASI>DBIC = 1 EASI>MASK = EASI>NUMFORM = “GEOCOORD, Suppress EASI>SAMPLING= 50, 50 EASI>TFILE = “C:\MyProject\Elevation_adjustment.txt EASI> run NUMWRIT

!initialize last vector segment to 0 !use coarse DEM to generate adjustment points !default !use first channel of input DEM !use all of the DEM, do not mask out any area !write points in geocoded format and suppress NO DATA VALUE th th !write text for every 50 pixel and 50 line !name of output text file

The output file (TFILE) will contain three columns. The first two are the calculated geocoded position. The third column contains the image value at the sampled point. Since the input file is a DEM, the image value is the elevation. An example of the output text generated by NUMWRIT is illustrated by: 5.809712590994443E+05 8.106017083988409E+06 3.6708078E+02 5.772212590994443E+05 8.103517083988409E+06 4.2518369E+02 5.784712590994443E+05 8.103517083988409E+06 4.5190347E+02 5.797212590994443E+05 8.103517083988409E+06 4.9227872E+02 ……………………………………………………………………………………………………………….. 5.809712590994443E+05 8.103517083988409E+06 4.1440872E+02

The next step is to copy the adjustment points into the vector segment of the unwrapped topographic product. The input is the text file generated by NUMWRIT. The format (in this

Version 1.1

April 2017

45

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

example) of the input is Universal Transverse Mercator UTM. We want to copy the vector point segment to our unwrapped topographic file. We can (optionally) assign a label to the vector segment and a description. An example of how to run VREAD and copy the adjustment points to the unwrapped topographic output is given below. EASI>FILV= “C:\MyProject\Elevation_Adjustment.txt EASI>VECUNIT = “UTM 16 X D000 EASI>FILE = “C:\MyProject\2014_07_17_2014_09_03_topo_unwrapped.pix EASI>DBSN=”ELEV EASI>DBSD=”Elevation Adjustment Points EASI> run VREAD

!text file with adjustment points !formatted as UTM, can add zone, row and datum !name unwrapped topo file !my name of segment (can be defaulted) !my description of segment !run program to read and copy vector segment

If you do not know which vector segment number (DBVS) has been assigned to the output, you will need to note the value (use FOCUS). In a script you can use the LASC value (Last Segment Created) which was modified by VREAD. Since I do not want to overwrite my original data, I will need to add a channel to the unwrapped topographic product. To do this, run PCIMOD before running DEMADJUST. EASI>FILE = “C:\MyProject\2014_07_17_2014_09_03_topo_unwrapped.pix !name unwrapped topo file EASI>PCIOP=”ADD !operation is to add a channel EASI>PCIVAL=0,0,0,1 !add one floating point channel EASI>run PCIMOD !run the program to add one channel. EASI>FILV=FILE EASI>DBVS=LASC EASI>FLDNME=”ZCOORD EASI>DBEC =2 EASI>DBOC=3 EASI>ADJMETHO = EASI> run DEMADJUST

!the vector segment is in the unwrapped topo file !this is the “recently added” segment number !use the elevation (zcoord) !the input is the uncalibrated unwrapped topo !use channel just added by PCIMOD for output !use all the adjustment points !run the topographicadjustment program

Figure 26 shows the original uncalibrated topographic output (left), the calibrated topographic output (middle) and the applied correction (right). The applied correction compensates for residual errors in elevation, atmospheric effects and noise.

Figure 26: Original Topographic Product (left), Calibrated Topographic Product (middle) and the Applied Topographic Correction (right)

Version 1.1

April 2017

46

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Figure 27 shows the original (unwrapped) topographic output (left), the calibrated topographic output (middle) and the applied correction (right) in DEMEDIT mode.

Figure 27: Original Topographic Product (left), Calibrated Topographic Product (middle) and Topographic Correction (right) shown in DEM EDIT mode

GEOCODEDEM (Convert DSM to Geocoded Product) The calibrated DSM output from DEMADJUST is written in slant range. The final step is to convert the DSM from slant range to a geocoded projection. This is accomplished using the GEOCODEDEM module. This module requires the name of the calibrated slant range topographic product which is specified by the MFILE parameter. The MAPUNITS parameter is a character string that identifies the projection of the output file while the channel containing the slant range DSM is specified by the DBEC parameter. The sample size of the output is specified by the PXSZOUT parameter. The upper left and lower right positions of the output file are given by UPLEFT and LORIGHT respectively. The name to be given to the output file is specified by the FILO parameter. An example of an EASI script to convert from slant range to ground range is: EASI>MFILE=“C:\MyProject\2014_07_17_2014_09_03_topo_unwrapped.pix EASI>DBEC =3 EASI>MAPUNITS=”UTM 16 X D000” EASI>UPLEFT=576025, 8106260 EASI>LORIGHT=581460, 8101360 EASI>PXSZOUT = 2.5 EASI>FILO =“C:\MyProject\2014_07_17_2014_09_03_topo_final.pix EASI> Run GEOCODEDEM

Version 1.1

April 2017

!unwrapped topo file with calibrated elevation !the input is the calibrated channel !want to write the output as UTM !want to define upper left UTM position !want to define lower right UTM position ! want a sample size of 2.5 meters !specify name of output file !convert from slant range to orthorectified

47

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Additional information and some examples on how to run the GEOCODEDEM module can be found in Appendix B.5 Figure 28 shows an example of the original DSM in slant range on the left and the geocoded DSM output on the right. The blank outline denotes the identical screen display. Note the shift of higher elevations away from the radar look direction (which is from the right).

Figure 28a: DSM in slant range projection

Figure 28b: DSM in geocoded projection

INSCALDEFO (Adjustment of Deformation Products) The output from the INSUNWRAP module for deformation products is a product which shows subtle changes in elevation. As with the topographic product, the random offset is an artefact of the phase unwrapping process since “zero” phase is random and can represent any deformation. In addition, there may be an error in the output of the deformation since a “tilt” may arise from an error in the unwrapping process. Each missed cycle will add or subtract an subsidence corresponding to one half of the radar wavelength (see INSINFO). The correction of the deformation values requires some “adjustment” points. If the locations of specific points which are “not moving” are known, it is possible to project the location of the stationary points into the image projection, and apply an interpolated correction for the deformation at each pixel. The module requires as input the unwrapped deformation file (FILE) and the input and output channels (DBIC and DBOC respectively). If the output channel (DBOC) is the same as the input channel (DBIC), the input channel is overwritten with calibrated displacement values. The pixel

Version 1.1

April 2017

48

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

and line locations of the stationary points can be provided in a text file (TFILE).If a single point is entered, it becomes the reference point and all deformation estimates are relative to the point. If multiple points are specified, the correction values over the image are interpolated. Since the deformation product may be “noisy”, the option to average the deformation estimate over a window is provided. The size of the averaging window (CALWIN) must be odd and is independent in the pixel and line direction. The default window size is 1,1 which uses the stationary point only. If CALWIN is set, the “correction value for the stationary point” is based upon the average value within the window size. Alternatively geocoded positions of stationary points can be added as a vector segment using NUMWRIT, and VREAD in a manner similar to that described in the DEMADJUST section. If a geocoded vector segment is available, the FILV, DBVS and FLDNAME parameters can be optionally used. The user must specify one (or both) of TFILE and FILV. The method to be used to compute the adjustment is given by the ADJMETHOD parameter. The available options for the ADJMETHOD parameter are described in the Help files. Since the deformation estimate offset can be large, it is recommended that all of the input points be used. Figure 29: Shows the histogram of the deformation over an area known to be stationary. The histogram indicates a deformation bias of about 2.35 centimeters (0.0235 m). To remove this bias, TFILE can be used to identify the pixel and line coordinates of a reference point which is assumed to be stationary (i.e. deformation = 0). The pixel and line location is given by 2627.5, 1692.5 and is the only point in my adjustment file.

Figure 29: Histogram and Statistics over a “stable” area. Version 1.1

April 2017

49

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

If you do not want to overwrite the original data, you will need to add a channel to the unwrapped topographic product. To do this, run PCIMOD before running INSCALDEFO. EASI>FILE = “C:\MyProject\2014_07_17_2014_09_03_defo_unwrapped.pix !name unwrapped deformation file EASI>PCIOP=”add !operation is to add a channel EASI>PCIVAL=0, 0, 0, 1 !add one floating point channel EASI>run PCIMOD !run the program to add one channel.

The following script will average 201 by 201 area centered on the single pixel and line location given in the text file and use the window average as the offset to be removed for the entire file. EASI>FILE = “C:\MyProject\2014_07_17_2014_09_03_defo_unwrapped.pix EASI>DBIC =2 EASI>DBOC=3 EASI>TFILE = “C:\MyProject\Stationary_points.txt EASI>FILV= EASI>DBVS= EASI>FLDNME= EASI>CALWIN=201,201 EASI>ADJMETHO = EASI> run INSCALDEFO

!name unwrapped deformation file !input is the uncalibrated unwrapped deformation !added by PCIMOD for calibrated output !text file with pixel, line location of stationary pts !no geocoded file is available !no geocoded vector segment is available !no elevation information is used !average around points to reduce noise !use all the adjustment points !run the topographic adjustment program

The histogram of the calibrated file (with offset removed) is given by Figure 30.

Figure 30: Histogram and Statistics over the identical “stable and calibrated” area after random offset is removed

Version 1.1

April 2017

50

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Figure 31 shows the original unwrapped deformation (left) and the calibrated deformation adjusted for offset (right)

Figure 31: Original Deformation Product (left) and Deformation Product Adjusted for Stable Target Constant Offset

Additional information and some examples on how to run the INSCALDEFO module can be found in Appendix B.6

Version 1.1

April 2017

51

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Temporal Analysis of Deformation Products Once a number of calibrated deformation products have been generated, it is possible to stack them in chronological order for temporal analysis. Before stacking the deformation products, it is necessary to coregister all of the data sets of interest. The deformation maps may be derived from one or more reference images, but prior to stacking (see INSSTACK), they must be perfectly coregistered with identical image characteristics (number of pixel and lines), projection and sample sizes.

ORTHO (Orthorectification of Deformation Products) The ORTHO module can be used to coregister all of the deformation maps into a common projection. A complete description of how to run ORTHO is beyond the scope of this manual, but a brief description of the important parameters follows. ORTHO requires the following information. The first is the name of the input file(s) (MFILE). MFILE can be a single file, or a text file listing the names of multiple files which can be processed simultaneously. If the MFILE represents multiple files the output is automatically named. The MFILE parameter has the same characteristics described in the INSINFO section. The DBIC parameter indicates which channels are to be orthorectified. The deformation product generated by INSUNWRAP has two channels. The first is the complex valued unwrapped interferogram and the second is the unwrapped deformation (in meters). The INSCALDEFO module (optionally) added a third real valued channel to correct for the random offset. For our example we will assume this configuration represents the input. The SRCBGD parameter specifies which input pixels are considered as “NO DATA”. Since the interferogram will likely have a portion of the image which contains NO DATA values, we will use the input file’s value. The FILO parameter is the name to be given to the orthorectified output. We can specify the value to be used as the NO DATA VALUE (OUTBGD) in the output. The value -32768 (which is -215) is traditionally used in PCIDSK. The geocoded coordinates of the upper left and lower right of the output are given by ULX, ULY, LRX and LRY respectively. The MAPUNITS parameter defines the projection of the output. The BXPXSZ and BYPXSZ parameters specify the size (in meters) of the output pixels. The FILEDEM is the name of the file containing the DEM and DBEC indicates the DEM channel to be used. The sampling interval (SAMPLING) represents the spacing at which the projection model is computed. The RESAMPLE parameter describes the resampling method. The user should refer to the ORTHO Help File for a more detailed description of the parameters. In this example we will orthorectify a single file. This gives us control over the name of the output. However; it is worth noting that many files can be orthorectified simultaneously. We are going to orthorectify all 3 channels (1 complex and two real) of the input using the input file’s NO DATA VALUE. We explicitly define the map projection and UTM positions of the image corners

Version 1.1

April 2017

52

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

to ensure the output data sets are perfectly aligned. We will write the output at 5 meter spacing in X and Y and compute a position for each pixel using the 16 point SINC interpolator. The ORTHO script is given by the following: EASI>MFILE = “C:\MyProject\2014_07_17_2014_09_03_defo_unwrapped.pix !name unwrapped defo file EASI>DBIC=1, 2, 3 !orthorectify all 3 channels of input file EASI>SRCBGD=”FILE” !recognize input file’s NO DATA VALUES EASI>FILO =“C:\MyProject\2014_07_17_2014_09_03_ortho_defo.pix !name to be given to output file EASI>OUTBGD= -32768 !use standard NO DATA VALUE EASI>ULX = “573650” !set upper left longitude or easting EASI>ULY = “8110150” !set upper left latitude or northing EASI>LRX = “585515” !set lower right longitude or easting EASI>ULY = “8099730” !set lower right latitude or northing EASI>MAPUNITS =”UTM 16 X D000” !set the map projection EASI> BXPXSZ=”5.0” !set output pixel spacing in X to 5.0 meters EASI> BYPXSZ=”5.0” !set output pixel spacing in Y to 5.0 meters EASI>FILEDEM =“C:\MyProject\Project_DEM.tif” !use DEM to correct for topography EASI>DBEC = 1 !use first channel of DEM file EASI>SAMPLING=1 !compute a correction for every pixel EASI>RESAMPLE=”SINC16” !use 16 point sinc function to resample EASI>run ORTHO !run Orthorectification program

For subsequent runs, only the input filename (i.e. MFILE) and output filename (FILO) should be modified for all deformation products to be used for temporal analysis. Figure 32 shows the original calibrated deformation map in slant range (left) and orthorectified deformation (right)

Figure 32: Slant Range Deformation Product (Left) and Orthorectified Deformation Product (right)

Version 1.1

April 2017

53

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

INSSTACK (Temporal Stacking of Deformation Products) The INSSTACK module converts a number of files containing either interferograms and/or terrain displacement (as specified by MFILE) into a single file of chronologically ordered layers. The input channels to be used are controlled by the DBIC parameter. If the DBIC channel points to a complex valued interferogram, the phase of the interferogram is used as input. If the DBIC channel indicates a real valued channel, it is assumed the real valued channel represents the corrected (i.e. edited) line of sight (LOS) displacement in meters which are subsequently converted to unwrapped phase. The output quantity (OQUANT) is specified by the user and indicates which interferometric parameter is generated. This can be one of interferometric phase (PHS), line of sight velocities (VEL), or the line of sight cumulative displacement (DISP). The selected quantity is written to the output file (FILO). The user has to option to apply temporal unwrapping (UNWRAP) to the output. It is important to note that the temporal unwrapping option always assumes the absolute phase difference between each temporally ordered layer does not exceed half a cycle (wavelength). If the input data has been edited or calibrated, the unwrapping option should not be used. The INSTACK module performs the following steps. The temporal midpoint of each layer is calculated. The temporal midpoint is defined as the average time of the reference and dependent file acquisitions. Based upon the midpoint, the input layers are automatically placed in chronological order. If the input layer represents displacement, the displacement is converted to LOS phase otherwise; the interferometric phase is used directly. If the acquisition date of the reference file occurs after the dependent file, the sign of the phase is automatically modified which ensures the output is always written in the time forward direction. If the unwrapping option is selected, the input phase may be adjusted by a number of cycles to ensure the difference in phase from the previous interferometric layer does not exceed half a cycle. (i.e. no phase difference between layers exceeds 180°). If the user selects the phase option PHS as the output, the process is complete when all of the temporally ordered phase values for each pixel have been adjusted for time direction and optionally unwrapped. If the user selects the velocity option (VEL) for the output, the INSSTACK module calculates the line of sight velocity of each pixel. This is accomplished by taking the (possibly adjusted) phase values and converting to a LOS displacement. The velocity is computed by dividing the displacement by the difference in time between the temporal midpoints of the layers. The user has control over both the distance units (DUNITS) and the time units (TUNITS) for the velocity calculations. If the displacement option (DISP) is selected, the cumulative displacement for each layer is written as the integral of the average velocities over the period between temporal midpoints.

Version 1.1

April 2017

54

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

An example illustrating how the INSSTACK output is calculated Assume a sensor with a wavelength of ten centimeters generates four data sets. Suppose the first interferogram covers a sixty day interval with the reference data set acquired on day zero and the dependent data set acquired on day sixty. Suppose the second interferogram covers an eighty day interval and is temporally inverted with a reference acquisition date corresponding to day one hundred and the dependent acquisition occurring on day twenty. For a given coregistered pixel, assume the phase from the first interferogram is 45° and the phase from the second (and temporally inverted) interferogram is -270°. If the UNWRAP parameter is set to “NO” and the user specifies the output as phase, the output for the first and second channels will be 45°, 270° respectively. Note the sign of the phase of the second channel was automatically changed to account for the inverted temporal order of the reference and dependent files. Even though the phase jump is greater than 180°, we assume the phase is correct because the user indicated no unwrapping is to be applied. If the user selects the velocity option with time units of years and distance units of centimeters, the velocity at each midpoint will be calculated in centimeters per year. The output for the velocities is calculated by determining the total displacement for each interferometric layer and dividing by the time interval covered by the interferogram. Note: the temporal intervals between acquisitions do not have to be equally spaced and can overlap in time. The first step is to compute the temporal midpoints of each interferogram. For the two interferograms, these are computed as day thirty (i.e. (0 + 60) / 2) and day sixty (i.e. (100 + 20) / 2) respectively. The total displacement that occurs during the interferometric interval of the first data set is given by the change in phase (i.e. 45°), multiplied by the radar wavelength (10 cm.) divided by 4 (720°). Similarly the total displacement for the interval given by the time between acquisitions for the second interferogram is given by the forward time change in phase (i.e. 270°) multiplied by the radar wavelength (10 cm.) divided by 4 (720°). The computed total displacements in centimeters for the two interferograms are 0.625 cm. (i.e. (45°x 10 cm.) / 720°) and 3.75 cm. (i.e. (270°x 10 cm.) / 720°) respectively. The displacement for the first interferogram covers an interval of 60 days (i.e. days 0-60), while the second covers an interval of 80 days (i.e. days 20-100). To compute to velocity of the first interferogram we convert the displacement over 60 days to a velocity of 3.80 centimeters per year (the units specified by the user). Similarly, the computed velocity for the second output layer is given as 17.1 centimeters per year which is the total displacement of 3.75 centimeters during the 80 day interval converted to centimeters per year. If the cumulative displacement (DISP) option is selected, the displacement computed as the integral of the average velocities multiplied by the time difference between the interferometric midpoints of each layer. Note it is assumed the cumulative displacement is initially zero. In our example the first cumulative displacement representing day 30 would be 0.156 cm. This is derived by computing the average velocity for the interval between temporal midpoints (day 0 to 30) and converting to centimeters. (i.e. (0+3.80)/2 * 30/365.25). The second layer representing

Version 1.1

April 2017

55

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

the integrated displacement up to day 60; would be the initial displacement (i.e. 0.156 for the first thirty days) plus the displacement between days 30 and 60. This is computed as 0.859 cm. which is calculated as (3.80 + 17.1)/2 centimeters per year multiplied by the difference in midpoints (i.e. (60 - 30) / 365.25 days/year). The cumulative displacements are written as 0.156 cm. for day 30 and 1.015 cm. (i.e. 0.156 + 0.859) at day 60. Using the same example, with the UNWRAP parameter set to “YES” with the OQUANT parameter set to phase, the output for the first and second channels will be 45°, -90° respectively. This occurs because original phase value is -270°. Because the input interferogram is temporally inverted, we change the sign (i.e. 270°) to be in forward time. However; the phase difference between the two layers exceeds the 180° limit and is lowered by one cycle (from 270° to -90°) to be within the 180° limit. This phase change will affect all of the subsequent calculations. The displacement for the second interferogram is now calculated using the modified phase as (-90° *10 cm) / 720° = -1.25 cm. Since this displacement occurs over an interval of 80 days, the velocity for the second interferogram is now computed as -1.25 cm /80 days * 365.25 days/year = -5.70 cm./year. The displacement over the interval of interest for the second layer is now computed as (3.80 - 5.70) / 2 * (60 - 30) /365.25 which gives -0.078 cm. The cumulative displacement for the second output layer is the sum of the displacement for the first channel (0.156 cm.) and the second channel (-0.078 cm.) giving a total displacement of 0.078 cm. The previous example is summarized by the table below. Interferogram Layer

Reference Date (days)

Dependent Date (days)

Temporal Midpoint (days)

Original Phase (degrees)

Select Unwrap Option

Temporally Adjusted Phase 45 270

Velocity cm./year at midpoint 3.80 17.1

Cumulative Displacement cm. at midpoint 0.156 1.105

1 2

0 100

60 20

30 60

45 -270

No No

1 2

0 100

60 20

30 60

45 -270

Yes Yes

45 -90

3.80 -5.70

0.156 0.078

An example of the content of a text file (STACK.TXT) which can be used as input to the INSSTACK module is shown below. In our example, we have calibrated the deformation products generated by INSUNWRAP and calibrated using INSCALDEFO. The input files listed in the text file have the original unwrapped interferogram in the first channel, the uncalibrated deformation in the second channel, and the calibrated deformation in the third channel. To ensure there are no large “jumps” between layers of deformation, we can set the temporal unwrapping to “YES”. Temporal unwrapping ensures that the deformation between layers does not exceed half the radar wavelength. If the time difference between deformation layers is large or you are very confident in the calibrated results, you can set the unwrap parameter to “NO” which will take each deformation estimate “as is”. We also want to use the calibrated deformation to generate our stacked deformation products. We can override the DBIC parameter by placing the channel to be used (i.e. 3) after the input file names separated by a semicolon.

Version 1.1

April 2017

56

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

STACK.TXT 2014_06_23_2014_07_17_defo_unwrapped_ortho.pix;3 2014_06_23_2014_08_10_defo_unwrapped_ortho.pix;3 2014_07_17_2014_08_10_defo_unwrapped_ortho.pix;3 2014_07_17_2014_09_03_defo_unwrapped_ortho.pix;3 2014_08_10_2014_09_03_defo_unwrapped_ortho.pix;3

Assume we want to estimate the deformation velocities in centimeters per year; the input to the INSSTACK is given by: EASI>MFILE =”C:\MyProject\Stack.txt EASI>DBIC = EASI>FILO =“C:\MyProject\Velocity_stack.pix EASI>OQUANT =“vel” EASI>UNWRAP= “yes” EASI>DUNITS= “cm” EASI>TUNITS= “year” EASI>REPORT= C:\MyProject\Velocity_stack.rpt”

!text file containing list of coregistered deformation maps !don’t need to set since channel of interest is given in MFILE !name to be given to output file (will contain 5 layers) !want to estimate velocity values at temporal midpoints !temporally unwrap the data !set distance units to cms !set time units to years !name of report file

An example of the report file is shown below. The report consists of the following information. The report file gives some preliminary information such as number of files used and the time used as the time offset (i.e. the time representing time zero.) The report also lists the information for each FILE in chronological order based upon the temporal midpoint of the interferogram DATE. The TimeTU is the midpoint time since the start time in the requested units (e.g. years) and should increase for each layer. The RangeTU is the time range (in years) between the acquisition times of the reference and dependent files. LambdaDU is the radar wavelength in the units requested by the user. Channel is the channel used as input. Format indicates whether the complex valued input was represented as Cartesian or Polar.

Images: 5 Start: 2014-06-23T12:34:31.897981Z

!five images used as input !this represents time zero

File: 2014_06_23_2014_07_17_defo_unwrapped_ortho.pix Date: 2014-07-05T12:34:31.766221Z TimeTU: 0.032854 RangeTU: 0.065708 LambdaDU: 5.546580 Channel: 3 Format: Polar No Data: -32768. ImageChn: N/A

!first file in chronological order !temporal midpoint of the file !number of years between layer midpoint and time zero !time in years between reference and dependent acquisitions !radar wavelength in centimeters !input channel for this file !input was complex data in polar format !value assigned to NO DATA areas

File: 2014_06_23_2014_08_10_defo_unwrapped_ortho.pix Date: 2014-07-17T12:34:31.781343Z TimeTU: 0.065708 RangeTU: 0.131417 LambdaDU: 5.546580 Channel: 3 Format: Polar No Data: -32768. ImageChn: N/A

!second file in chronological order

Version 1.1

April 2017

57

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

File: 2014_07_17_2014_08_10_defo_unwrapped_ortho.pix Date: 2014-07-29T12:34:31.649583Z TimeTU: 0.098563 RangeTU: 0.065708 LambdaDU: 5.546580 Channel: 3 Format: Polar No Data: -32768. ImageChn: N/A

!third file in chronological order

File: 2014_07_17_2014_09_03_defo_unwrapped_ortho.pix !fourth file in chronological order Date: 2014-08-10T12:34:31.510091Z TimeTU: 0.131417 RangeTU: 0.131417 LambdaDU: 5.546580 Channel: 3 Format: Polar No Data: -32768. ImageChn: N/A File: 2014_08_10_2014_09_03_defo_unwrapped_ortho.pix Date: 2014-08-22T12:34:31.525214Z TimeTU: 0.164271 RangeTU: 0.065708 LambdaDU: 5.546580 Channel: 3 Format: Polar No Data: -32768. ImageChn: N/A

!fifth file in chronological order

An example showing the 5 layers representing the temporal midpoint surface velocities is shown in Figure 33.

July 5

July 17

July 29

Aug 10

Aug 22

Figure 33: Surface Velocities at Specified Dates

If we want to estimate the cumulative surface deformation (in centimeters) the input to the INSTACK is given by: EASI>MFILE =”C:\MyProject\Stack.txt EASI>DBIC =

Version 1.1

!text file containing list of coregistered deformation maps !don’t need to set since channel of interest is given in MFILE

April 2017

58

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

EASI>FILO =“C:\MyProject\Displacement_stack.pix EASI>OQUANT =“displ” EASI>UNWRAP= “no” EASI>DUNITS= “cm” EASI>TUNITS= “year” EASI>REPORT= C:\MyProject\Displacement_stack.rpt”

!name to be given to output file (will contain 5 layers) !estimate cumulative displacement at temporal midpoints !use the calibrated displacements !set distance units to cms !set time units to years !name of report file

An example of the cumulative displacement (since the initial time of June 23) calculated at each temporal midpoint is shown in Figure 34.

July 5

July 17

July 29

Aug 10

Aug 22

Figure 34: Cumulative Displacement at Temporal Midpoints since time zero (June 23) Additional information and some examples on how to run the INSSTACK module can be found in Appendix C.1

Temporal Analysis of Interferometric Stacks The final step is the temporal analysis of the interferometric stack. This can be accomplished by loading the file created by the INSSTACK module into FOCUS. The output file from INSSTACK contains the metadata required for the analysis. This includes the start time of the stack, the temporal midpoints, the distance units (DUNITS) and the temporal units (TUNITS). If ancillary image files are to be used for spatial context, they should be loaded into FOCUS as well. To invoke the Temporal Analysis Tool (TAT) the output from INSSTACK must be loaded into FOCUS. To “select” the stack file, click on the stack file name in the FOCUS MAP tab. Once you have selected the file, go to ANALYSIS and select the InSAR Temporal Chart option. The TAT will appear. The TAT displays the X axis in the stack’s time units. These are extracted from the metadata and were specified by TUNITS in INSSTACK). The output parameter (i.e. one of phase, velocity or cumulative displacement) is displayed in the Y direction in units specified by DUNITS or DUNITS/TUNITS in INSSTACK. The axis labels appear along the left edge and at the bottom of the chart. The name of the selected stack file and name of the Y parameter is shown at the top of the chart. The user can modify the Y range of values by clicking on the Min Y and Max Y boxes and typing in the desired range. The chart will display the parameter values

Version 1.1

April 2017

59

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

corresponding to the cursor position. The geocoded position corresponding to the cursor location is also displayed.

An example of the SAR Chart output (cumulative displacement) is shown in Figure 35

File Name and Parameter Type Graph of Temporal Output Values

Parameter in user specified units

Time in user specified units

User Selected Minimum and Maximum range of Y Values

Geocoded Location of Cursor

Figure 35: Temporal Analysis Tool

Version 1.1

April 2017

60

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Appendix A A.1 InSAR Data Analysis Module Name:

INSARINFO

Processing Requirement:

Optional

Purpose: INSARINFO is a module designed to provide pertinent interferometric parameters such as acquisition dates, baseline distances, incident angles for the reference file and the dependent file(s). If multiple dependent files are to be evaluated, the file names can be listed in a text file. The module supports the original vendor data or files which have been imported using SARINGEST.

User Interface: 

FILREF (required) o Name of the reference file



MFILE: (required) o List of dependent SAR data sets



TFILE: (default = Interferometric_Report.txt”) o Name of Text File to be generated as output

Input Example: FILREF= “C:\MyProject\Reference_data.xml” MFILE = “C:\MyProject\Insar_Files.txt” TFILE = “C:\MyProject\Report.txt” Run InSARInfo where the text file “Insar_Files.txt” contains the reference file (listed first) and a list of dependent files such as C:\MyProject\2014_06_21\product.xml C:\MyProject\2014_07_15\product.xml C:\MyProject\2014_08_08\product.xml …………………………………………………………….. C:\MyProject:\2014_09_01\product.xml

Version 1.1

April 2017

(dependent file #1) (dependent file #2)

(dependent file #N)

61

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Example of Output: C:\MyProject\2013_06_26\product.xml Acquisition Time: 2013-06-26T13:32:42.217153Z Wavelength (m): 0.0554658 Beam Mode: U6 Polarizations: HH Look Direction: Right Pass Direction: Descending Near Incidence Angle (degrees): 34.073978 Far Incidence Angle (degrees): 35.356373 First Pixel Slant Range (m): 959043.270900 Last Pixel Slant Range (m): 946535.863900 Pulse Repetition Frequency (Hz): 1606.794556 Processing Elevation (m): 44.190601 Critical Baseline (m): 12217.054929 ECEF Satellite Position (m): 272016.955367 -2025974.934942 6864205.134259 C:\MyProject\2013_06_26\product.xml C:\MyProject\2014_06_21\product.xml Acquisition Time: 2014-06-21T13:32:37.331065Z Wavelength (m): 0.0554658 Beam Mode: U6 Polarizations: HH Look Direction: Right Pass Direction: Descending Near Incidence Angle (degrees): 34.072884 Far Incidence Angle (degrees): 35.355320 First Pixel Slant Range (m): 959043.270900 Last Pixel Slant Range (m): 946535.863900 Pulse Repetition Frequency (Hz): 1606.794556 Processing Elevation (m): 44.246700 Critical Baseline (m): 12216.565824 ECEF Satellite Position (m): 272000.416113 -2025939.974206 6864225.330988 Time Difference: 359:23:59:55.113912 Percent Overlap: 99.702718% Height Ambiguity (m): 419.790527 First Line Baseline Distance (m): 37.102894 First Line Parallel Baseline Distance (m): 8.637146 First Line Perpendicular Baseline Distance (m): 36.083576 Middle Line Baseline Distance (m): 37.928850 Middle Line Parallel Baseline Distance (m): 8.812191 Middle Line Perpendicular Baseline Distance (m): Last Line Baseline Distance (m): 38.768223 Last Line Parallel Baseline Distance (m): 8.987138 Last Line Perpendicular Baseline Distance (m): 37.712153 Note: The dependent file information is repeated for each file in the MFILE list.

Version 1.1

April 2017

62

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

A.2 InSAR Data Ingest Module Name:

SARINGEST

Processing Requirement:

Optional

Purpose: SARINGEST is a generic module designed to import the specified SAR data set into the PCIDSK image format. The PCIDSK format is a standardized image structure which includes additional InSAR metadata describing acquisition parameters, orbit, calibration and ancillary geocoding information. For single look complex (SLC) data, SARINGEST stores each polarization as a single layer of complex values for subsequent processing.

User Interface: (recommended settings *)

Version 1.1



FILI (required) o Input File Name



DBIW (default is all) o Processing Window (specified in raster units)  XOFFSET, YOFFSET, XSIZE,YSIZE



POPTION (default is averaged) o Pyramid options  Off  Nearest Neighbor  Averaged *



DBLAYOUT (default is band) o Internal data organization  Pixel  Band  Tiled256*



CALIBTYP (default is none) o Calibration Type  Sigma*  Beta  Gamma  None

April 2017

63

PCI Geomatics Enterprises Inc. 

InSAR User Guide Geomatica 2017

FILO (required) o Output File Name

Input Example: FILI = “C:\InSAR\product.xml” FILO = “C:\InSAR\output.pix” POPTION = “Average” DBLAYOUT = “Tiled256” CALIBTYP = “Sigma” Run SARINGEST

Example of Output:

Geocoded Image File

File metadata

Version 1.1

April 2017

64

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

A.3 Coregistration of Reference and Dependent Images Module Name:

INSCOREG

Processing Requirement:

Mandatory

Purpose: INSCOREG is a module designed to automatically resample and co-register the dependent file to the reference file. The module automatically acquires control points, removes outliers, and resamples the dependent file to match one-to-one with the reference file. The metadata of the resampled file is updated to reflect the modified geolocation information. Non-overlapping areas are flagged as “NO DATA” values. The (single) channel specified by DBIC_REF and DBIC_DEP are the channels to be used from the reference and dependent files for the image matching and generation of control points. The control points are used to derive the low order polynomial mapping which is applied to all channels in the DBIC list. The default values for DBIC_REF and DBIC_DEP is 1. The resampled output is written in the order specified by the DBIC parameter. Once the match points have been found, INSCOREG rejects any candidate point which is not within 0.1 pixels (in both the X and Y direction) of the expected location. The correction is a second order polynomial which requires a minimum of 10 match points. The resampling kernel is a 2 dimension 16 point complex valued sinc (i.e. sin (x)/x) function

User Interface:

Version 1.1



FILI_REF: (required) o Reference SAR data set



DBIC_REF (required, default = 1) o The channel of reference file to be used for the image to image matching



FILI: (required) o Dependent SAR data set



DBIC_DEP (required, default = 1) o The channel of dependent file to be used for the image to image matching



DBIC (required) o Specifies the order in which the dependent values are to be written. o The reference DBIC_REF channel is compared to the dependent DBIC_DEP channel to compute the one-to-one mapping.

April 2017

65

PCI Geomatics Enterprises Inc. o

InSAR User Guide Geomatica 2017

The mapping is applied to all channels listed in DBIC and written sequentially to the output file (FILO)



NUMPTS (default = 500) o This represents the number of candidate points to generate. o The candidate points are evenly distributed over the reference file.



MINSCORE (default =0.75) o This represents the minimum acceptable correlation score. o The range of permitted values is 0.72 to 1.0



SEARCHR (default =50) o This represents the maximum search radius (in pixels) to look for an image match.



FILO: (required) o File name to be given to resampled dependent data set o This output file will contain the number of channels specified by DBIC o This file cannot already exist



REPORT (optional) o Name of text file containing summary of GCP information o Default = INSAR_GCP.rpt

Input Example (single polarization): FILI_REF = “C:\MyProject\Reference_file.pix” DBIC_REF = 1 DBIC_DEP = 1 DBIC = 1 NUMPTS = 300 MINSCORE = 0.85 SEARCHR = 20 FILI = “C:\MyProject\Dependent_file.pix” FILO =”C:\MyProject\Resampled_Dependent_file.pix” REPORT =”C:\MyProject\Coregistration_Report.txt” Run INSCOREG

Version 1.1

April 2017

66

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Input Example (multiple polarizations): FILI_REF = “C:\MyProject\Multi_Pol_Reference_file.pix” DBIC_REF = 4 DBIC_DEP = 4 NUMPTS = 400 MINSCORE = 0.8 SEARCHR = 100 FILI = “C:\MyProject\Mult_Pol_Dependent_file.pix” DBIC = 1, 2, 3, 4 FILO = “C:\MyProject\Mult_Pol_Resampled_Dependent_file.pix” REPORT = “C:\MyProject\Multi_Pol_Coregistration_Report.txt” Run INSCOREG

Example of Report Output: Residual Report for 110 GCPs after refinement: GCP ID X Y -------------------------- ----------12C1CF5A_4314_016 -0.02 0.05 12C1CF5A_4314_018 0.03 -0.04 12C1CF5A_4314_019 0.01 -0.07 12C1CF5A_4314_020 0.01 0.00 12C1CF5A_4314_032 0.09 0.08 12C1CF5A_4314_033 0.06 0.05 12C1CF5A_4314_038 -0.06 -0.08 12C1CF5A_4314_039 0.05 0.03 12C1CF5A_4314_043 -0.06 0.03 12C1CF5A_4314_044 0.07 -0.01 12C1CF5A_4314_046 -0.08 0.06 12C1CF5A_4314_048 0.07 -0.07 12C1CF5A_4314_053 0.08 -0.06 12C1CF5A_4314_054 -0.09 0.06 12C1CF5A_4314_055 -0.05 -0.09 12C1CF5A_4314_057 -0.10 -0.00 12C1CF5A_4314_058 -0.01 -0.04 12C1CF5A_4314_060 -0.03 -0.07 12C1CF5A_4314_061 -0.04 -0.04 12C1CF5A_4314_314 -0.03 -0.05 12C1CF5A_4314_316 0.06 0.03 12C1CF5A_4314_317 -0.04 0.05 12C1CF5A_4314_319 -0.01 0.01 12C1CF5A_4314_320 0.05 0.02 12C1CF5A_4314_321 0.09 -0.08 RMS:

Version 1.1

0.05

0.05

April 2017

67

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Example of (unaltered) Reference and Coregistered Resampled File

Unaltered Reference File

The resampled dependent file is automatically co-registered pixel by pixel with the reference file. Matching points will have identical raster coordinates and geolocation information. Version 1.1

April 2017

68

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

A.4 Generation of Raw (or Filtered) Interferogram Module Name:

INSRAW

Processing Requirement:

Mandatory

Purpose: INSRAW is the module required to produce the raw interferogram using the specified reference file and resampled dependent file generated by the INSCOREG module. The user has the option to filter the interferogram using an adaptive Lee filter by specifying a filter size (odd) in pixels for the FLSZ parameter. If a filter size is specified it must be five or bigger. Note: If the default (i.e. no filtering) FLSZ is used, the resulting interferogram should be displayed as either “Amplitude and Phase” or interpreted as “Phase (in radians or degrees) “ in FOCUS because the default complex display (intensity) of the interferometric output will be constant (one) and will appear in FOCUS as solid black.

User Interface:

Version 1.1



FILREF: (required) o Reference SAR data set



DBIC_REF (required) o Channel(s) of reference file containing reference SLC data



FILI: (required) o Resampled dependent SAR data set (generated by INSCOREG)



DBIC (required) o Channel(s) of resampled dependent file containing reference SLC data



FILO: (required) o Name of output file containing the raw interferogram(s) o The file cannot already exist



FLSZ (default = none) o If specified, filter size in pixels (N x N where N is odd and five or bigger) of adaptive Lee filter

April 2017

69

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Input Example (no filtering): FILREF = “C:\MyProject\Reference_file.pix” DBIC_REF = 1 FILI = “C:\MyProject\Resampled_Dependent_file.pix” DBIC = 1 FILO =“C:\MyProject\Raw_interferogram.pix” FLSZ = Run INSRAW

Input Example (with adaptive filtering): FILREF = “C:\MyProject\Reference_file.pix” DBIC_REF = 1 FILI = “C:\MyProject\Resampled_Dependent_file.pix” DBIC = 1 FILO =“C\MyProject:\Filtered_interferogram.pix” FLSZ = 7 Run INSRAW

Example of Outputs:

Unfiltered Raw Interferogram Displayed as Amplitude and Phase

Version 1.1

April 2017

Filtered (7x7) Raw Interferogram Displayed as Amplitude and Phase

70

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Appendix B B.1 Flat Earth and Topographic Correction Module Name:

INSTOPO

Processing Requirement:

Mandatory

Purpose: The INSTOPO module extracts the reference and dependent file orbit data and calculates the required phase adjustments to correct for systematic flat earth effects projected to the specified earth model (which is generally the WGS-84 ellipsoid). Variations in phase caused by topography are also corrected if a digital elevation model is supplied. The flat earth correction is applied for the ellipsoid adjusted to the processing elevation (as extracted from the metadata). If a digital elevation model is specified (optional), the DEM elevations are automatically adjusted to compensate for the nominal processing elevation which is extracted from the metadata. The modified elevations are projected to the reference file slant range prior to calculating the topographic phase correction. When no DEM is specified, only the systematic flat earth phase is removed and the output contains fringes due to topographic effects. These fringes can be converted to elevation using the INSUNWRAP module. . Residual fringes (if any) due to orbital error can be corrected using INSADJUST.

User Interface:

Version 1.1



FILI: (required) o Raw interferogram (from INSRAW)



DBIC: (required) o List of channels to be corrected



FILDEM: (optional) o The name of the file containing the digital elevation model. The elevation is automatically adjusted to compensate for the nominal processing elevation given by the metadata parameter “NominalLocation_Height.



FILO: (required) o Name of output file containing the flat earth or topographically corrected interferogram o The file cannot already exist

April 2017

71

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Input Example: Flat Earth Correction Only: FILI = “C:\MyProject\Raw_Interferogram.pix” DBIC = 1 FILDEM = FILO =”C:\MyProject\Flat_Earth_Correction_Only.pix” Run INSTOPO

Example of interferogram with only flat earth phase removed.

Version 1.1

April 2017

72

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Input Example: Flat Earth and Topographic Correction: FILI = “C:\MyProject\Raw_Interferogram.pix” DBIC = 1 FILDEM = “C:\MyProject\Digital_Elevation_Model.tif” FILO = “C:\MyProject\Topographic_Correction.pix” Run INSTOPO

Example of interferogram with flat earth and topographic phase removed. Remaining fringes are due to orbital error and atmospheric effects.

Version 1.1

April 2017

73

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

B.2 Orbital Correction and Residual Fringe Removal Module Name:

INSADJUST

Processing Requirement:

Optional

Purpose: The INSADJUST module automatically removes low frequency residual phase caused by orbital drift and/or baseline errors. The module iteratively estimates both the horizontal and vertical spectral residuals and applies the correction based upon the viewing geometry for the interferometric pair. Note: INSADJUST uses large amounts of memory; systems with limited memory may perform poorly when processing large images. The input file must be in slant range geometry and the incident angle array segment must be available in the input file’s metadata. The Fast Fourier Transform (FFT) is derived for the whole image, and its X and Y sizes are set to the smallest power-of-two equal to or higher than the input interferogram file pixel and line dimensions. The transform is stored in a temporary file, and is deleted at the end of the run. The algorithm iterates to detect the residual low-frequency signal in the FFT domain, and then removes it in the image domain. The maximum number of iterations is specified in the MAXITER parameter. The iterative correction process is monitored, and may terminate earlier if the current-iteration adjustment is smaller than 0.001 cycles per image (checked separately for the horizontal and vertical image dimensions). The process may also terminate for an oscillating correction; this condition is detected when the derived current iteration adjustment is larger than the adjustment from the previous iteration. For oscillating solutions the output image reflects the correction applied in the iteration prior to the one that diverged. INSADJUST assumes that most of the deterministic effects are already removed from the interferogram. They include the systematic flat Earth and terrain effects. The remaining lowfrequency fringes are usually caused by inaccuracies in satellite orbits and by slow changes in the interferometric baseline, as well as by inaccuracies in the already applied corrections. In typical processing scenarios, the INSADJUST module for orbital correction (if it is to be applied) should be used after executing INSTOPO. The orbit adjusted interferogram should have no, or strongly attenuated, low-frequency noise (residual fringes), and average phase of the output should be a value close to zero radians.

Version 1.1

April 2017

74

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

User Interface: 

FILI: (required) o Specifies the name of the PCIDSK file from which interferogram data is to be read



DBIC_OCS (required) o Specifies the input channel to be used to determine the number of residual fringes. o DBIC_OCS must be a complex-valued channel



DBIC (required) o Specifies the input interferometric channels to be corrected. The channels specified by DBIC must be complex-valued.



FILO: (required) o Specifies the name of the PCIDSK file to receive the interferometric data with the adjusted phase values. o The file cannot already exist



MAXITER (default = 10) o Specifies the maximum number of iterations for removing the residual fringes, to a maximum of 20. The iterative correction process is monitored and may terminate earlier, if convergence or oscillation is detected. The default value is 10, but in most cases the process converges in fewer steps.



REPORT (optional) o Name of text file containing summary of orbital correction information

Input Example: FILI = “C:\MyProject\Topo_Corrected_file.pix” DBIC_OCS =1 DBIC = 1 FILO = “C:\MyProject\Orbit_Corrected_file.pix” MAXITER = 20 REPORT =“C:\MyProject\Orbit_Correction.rpt” Run INSADJUST

Version 1.1

April 2017

75

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Example of Output: Iteration 1 of 20 Detected residual fringes in cycles: horizontal 1.3710, vertical -4.9419 Iteration 2 of 20 Detected residual fringes in cycles: horizontal 0.0313, vertical -0.0186 Iteration 3 of 20 Detected residual fringes in cycles: horizontal 0.0015, vertical 0.0010 Iteration 4 of 20 Detected residual fringes in cycles: horizontal 0.0001, vertical -0.0000 Correction converged in 4 iterations. Mode adjustment of interferogram phases: 66.50 degrees, 1.16064 radians. Cumulative adjustments in cycles: horizontal 1.4039, vertical -4.9595 Equivalent orbital adjustments: Midpoint correction for orbital offset: -0.091254 meters Correction for orbital drift -0.010891 meters per second

Example of an interferogram with residual fringes due to errors in orbital position

Version 1.1

April 2017

Example of an interferogram with residual fringes automatically removed by INSADJUST

76

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

B.3 Phase Unwrapping Module Name:

INSUNWRAP

Processing Requirement:

Mandatory

Purpose: The INSUNWRAP module unwraps the interferometric phase by performing a two dimensional integration of the wrapped ( ,  ] phase values to generate the unambiguous phase. The unwrapping is based upon the Statistical cost, Network flow, Algorithm for Phase Unwrapping (SNAPHU) open source code and is described in Ref [1]. The algorithm computes the most likely unwrapped solution based upon the observable data (coherence). The INSUNWRAP module extracts the pertinent information from the PCIDSK metadata to generate a configuration file to be used by SNAPHU. The INSUNWRAP module automatically performs all the steps necessary to launch SNAPHU. The unwrapped output is written as complex polar values (i.e. amplitude and unwrapped phase) rather than Cartesian co-ordinates (i.e. real and imaginary). If the topographic phase has already been removed in INSTOPO, the metadata of the unwrapped output is tagged as a deformation product. If no DEM is specified in INSTOPO, the output is flagged as a topographic product. A vector mask can be used to identify areas of known low coherence. The areas within the geocoded vector polygons are marked as regions which should not be unwrapped. The output phase for these regions is interpolated. Note: The elapsed time required by SNAPHU is affected by  Size of image to be unwrapped  Size and distribution of non-coherent areas  Amount of available computer memory [Ref 1]

C. W. Chen, H. A. Zebker, “Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization,'' Journal of the Optical Society of America A, vol. 18, pp. 338-351 (2001).

User Interface: 

Version 1.1

FILI (required) o Specifies the name of the PCIDSK file from which interferogram data is to be read.

April 2017

77

PCI Geomatics Enterprises Inc.

Version 1.1

InSAR User Guide Geomatica 2017



DBIC (required) o Specifies the input interferogram to unwrap. DBIC must be a complexvalued channel, and have the InSARContent metadata tag set to Interferogram.



MASK (optional) o Specifies the bitmask segment or vector layers to be used for masking. If the MASKFILE parameter is not specified, the mask layer must exist in the input file.



MASKFILE (optional) o Specifies the name of the file containing the geocoded vector polygons. If the MASKFILE parameter is specified, the mask value must be set.



FILO (required) o The name of a PCIDSK (.pix) output file containing the unwrapped interferogram in the amplitude and unwrapped phase format. o The output file will contain only one channel. This is a limitation of SNAPHU.



INITALGO (required) o SNAPHU cost initialization algorithm to use, MST or MCF for the minimum spanning tree or minimum cost flow algorithm, respectively. Please consult the SNAPHU manual for more details.



PROCMODE (required) o INSUNWRAP processing mode. The following modes are recognized: 

RUN: the complete unwrapping process is executed. It includes creating the input and output raster files, setting up the configuration file with parameters for SNAPHU, executing SNAPHU, and ingesting its results into the output file in FILO.



PREP: the data are prepared for the execution of SNAPHU. The input and output raster files and the configuration file with parameters for SNAPHU are created, the working directory and the command to execute are written to the REPORT destination, and the module terminates. This mode allows adjustments to be made to the automatically generated configuration file (except for the file names and formats), but it requires a subsequent manual execution of SNAPHU.



FINISH: this mode assumes that SANPHU was set up in the PREP mode and then manually executed. The phase

April 2017

78

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017 unwrapped by SNAPHU is verified and then ingested into the FILO. This mode can be also used if SNAPHU fails for some reason, but later succeeds after adjustments to its configuration parameters.

The processing mode can be specified by its first letter only. Note: the adjustments to SNAPHU parameters must not modify the input and output file names, their line length, and formats. The INSUNWRAP processing parameters in the FINISH mode must be the same as in the previous execution of INSUNWRAP in the RUN or PREP mode, except for the PROCMODE parameter itself. 

MEMSIZE [input] o Specifies the amount of physical memory, in megabytes, to use for image buffers and for SNAPHU. The default is 50 percent of the total system memory. An error occurs if the amount requested exceeds the total physical memory available in the system. o The tiling scheme in SNAPHU is derived according to available memory, and the number of bytes per pixel (BPP) used by it. When the full image is processed at once, SNAPHU requires about 350 BPP. When tiling is used, SNAPHU requires about 2100 BPP.

Input Examples: Topographic Example using Vector Masks: Note: The topographic fringes in FILO were not removed by INSTOPO FILI = “C:\MyProject\Wrapped_Topographic_Interferogram.pix” DBIC = 1 MASK = 3 MASKFILE = “C:\MyProject\Masked_Out_Areas.pix” FILO = “C:\MyProject\Unwrapped_Topographic_Interferogram.pix INITALGO = “MCF” PROCMODE = “Run” MEMSIZE = 20000 REPORT =“C:\MyProject\Unwrapped_Topographic_Interferogram.rpt Run INSUNWRAP

Version 1.1

April 2017

79

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Example of a wrapped interferogram containing topographic information

Example showing the unwrapped topographic interferogram

Input Example: Deformation Example without Vector Masks: Note: The topographic fringes in FILO have already been removed by INSTOPO FILI = “C:\MyProject\Wrapped_Deformation_Interferogram.pix” DBIC = 1 MASK = MASKFILE = FILO = “C:\MyProject\Unwrapped_Deformation_Interferogram.pix INITALGO = “MCF” PROCMODE = “Run” MEMSIZE = 20000 REPORT =“C:\MyProject\Unwrapped_Deformation_Interferogram.rpt Run INSUNWRAP

Version 1.1

April 2017

80

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Example of a wrapped interferogram containing deformation information

Version 1.1

April 2017

Example showing the unwrapped deformation information

81

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

B.4 Adjustment of Topographic Products Module Name:

DEMADJUST

Processing Requirement:

Optional

Purpose: This module applies a smooth correction to the unwrapped topographic output produced by INSUNWRAP. The unwrapping process introduces a random offset to the topographic product. This module adjusts the raster DEM to minimize the difference between the unwrapped output and the control point elevations stored in the vector layer. The adjustment algorithm offers a number of mechanisms to remove statistical outliers based upon the absolute difference or standard deviation of the difference between control point and input.

User Interface:

Version 1.1



FILV: (required) o Name of input file containing the vector segment



DBVS (required) o Layer containing the vector segment (created by NUMWRIT and VREAD)



FLDNME: (default = ZCOORD) o Specifies the Field Name containing the elevation values



FILEDEM (required) o Name of DEM file to be adjusted o Usually the same as FILV



DBEC: (required) o Channel containing the DEM to be adjusted



DBOC (required) o Channel containing the adjusted DEM o This can be equal to DBEC



ADJMETHO o Use all elevation points (default) o “Under X” (Only use elevation points with differences less than X)

April 2017

82

PCI Geomatics Enterprises Inc. o o o

InSAR User Guide Geomatica 2017

“CLAMP X” (Use all elevation points but clamp differences to a maximum of X “UNDERSIGMA X” (Only use elevation points with differences under X sigma “CLAMPSIGMA X” (Use all elevation points but clamp differences to X sigma

Input Example FILV = “C:\MyProject\Unwrapped_Topographic_Interferogram.pix” DBVS = 5 FLDNME = FILEDEM= C:\MyProject\Unwrapped_Topographic_Interferogram.pix” DBEC= 2 DBOC=3 ADJMETHO =”Undersigma 2.3” Run DEMADJUST

Output Example

Uncalibrated Topographic Product

Version 1.1

April 2017

Calibrated Topographic Product

83

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

B.5 Conversion from Slant Range to Orthorectified Product Module Name:

GEOCODEDEM

Processing Requirement:

Optional

Purpose: This module converts the digital surface model from slant range to geocoded output. The user has control over the sample size and projection of the output. The process is similar to orthorectification with the exception that the elevation values used as input are extracted from the available file.

User Interface: 

MFILE: (required) o Name of input file(s) to be converted from slant range to ground range



DBEC o Layer containing the (calibrated) elevations



DBIW o Processing Window (specified in raster units)  XOFFSET, YOFFSET, XSIZE,YSIZE



FILO: (required) o Name of output file containing the geocoded output.



MAPUNITS o Specifies the projection of the output file. If this is not specified the projection of the input is used. o The standard definitions are “Long/Lat” or “UTM mm r Dnnn” where mm is the zone number, r is the zone row and Dnnn is the earth model. PXSZOUT o The resolution of the output file. If the output is in LONG/LATs the resolution is specified in degrees for UTM output is specified in meters.



Version 1.1

April 2017

84

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Input Example MFILE = “C:\MyProject\Unwrapped_Topographic_Interferogram.pix” DBEC = 3 MAPUNITS = “UTM 16 X D000” FILO= “C:\MyProject\Final_Topographic_Product.pix” Run GEOCODEDEM

Output Example Note: Higher elevations are shifted to the right

DSM in slant range

Version 1.1

Geocoded DSM

April 2017

85

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

B.6 Adjustment of Deformation Products Module Name:

INSCALDEFO

Processing Requirement:

Optional

Purpose: This module adjusts all the displacement values (from INSUNWRAP) to be zero for all of the points specified by the TFILE (in pixel and line co-ordinates) and/or the geocoded ground positions specified in the vector file (FILE). If the window size (CALWIN) is defaulted, the exact value at the specified point is used to estimate the correction; otherwise the average value of the window centered on the point is used to compute the correction factor.

User Interface:

Version 1.1



FILE: (required) o The name of the file containing the displacement values



DBIC (required) o The channel number of the uncalibrated displacement values



DBOC (required) o The channel number of the calibrated displacement values o Can be equal to DBOC



TFILE: (optional) o Text file containing the image coordinates of the adjustment points



FILV (optional) o Name of the input file containing the geocoded vector point layer



DBVS (optional) o Segment number of the geocoded vectors



FLDNME: (default = ZCOORD) o Field name of elevation values



CALWIN(default = 1,1) o Window size of points to be used for adjustment



ADJMETHO o Use all elevation points (default)

April 2017

86

PCI Geomatics Enterprises Inc. o o o o

InSAR User Guide Geomatica 2017

“Under X” (Only use elevation points with differences under X) “CLAMP X” (Use all elevation points but clamp differences to a maximum of X “UNDERSIGMA X (Only use elevation points with differences under X sigma “CLAMPSIGMA X” (Use all elevation points but clamp differences to X sigma

Input Example FILE = “C:\MyProject\Unwrapped_Deformation_Interferogram.pix” DBIC = 2 DBOC= 3 TFILE = “C:\MyProject\Image_Adjustment_Points.txt” FILV = “C:\ MyProject\Geocoded_Adjustment_Points.pix DBVS = 5 FLDNME = CALWIN = 201,201 ADJMETHOD = “ClampSigma 2.5” Run INSCALDEFO

Output Example

Uncalibrated Deformation Product

Version 1.1

April 2017

Calibrated Deformation Product

87

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Appendix C C.1 Create a Stack of Displacements or Velocities Module Name:

INSSTACK

Processing Requirement:

Optional

Purpose: The module converts a number of co-registered SAR interferograms into a time ordered stack of phase values, accumulative displacements or velocities. It also produces a file of statistical values including average coherence and best fit linear velocity. The stacks can be produced from slant range interferograms or orthorectified interferograms. The input interferograms must be coregistered on a one-to-one basis and in the same projection. The user has the option to apply temporal unwrapping. Temporal unwrapping limits the maximum phase jump between temporal layers (i.e. in the time direction) to   . The user has control over the time and distance units of the temporally ordered interferometric stack. The module also produces ancillary files containing the mean and standard deviations of the output products.

User Interface:

Version 1.1



MFILE: (required) o Name of text file containing list of input images



DBIC: (required) o Channel number to be processed from list of input images o Note: the channel number can be overridden in the text file specified by MFILE



FILO: (required) o The file name of the output file containing the temporally ordered phase, velocities or cumulative displacements



OQUANT o Output Quantity one of  Phase, Velocity, Cumulative Displacement (default)



UNWRAP o Unwrap phases

April 2017

88

PCI Geomatics Enterprises Inc.  

InSAR User Guide Geomatica 2017

Yes (default) implies no phase jump exceeding 180°is permitted between layers. No implies take the unwrapped value “as is”. Note: A single-interferogram is not unwrapped, regardless of the setting of this parameter.



DUNITS (required) o Distance Units one of  mm, cm, m, km



TUNITS (required) o Time Units one of  sec, min, hour, day, month, year



REPORT o Name of report file

Input Displacement Example: MFILE = “C:\MyProject\List_of_interferograms.txt” DBIC =1 FILO = “C:\MyProject\Interferometric_Stack_of_ Displacements.pix” DUNITS = “cm” TUNITS =”year” OQUANT= “disp” UNWRAP= “Yes” Run INSSTACK

Input Velocity Example: MFILE = “C:\MyProject\List_Of_interferograms.txt” DBIC = 1 FILO = “C:\MyProject\Interferometric_Stack_of_ Velocities.pix” DUNITS = “cm” TUNITS = “year” OQUANT= “vel” UNWRAP= “Yes” Run INSARSTACK

Version 1.1

April 2017

89

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Examples of Output:

Sequential Series of Cumulative Displacements computed at Interferometric Temporal Midpoints

Sequential Series of Velocities computed at Interferometric Temporal Midpoints

Ancillary Temporal Products

Best Fit and Standard Deviation of Cumulative Displacement

Version 1.1

April 2017

Average Velocity and Standard Deviation of Velocity

90

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

C.2 Temporal Analysis Tool Module Name: Maps TabSelect FileAnalysis InSAR Temporal Chart Processing Requirement:

Optional

Purpose: This is a new feature for FOCUS that supports temporal InSAR analysis. The X axis values represent the temporal information and are generated by reading the temporal information (date and time) stored in the layer’s metadata. The Y values are extracted by burrowing into the raster position selected by the user.

User Interface: 

FILI: (required) o Name of file containing interferometric stack of data

Example of Velocity Output (centimeters per year)

Example of output from Temporal Analysis Tool

Version 1.1

April 2017

91

PCI Geomatics Enterprises Inc.

InSAR User Guide Geomatica 2017

Complex Value Display in FOCUS Module Name:

FOCUS

Processing Requirement:

Optional

Purpose: This is a new feature for FOCUS that supports visualization of complex values. The amplitude of the complex layer is displayed in black and white while the phase information is displayed as a semi-transparent (30%) color ramp.

User Interface: 

FILI: (required) o Name of file containing complex valued data.

Example of Output:

Visualization of Complex Valued Data Showing Magnitude and Phase Simultaneously

Version 1.1

April 2017

92