Заголовок
Источник: http://picxxx.info
Ссылка на PDF: http://picxxx.info/pml.php?action=GETCONTENT&md5=7e99343e11040e6b4322e83e605d21d0
Конец заголовка
Landsat 7
Image Assessment System (IAS)
Geometric Algorithm Theoretical Basis
Document
IAS Geometric ATBD Version 3
TABLE OF CONTENTS
1. INTRODUCTION
1
1.1 Purpose
1
1.2 Applicable Documents
2
2. OVERVIEW AND BACKGROUND INFORMATION
6
2.1 Experimental Objective
6
2.2 Historical Perspective
6
2.3 Instrument Characteristics
7
2.4 Ancillary Input Data
7
3. ALGORITHM DESCRIPTION
8
3.1 Theoretical Description
3.1.1 Landsat 7 ETM+ Viewing Geometry Overview
3.1.2 Coordinate Systems
3.1.3 Coordinate Transformations
3.1.4 Time Systems
3.1.5 Mathematical Description of Algorithms
3.1.5.1 ETM+ Instrument Model
3.1.5.1.1 Detector/Focal Plane Geometry
3.1.5.1.2 Scan Mirror Model
3.1.5.1.3 Scan Line Corrector Mirror Model
3.1.5.2 Landsat 7 Platform Model
3.1.5.2.1 Sensor Pointing (Attitude and Jitter) Model
3.1.5.2.2 Platform Position and Velocity
3.1.5.3 Level 1 Processing Algorithms
3.1.5.3.1 PCD Processing
3.1.5.3.2 MSCD Processing
3.1.5.3.3 Platform/Sensor Model Creation
3.1.5.3.3.1 Smoothing Low Frequency Attitude Information
3.1.5.3.3.2 Combining Low and High frequency Attitude Information
3.1.5.3.4 Line of Sight Generation and Projection
3.1.5.3.5 Correction Grid Generation
3.1.5.3.6 Resampling
3.1.5.3.6.1 Detector Delays
3.1.5.3.7 Precision Correction
3.1.5.3.7.1 Precision Correction Outlier Detection
3.1.5.3.7.2 Inverse Mapping
3.1.5.3.7.3 Precision Ephemeris Update
3.1.5.3.8 Terrain Correction
3.1.5.4 Geometric Calibration Algorithms
3.1.5.4.1 Sensor Alignment Calibration
3.1.5.4.2 Focal Plane Calibration
3.1.5.4.3 Scan Mirror Calibration
3.1.5.5 Geometric Characterization Algorithms
Version 3.2
ii
8
8
10
16
19
22
22
22
23
24
24
24
25
26
28
32
34
48
56
65
73
78
91
94
105
109
110
112
116
116
123
127
135
7/9/98
IAS Geometric ATBD Version 3
3.1.5.5.1 Band to Band Registration
3.1.5.5.2 Image to Image Registration
3.1.5.5.3 Geodetic Accuracy
3.1.5.5.4 Geometric Accuracy
3.1.6 Variance or Uncertainty Estimates
3.1.6.1 Level 1 Processing Accuracy
3.1.6.2 Geometric Calibration Accuracy
3.1.6.2.1 Sensor Alignment Calibration Accuracy
3.1.6.2.2 Scan Mirror Calibration Accuracy
3.1.6.2.3 Band Placement Calibration Accuracy
3.1.6.3 Geometric Characterization Accuracy
3.1.6.3.1 BandtoBand Registration Accuracy
3.1.6.3.2 ImagetoImage Registration Accuracy
4. CONSTRAINTS, LIMITATIONS, ASSUMPTIONS
Version 3.2
iii
135
143
150
151
154
154
158
158
160
163
165
165
167
170
7/9/98
IAS Geometric ATBD Version 3
TABLE OF FIGURES
FIGURE 3.1.11 ETM+ GEOMETRY OVERVIEW
FIGURE 3.1.12 ETM+ DETECTOR ARRAY GROUND PROJECTION
FIGURE 3.1.21 ETM+ FOCAL PLANE COORDINATE SYSTEM
FIGURE 3.1.22 ETM+ SENSOR COORDINATE SYSTEM
FIGURE 3.1.23 ORBITAL COORDINATE SYSTEM
FIGURE 3.1.24 EARTH CENTERED INERTIAL (ECI) COORDINATE SYSTEM
FIGURE 3.1.25 EARTH CENTERED ROTATING (ECR) COORDINATE SYSTEM
FIGURE 3.1.26 GEODETIC COORDINATE SYSTEM
FIGURE 3.1.5 1 MODEL INITIALIZATION
FIGURE 3.1.5 2 RECTIFICATION AND RESAMPLING
FIGURE 3.1.5 3 PRECISION CORRECTION SOLUTION
FIGURE 3.1.5 4 PRECISION/TERRAIN CORRECTION
FIGURE 3.1.5 5 SCAN MIRROR PROFILE DEVIATIONS
FIGURE 3.1.5 6 EFFECT OF ROLL JITTER ON LINE OF SIGHT DIRECTION
FIGURE 3.1.5 7 MAGNITUDE RESPONSE OF GYRO AND ADS
FIGURE 3.1.5 8 MAGNITUDE RESPONSE OF GYRO PLUS ADS
FIGURE 3.1.5 9 MAGNITUDE RESPONSE OF GYRO + MAGNITUDE RESPONSE OF ADS
FIGURE 3.1.5 10 ATTITUDE PROCESSING NETWORK
FIGURE 3.1.5 11 INPUT IMAGE GRIDDING PATTERN
FIGURE 3.1.5 12 EXTENDED PIXELS AND SCAN ALIGNMENT
FIGURE 3.1.5 13 CALCULATION OF SCAN GAP
FIGURE 3.1.5 14 SCAN MISALIGNMENT AND GAP
FIGURE 3.1.5 15 EXTENDED SCAN LINES
FIGURE 3.1.5 16 CUBIC SPLINE WEIGHTS
FIGURE 3.1.5 17 INVERSE MAPPING WITH “ROUGH” POLYNOMIAL
FIGURE 3.1.5 18 “ROUGH” POLYNOMIAL – FIRST ITERATION
FIGURE 3.1.5 19 RESULTS MAPPED BACK TO INPUT SPACE
FIGURE 3.1.5 20 NEAREST NEIGHBOR RESAMPLING
FIGURE 3.1.5 21 DETECTOR DELAY DEFINITION
FIGURE 3.1.5 22 RESAMPLING WEIGHT DETERMINATION
FIGURE 3.1.5 23 DEFINITION OF ORBIT REFERENCE SYSTEM
FIGURE 3.1.5 24 LOOK VECTOR GEOMETRY
FIGURE 3.1.5 25 TERRAIN CORRECTION GEOMETRY
FIGURE 3.1.5 26 DOQ QUARTERQUAD IMAGES COVERING ABOUT 20 TM SCANS
FIGURE 3.1.5 27 LEGENDRE DIFFERENCES
FIGURE 3.1.5 28 DIFFERENCE BETWEEN ACTUAL MIRROR PROFILES
FIGURE 3.1.5 29 CORRELATION ERRORS
FIGURE 3.1.5 30 CORRELATION STANDARD DEVIATIONS
FIGURE 3.1.5 31 ROAD CROSS SECTION EXAMPLE
FIGURE 3.1.5 32 CORRELATION ERRORS
FIGURE 3.1.6.2.1 1 ALIGNMENT ESTIMATION ACCURACY VS. NUMBER OF SCENES
FIGURE 3.1.6.2.3 1 BANDTOBAND TEST POINT NUMBER VS. ACCURACY
Version 3.2
iv
9
10
11
12
13
14
15
16
26
27
27
28
44
47
57
58
58
59
77
86
87
88
88
89
89
90
90
90
92
94
95
99
115
133
134
135
146
147
148
149
162
165
7/9/98
IAS Geometric ATBD Version 3
1. Introduction
1.1 Purpose
This document describes the geometric algorithms used by the Landsat 7 Image
Assessment System (IAS). These algorithms are implemented as part of the IAS Level 1
processing, geometric characterization, and geometric calibration software components.
The overall purpose of the IAS geometric algorithms is to use Earth ellipsoid and terrain
surface information in conjunction with spacecraft ephemeris and attitude data, and
knowledge of the Enhanced Thematic Mapper Plus (ETM+) instrument and Landsat 7
satellite geometry to relate locations in ETM+ image space (band, scan, detector, sample)
to geodetic object space (latitude, longitude, and height). These algorithms are used for
purposes of creating accurate Level 1 output products, characterizing the ETM+ absolute
and relative geometric accuracy, and deriving improved estimates of geometric calibration
parameters such as the sensor to spacecraft alignment. This document presents
background material describing the relevant coordinate systems, time systems, ETM+
sensor geometry, and Landsat 7 spacecraft geometry as well as descriptions of the IAS
processing algorithms.
The Level 1 processing algorithms include:
1. Payload Correction Data (PCD) processing;
2. Mirror Scan Correction Data (MSCD) processing;
3. ETM+/Landsat 7 sensor/platform geometric model creation;
4. sensor line of sight generation and projection;
5. output space/input space correction grid generation;
6. image resampling;
7. geometric model precision correction using ground control;
8. terrain correction.
These algorithms are discussed in detail in section 3.1.5.3.
The geometric calibration algorithms, discussed in section 3.1.5.4, include:
1. ETM+ sensor alignment calibration;
2. focal plane calibration (focal plane bandtoband alignment);
3. scan mirror calibration.
The geometric characterization algorithms, discussed in section 3.1.5.5, include:
1. band to band registration;
2. image to image registration;
3. geodetic accuracy assessment (absolute external accuracy);
4. geometric accuracy assessment (relative internal accuracy).
Version 3.2
1
7/9/98
IAS Geometric ATBD Version 3
1.2 Applicable Documents
Reference 1
Landsat 7 System Data Format Control Book (DFCB) Volume IV  Wideband Data,
Revision H, prepared by Martin Marietta Astro Space, document number 23007702IVH,
dated 26 February 1998.
Reference 2
Landsat 7 System Program Coordinates System Standard, Revision B, prepared by Martin
Marietta Astro Space, document number PS23007610B, dated 2 December 1994.
Reference 3
Landsat 7 System Image Assessment System (IAS) Element Specification (Revision 1),
NASA Goddard Space Flight Center, document number 43015010011, dated February
1998.
Reference 4
Theoretical Basis of the Science Data Processing Toolkit Geolocation Package for the
ECS Project, prepared by the EOSDIS Core System Project, document number 445TP002002, dated May 1995.
Reference 5
Landsat 7 System Specification, Revision H, Goddard Space Flight Center, document
number 430L0002H, dated June 1996.
Reference 6
Press, Teukolsky, Vetterling, and Flannery, Numerical Recipes in C, 2nd edition,
Cambridge University Press, 1992.
Reference 7
DMA TR 8350.2A, DMA Technical Report, Supplement to Department of Defense
World Geodetic System 1984 Technical Report, prepared by the Defense Mapping
Agency WGS84 Development Committee, dated December 1, 1987.
Reference 8
Snyder, John P., Map Projections  A Working Manual, United States Geological Survey
Professional Paper 1395, U. S. Government Printing Office, Washington, 1987.
Reference 9
Landsat 7 Space Segment Technical Description Document  Enhanced Thematic
Mapper Plus, prepared by Hughes Santa Barbara Research Center, May 1994, updated by
NASA Goddard Space Flight Center, September 1995.
Reference 10
Sehn, G. J. and S. F. Miller, “LANDSAT Thematic Mapper Attitude Data Processing”,
presented at the AIAA/AAS Astrodynamics Conference, Seattle, Washington, August 2022, 1984.
Reference 11
Landsat to Ground Station Interface Description, Revision 9, Goddard Space Flight
Center, Greenbelt, Maryland, January 1986.
Version 3.2
2
7/9/98
IAS Geometric ATBD Version 3
Reference 12
Park, S.K., R.A. Schowengerdt, Image Reconstruction by Parametric Cubic Convolution,
Computer Vision, Graphics and Image Processing, v.23, no.3, September 1983.
Reference 13
Rosborough, G.W., D.G. Baldwin and W.J. Emery, "Precise AVHRR Image Navigation",
IEEE Transaction on Geoscience and Remote Sensing, Vol. 32, No. 3., May 1994.
Reference 14
Rao, C.R., "Linear Statistical Inference and Its Applications", John Wiley and Sons, Inc.,
1973.
Reference 15
Sahin, M., P.A. Cross, and P.C. Sellers, "Variance Component Estimation Applied to
Satellite Laser Ranging", Bulletin Geodesique, 66: 284295, 1992.
Reference 16
Yuan, Dahning, "Optimum Determination of the Earth's Gravitational Field from
Satellite Tracking and Surface Measurements", Dissertation, The University of Texas at
Austin, 1991.
Reference 17
Snyder, John P., Space Oblique Mercator Projection Mathematical Development, United
States Geological Survey Bulletin 1518, U. S. Government Printing Office, Washington,
1981.
Reference 18
LAS 6.0 Geometric Manipulation Package Overview Document, USGS EROS Data
Center.
Reference 19
LAS 6.0 Programmer’s Document, USGS EROS Data Center
Reference 20
Escobal, Pedro R., Methods of Orbit Determination, 2nd edition, John Wiley & Sons,
Inc., New York, 1976.
Reference 21
Dvornychenko, V.N., Bounds on (Deterministic) Correlation Functions with Application
to Registration, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.5,
n.2, Mar 1983
Reference 22
Steinwand, D. and Wivell, C., Landsat Thematic Mapper Terrain Corrections in LAS,
HSTX InterOffice Memo OAB821, USGS EROS Data Center, August 12, 1993.
Reference 23
Brown, R.G., Introduction to Random Signal Analysis and Kalman Filtering, John Wiley
and Sons, New York, 1983.
Reference 24
Version 3.2
3
7/9/98
IAS Geometric ATBD Version 3
Landsat 7 System Calibration Parameter File Definition, Revision 1, 43015010020
NASA
Goddard Space Flight Center, Greenbelt, Maryland, January 1998.
Reference 25
Burt P. and Adelson E., The Laplacian Pyramid as a Compact Image Code, IEEE
Transactions on Communications, vol. Com31, no. 4, April 1983.
Reference 26
Prakash, A. and Beyer, E.P., Landsat D Thematic Mapper Image Resampling for Scan
Geometry Correction.
Reference 27
Beyer, E.P., An Overview of the Thematic Mapper Geometric Correction System
Reference 28
Wolberg, G., Digital Image Warping, IEEE Computer Science Press, 1990.
Reference 29
Gonzalez, R.C. and Wintz, P., Digital Image Processing, 2nd Edition, AddisonWesley
Publishing Company, 1987.
Reference 30
Fischel, D., "Validation of the Thematic Mapper Radiometric and Geometric Correction
Algorithms" IEEE Transactions on Geoscience and Remote Sensing, Vol. GE22 No. 3
May 1984.
Reference 31
Golub, G.H. and Van Loan, C.F., Matrix Computations (2nd ed.), John Hopkins University
Press, Baltimore, 1989.
Reference 32
Landsat 7 ETM+ Geometric Calibration Plan, EROS Data Center, December 2, 1996
Reference 33
Goodman, D.J. and Carey, M.J., Nine Digital Filters for decimation and Interpolation,
IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 25, April 1977
Reference 34
Haykin, Simon, “Adaptive Filtering Theory”, Prentince Hall, 1991.
Reference 35
Landsat 7 Celestial Sensor Assembly, Including Star Transit Processing Electronics and
Celestial Sensor Sun Shield, Martin Marietta Astro Space, PS23007911A, 9 September
1994.
Reference 36
Version 3.2
4
7/9/98
IAS Geometric ATBD Version 3
Cook, R. Dennis, “Detection of Influential Observations in Linear Regression”,
Technometrics, Volume 19, Number 1, February, 1977.
Reference 37
IAS Radiometry Algorithm 6.6a, “Outliers in a Regression”, Goddard Space Flight Center,
1996.
Reference 38
EOSAT Ground Segment Interface Control Document Volume 2D, G/SICD2040
Revision D, Nov, 1, 1996.
Reference 39
Gore, Robert, “Landsat 4/5 Attitude Jitter Processing Analysis”, Lockheed/Martin
Marietta, interoffice memo.
Reference 40
Gore, Robert, “Landsat 7 System ETM+ Process Algorithms (Attitude Jitter Processing)”,
contract No. NAS532633.
Version 3.2
5
7/9/98
IAS Geometric ATBD Version 3
2. Overview and Background Information
The Landsat 7 Enhanced Thematic Mapper Plus (ETM+) continues the Landsat program
of multispectral spaceborne Earth remote sensing satellites that began with the launch of
Landsat 1 in 1972. Landsat 7 will provide data consistent with the Landsat historical
record in the time period between the decommissioning of the stillfunctioning Landsat 5
and the introduction of successor instruments based on newer technology in the next
century. The basic sensor technology used in the ETM+ is similar to the Thematic
Mapper instruments flown in Landsats 4 and 5, and the Enhanced Thematic Mapper
built for Landsat 6, which suffered a launch failure. The geometric processing,
characterization, and calibration algorithms described in this document take into account
the new 15 meter panchromatic band (also present in Landsat 6) and the higher resolution
thermal band in adapting the processing techniques applied to Landsat 4/5 Thematic
Mapper data to the Landsat 7 ETM+.
The inclusion of an Image Assessment System as an integral part of the Landsat 7 ground
system illustrates the more systematic approach to instrument calibration and inorbit
monitoring and characterization necessitated by the more stringent calibration
requirements of the Landsat 7 ETM+ as compared to earlier Landsat missions. This is
especially true of the radiometric calibration requirements which include using partial and
fullaperture solar calibrations to monitor the stability and performance of the ETM+
detectors and onboard calibration lamps. The IAS also provides the platform for
systematic geometric performance monitoring, characterization and calibration as
described in the remainder of this document.
2.1 Experimental Objective
The objective of the Landsat 7 ETM+ mission is to provide high resolution (15 m
panchromatic, 30 meter multispectral, 60 meter thermal) imagery of Earth’s land areas
from near polar, sunsynchronous orbit. These data will continue and extend the
continuous data record collected by Landsats 15, provide higher spatial resolution in the
new panchromatic band, provide greater calibration accuracy to support new and
improved analysis applications, and provide a high resolution reference for the EOSAM1 MODIS, MISR, and ASTER instruments. The geometric algorithms described in
this document will be used by the Landsat 7 Image Assessment System and by the Level1 Product Generation System (LPGS) to ensure that the complex ETM+ internal
geometry is sufficiently modeled and characterized to permit the generation of Level 1
products which meet the Landsat 7 system geometric accuracy specifications.
2.2 Historical Perspective
The Landsat 7 ETM+ mission continues the evolutionary improvement of the Landsat
family of satellites which began with Landsats 1, 2, and 3 carrying the MultiSpectral
Version 3.2
6
7/9/98
IAS Geometric ATBD Version 3
Scanner (MSS), and builds on the Landsat 4 and 5 Thematic Mapper (TM) heritage. The
ETM+ adds a 15 meter GSD panchromatic band, similar to that used in the Enhanced
Thematic Mapper (ETM) instrument built for Landsat 6, improves the thermal band
spatial resolution from 120 meters (TM/ETM) to 60 meters (ETM+), and improves the
radiometric calibration accuracy from 10% (TM/ETM) to 5% (ETM+). The Landsat 7
system is also taking a more systematic approach to monitoring and measuring system
geometric performance in flight than was applied to Landsats 4 and 5, which were turned
over to commercial operations in October, 1985.
2.3 Instrument Characteristics
The fundamental ETM+ instrument geometric characteristics include:
• Nadir viewing +/ 7.5 degree instrument FOV;
• Bidirectional cross track scanning mirror;
• Along track scan line corrector mirrors;
• 8 bands, 3 resolutions, 2 focal planes;
• Even/odd detectors staggered on focal planes;
• Prelaunch ground calibrated nonlinear mirror profile;
• Scan angle monitor measures actual midscan and scan end times (as deviations
from nominal);
• High frequency (2125 Hz) angular displacement sensors to measure instrument
jitter.
2.4 Ancillary Input Data
The Landsat 7 ETM+ geometric characterization, calibration, and correction algorithms
are applied to the wideband (image plus supporting PCD and MSCD) data contained in
level 0R (raw) or 1R (radiometrically corrected) products. Some of these algorithms also
require additional ancillary input data sets. These include:
1. Precise ephemeris from the Flight Dynamics Facility  used to minimize
ephemeris error when performing sensor to spacecraft alignment calibration.
2. Ground control/reference images for geometric test sites  used in precision
correction, geodetic accuracy assessment, and geometric calibration algorithms.
3. Digital elevation data for geometric test sites  used in terrain correction and
geometric calibration.
4. Prelaunch ground calibration results including band/detector placement and
timing, scan mirror profiles, and attitude sensor data transfer functions (gyro and
ADS), to be used in the generation of the initial Calibration Parameter File.
5. Earth parameters including: static Earth model parameters (e.g., ellipsoid axes,
gravity constants) and dynamic Earth model parameters (e.g., polar wander
offsets, UT1UTC time corrections)  used in systematic model creation and
incorporated into the Calibration Parameter File.
Version 3.2
7
7/9/98
IAS Geometric ATBD Version 3
3. Algorithm Descriptions
This section presents the underlying theory and mathematical development of the IAS
geometric algorithms.
3.1 Theoretical Description
The supporting theoretical concepts and mathematics of the IAS geometric algorithms
are presented in the following subsections. Section 3.1.1 presents a review of the
Landsat 7 ETM+ viewing geometry to put the subsequent discussion in context. Sections
3.1.2 and 3.1.3 address the coordinate systems used by the algorithms and the
relationships between them, citing references where appropriate. Section 3.1.4 briefly
defines and discusses the various time systems used by the IAS algorithms. Section
3.1.5 presents the mathematical development of, and solution procedures for the Level 1
processing, geometric calibration, and geometric characterization algorithms. Section
3.1.6 discusses estimates of uncertainty (error analysis) associated with each of the
algorithms.
3.1.1 Landsat 7 ETM+ Viewing Geometry Overview
The ETM+ instrument detectors are aligned in parallel rows on two separate focal
planes: the primary focal plane, containing bands 14 and 8 (panchromatic), and the cold
focal plane containing bands 5, 6, and 7. The primary focal plane is illuminated by the
ETM+ scanning mirror, primary mirror, secondary mirror, and scan line corrector mirror.
In addition to these optical elements, the cold focal plane optical train includes the relay
folding mirror and the spherical relay mirror. This is depicted in Figure 1. The ETM+
scan mirror provides a nearly linear cross track scan motion which covers a 185 km wide
swath on the ground. The scan line corrector compensates for the forward motion of the
spacecraft and allows the scan mirror to produce usable data in both scan directions.
Within each focal plane, the rows of detectors from each band are separated in the alongscan (crosstrack) direction. The odd and even detectors from each band are also
separated slightly. The detector layout geometry is shown in Figure 3.1.12.. Samples
from the ETM+ detectors are packaged into minor frames as they are collected (in time
order) for downlink in the wideband data stream. Within each minor frame, the output
from the odd detectors from bands 15, and 7 and all band 8 detectors are held at the
beginning of the minor frame. The even detectors from bands 15, and 7 and all band 8
detectors are held at the midpoint of the minor frame. The even and odd detectors from
band 6 are held at the beginning of alternate minor frames. The bands are nominally
aligned during level 0R ground processing by delaying the even and odd detector
samples from each band to account for the alongscan motion needed to view the same
target point. These delays are implemented as two sets of fixed integer offsets for the
even and odd detectors of each bandone set for forward scans and one set for reverse
Version 3.2
8
7/9/98
IAS Geometric ATBD Version 3
scans. Two sets of offset delays are necessary since the band/detector sampling sequence
with respect to a ground target is inverted for reverse scans. Level 0R ground processing
also reverses the sample time sequence for reverse scans to provide nominal spatial
alignment between forward and reverse scans.
1
16
16
1
1
1
7
16
Cooled Focal
5 Plane Detector
16
8
16
6
Orientation
16
1
1
Orientation
1
1
1
32
4
3 Primary Focal
2 Plane Detector
8
Primary Mirror
Secondary Mirror
Scan Mirror
Relay
Folding
Mirror
Spherical
Relay
Mirror
Scan Line
Corrector
Figure 1 ETM+ Geometry Overview
The instrument's 15degree field of view is swept over the ETM+ focal planes by the
scan mirror. Each scan cycle, consisting of a forward and reverse scan, has a period of
Version 3.2
9
7/9/98
IAS Geometric ATBD Version 3
142.245 milliseconds. The scan mirror period for a single forward or reverse scan is
nominally 71.462 milliseconds. Of this time, approximately 60.743 milliseconds is
devoted to the Earth viewing portion of the scan with 10.719 milliseconds devoted to the
collection of calibration data and mirror turnaround. Within each scan, detector samples
(minor frames) are collected every 9.611 microseconds (with two samples per minor
frame for the 15 meter resolution panchromatic band). More detailed information on the
ETM+ instrument's construction, operation, and preflight testing is provided in
Reference 5.
Primary Focal Plane
Band 8
Band 1
Band 2
Cold Focal Plane
Band 3
Band 4
Band 7
Band 5
Band 6
Optical
Axis
10.4
25 IFOV
25 IFOV
25 IFOV
25 IFOV
45 IFOV
26 IFOV
34.75 IFOV
Figure 2 ETM+ Detector Array Ground Projection
The 15degree instrument field of view sweeps out a ground swath approximately 185
kilometers wide. This swath is sampled 6,320 (nominally) times by the ETM+ 30meter
resolution bands. Since sixteen detectors from each 30meter band are sampled in each
data frame, the nominal scan width is 480 meters. The scan line corrector compensates
for spacecraft motion during the scan to prevent excessive (more than 1 pixel) overlap or
underlap at the scan edges. Despite this, the varying effects of Earth rotation and
spacecraft altitude (both of which are functions of position in orbit) will lead to variations
in the scantoscan gap observed.
3.1.2 Coordinate Systems
The coordinate systems used by the Landsat 7 system are described in detail in
Reference 2. There are ten coordinate systems of interest for the Landsat 7 Image
Assessment System geometric algorithms. These coordinate systems are referred to
frequently in the remainder of this document and are briefly defined here to provide a
context for the subsequent discussion. They are presented in the logical order in which
they would be used to transform a detector and sample time into a ground position.
1. ETM+ Array Reference (Focal Plane) Coordinate System
The focal plane coordinate system is used to define the band and detector offsets from
the instrument optical axis which are ultimately used to generate the image space viewing
vectors for individual detector samples. It is defined so that the Z axis is along the optical
axis and is positive toward the scan mirror. The origin is where the optical axis intersects
the focal planes (see Figure 3). The X axis is parallel to the scan mirror's axis of rotation
Version 3.2
10
7/9/98
IAS Geometric ATBD Version 3
and is in the alongtrack direction after a reflection off the scan mirror with the positive
direction toward detector 1. The Y axis is in the focal plane's alongscan direction with the
positive direction toward band 8. This definition deviates from the description in the
Landsat 7 Program Coordinate System Standard document in that the origin is placed at
the focal plane/optical axis intersection rather than at an arbitrary reference detector
location.
8
1
32
2
16
3
16
4
16
7
16
5
16
6
16
8
Y
1
1
1
1
1
1
1
1
X
Figure 3 ETM+ Focal Plane Coordinate System
2. ETM+ Sensor Coordinate System
The scan mirror and scan line corrector models are defined in the sensor coordinate
system. This is the coordinate system in which an image space vector, representing
the line of sight from the center of a detector, and scan mirror and scan line corrector
angles are converted to an object space viewing vector. The Z axis corresponds to the Z
axis of the focal plane coordinate system (optical axis) after reflection off the scan mirror
and is positive outward (toward the Earth). The X axis is parallel to the scan mirror axis
of rotation and is positive in the alongtrack direction. The Y axis completes a right
handed coordinate system. Scan mirror angles are rotations about the sensor X axis and
are measured in the sensor YZ plane with the Z axis corresponding to a zero scan angle
and positive angles toward the Y axis (west in descending mode). Scan line corrector
angles are rotations about the sensor Y axis and are measured in the sensor XZ plane
with the Z axis corresponding to a zero scan line corrector angle and positive angles
toward the X axis (the direction of flight). The sensor coordinate system is depicted in
Figure 4. In this coordinate system, the scan line corrector performs a negative rotation
(from a positive angle to a negative angle) on every scan while the scan mirror performs a
positive to negative rotation for forward scans and a negative to positive rotation for
reverse scans.
Version 3.2
11
7/9/98
IAS Geometric ATBD Version 3
Y
X
Z
Figure 4 ETM+ Sensor Coordinate System
3. Navigation Reference Coordinate System
The navigation reference frame is the bodyfixed coordinate system used for spacecraft
attitude determination and control. The coordinate axes are defined by the spacecraft
attitude control system (ACS) which attempts to keep the navigation reference frame
aligned with the orbital coordinate system so that the ETM+ optical axis is always
pointing toward the center of the Earth. It is the orientation of this coordinate system
relative to the inertial coordinate system that is captured in the spacecraft attitude data.
4. Inertial Measurement Unit (IMU) Coordinate System
The spacecraft orientation data provided by the gyros in the inertial measurement unit are
referenced to the IMU coordinate system. This coordinate system is nominally aligned
with the navigation reference coordinate system. The actual alignment of the IMU with
respect to the navigation reference will be measured preflight as part of the attitude
control system calibration. This alignment transformation is used by the IAS to convert
the gyro data contained in the Landsat 7 Payload Correction Data (PCD) to the navigation
reference coordinate system for blending with the ACS quaternions and the Angular
Displacement Assembly (ADA) jitter measurements.
5. Angular Displacement Assembly (ADA) Coordinate System
The high frequency angular displacements of the sensor line of sight (jitter) measured by
the ADA are referenced to the ADA coordinate system. This reference frame is at a
nominal 20degree angle with respect to the sensor coordinate system, corresponding to a
20degree rotation about the ADA X axis, which is nominally coincident with the sensor
Version 3.2
12
7/9/98
IAS Geometric ATBD Version 3
X axis. The ADA output recorded in the Landsat 7 PCD is converted from the ADA
coordinate system to the sensor coordinate system, using the ADA alignment matrix, and
then to the navigation reference coordinate system (using the sensor alignment matrix) so
that the high frequency jitter data can be combined with the low frequency gyro data to
determine the actual ETM+ sensor lineofsight pointing as a function of time.
6. Orbital Coordinate System
The orbital coordinate system is centered on the satellite, and its orientation is based on
the spacecraft position in inertial space (see Figure 5). The origin is the spacecraft center
of mass, with the Z axis pointing from the spacecraft center of mass to the Earth center of
mass. The Y axis is the normalized cross product of the Z axis and the instantaneous
(inertial) velocity vector, and corresponds to the negative of the instantaneous angular
momentum vector direction. The X axis is the cross product of the Y and Z axes.
Z eci
Spacecraft
Position
Yorb
Spacecraft Orbit
Earth's Axis of Rotation
X
X
orb
Zorb
Yeci
Equator
eci
To Vernal
Equinox
Figure 5 Orbital Coordinate System
Version 3.2
13
7/9/98
IAS Geometric ATBD Version 3
7. Earth Centered Inertial (ECI) Coordinate System
The ECI coordinate system is spacefixed with its origin at the Earth's center of mass (see
Figure 6). The Z axis corresponds to the mean north celestial pole of epoch J2000.0. The
X axis is based on the mean vernal equinox of epoch J2000.0. The Y axis is the cross
product of the Z and X axes. This coordinate system is described in detail in Reference 7
and in Reference 2. Data in the ECI coordinate system will be present in the Landsat 7
ETM+ Level 0R product in the form of ephemeris and attitude data contained in the
spacecraft PCD.
Although the IAS uses J2000, there is some ambiguity in the existing Landsat 7 system
documents which discuss the inertial coordinate system. Reference 2 states that the
ephemeris and attitude data provided to and computed by the spacecraft onboard
computer are referenced to the ECI True of Date (ECITOD) coordinate system (which is
based on the celestial pole and vernal equinox of date rather than at epoch J2000.0). This
appears to contradict Reference 1 which states that the ephemeris and attitude data
contained in the Landsat 7 PCD (i.e., from the spacecraft onboard computer) is
referenced to J2000.0. The relationship between these two inertial coordinate systems
consists of the slow variation in orientation due to nutation and precession. This is
described in Reference 2 and 7.
Z
Earth's Axis of Rotation
Y
Equator
X
To Vernal
Equinox
Version 3.2
14
7/9/98
IAS Geometric ATBD Version 3
Figure 6 Earth Centered Inertial (ECI) Coordinate System
8. Earth Centered Rotating (ECR) Coordinate System
The ECR coordinate system is Earthfixed with its origin at the center of mass of the
Earth (see Figure 7). It corresponds to the Conventional Terrestrial System defined by
the Bureau International de l’Heure (BIH) which is the same as the U. S. Department of
Defense World Geodetic System 1984 (WGS84) geocentric reference system. This
coordinate system is described in Reference 7.
Z
Greenwich Meridian
Earth's Axis of Rotation
Y
Equator
X
Figure 7 Earth Centered Rotating (ECR) Coordinate System
9. Geodetic Coordinate System
The geodetic coordinate system is based on the WGS84 reference frame with coordinates
expressed in latitude, longitude, and height above the reference Earth ellipsoid (see Figure
8). No ellipsoid is required by the definition of the ECR coordinate system, but the
geodetic coordinate system depends on the selection of an Earth ellipsoid. Latitude and
longitude are defined as the angle between the ellipsoid normal and its projection onto the
equator, and the angle between the local meridian and the Greenwich meridian,
Version 3.2
15
7/9/98
IAS Geometric ATBD Version 3
respectively. The scene center and scene corner coordinates in the Level 0R product
metadata are expressed in the geodetic coordinate system.
Z
Greenwich Meridian
Height
Ellipsoid Normal
Latitude
Equator
Y
Longitude
X
Figure 8 Geodetic Coordinate System
10. Map Projection Coordinate System
Level 1 products are generated with respect to a map projection coordinate system, such
as Universal Transverse Mercator, which provides a mapping from latitude and longitude
to a plane coordinate system which is an approximation to a Cartesian coordinate system
for a portion of the Earth’s surface. It is used for convenience as a method of providing
the digital image data in an Earthreferenced grid that is compatible with other ground
referenced data sets. Although the map projection coordinate system is only an
approximation to a true local Cartesian coordinate system at the Earth’s surface, the
mathematical relationship between the map projection and geodetic coordinate systems is
precisely and unambiguously defined; the map projections the IAS uses are described in
Reference 8.
3.1.3 Coordinate Transformations
There are nine transformations between the ten coordinate systems used by the IAS
geometric algorithms. These transformations are referred to frequently in the remainder of
this document and are defined here. They are presented in the logical order in which a
detector and sample number would be transformed into a ground position.
Version 3.2
16
7/9/98
IAS Geometric ATBD Version 3
1. Focal Plane to Sensor
The relationship between the focal plane and sensor coordinate systems is a rapidly
varying function of time, incorporating the scan mirror and scan line corrector models.
The focal plane location of each detector from each band can be specified by a set of field
angles in the focal plane coordinate system, as described in section 3.1.5.1.1. These are
transformed to sensor coordinates based on the sampling time (sample/minor frame
number from scan start) which is used to compute mirror and scan line corrector angles
based on the models described in sections 3.1.5.1.2 and 3.1.5.1.3.
2. Sensor to Navigation Reference
The relationship between the sensor and navigation reference coordinate systems is
described by the ETM+ instrument alignment matrix. The transformation from sensor
coordinates to navigation reference coordinates is described in section 4.21 of Reference
2. It includes a threedimensional rotation, implemented as a matrix multiplication,
and an offset to account for the distance between the ACS reference and the instrument
scan mirror. The transformation matrix will initially be defined to be fixed (nontime
varying) with improved estimates being provided postlaunch. Subsequent analysis may
detect repeatable variations with time that can be effectively modeled, making this a
(slowly) timevarying transformation. The nominal rotation matrix is the identity matrix.
3. IMU to Navigation Reference
The IMU coordinate system is related to the navigation reference coordinate system by
the IMU alignment matrix which captures the orientation of the IMU axes with respect to
the navigation base. This transformation is applied to gyro data prior to its integration
with the ADA data and the ACS quaternions. The IMU alignment will be measured preflight and is nominally the identity matrix. This transformation is defined in section 4.22
of Reference 2.
4. ADA to Sensor
The angular displacement assembly (ADA) is nominally aligned with the sensor
coordinate system with a 20degree rotation about the X axis. The actual alignment will
be measured prelaunch and will be used to rotate the ADA jitter observations into the
sensor coordinate system, from which they can be further rotated into the navigation
reference system using the sensor to navigation reference alignment matrix.
5. Navigation Reference to Orbital
The relationship between the navigation reference and orbital coordinate systems is
defined by the spacecraft attitude. This transformation is a threedimensional rotation
matrix with the components of the rotation matrix being functions of the spacecraft roll,
Version 3.2
17
7/9/98
IAS Geometric ATBD Version 3
pitch, and yaw attitude angles. The nature of the functions of roll, pitch, and yaw
depends on the exact definition of these angles ( i.e., how they are generated by the
attitude control system). In the initial model, it is assumed that the rotations are
performed in the order rollpitchyaw as shown in section 4.17 of Reference 2. Since the
spacecraft attitude is constantly changing, this transformation is time varying. The
nominal rotation matrix is the identity matrix since the ACS strives to maintain
the ETM+ instrument pointing toward the center of the Earth.
6. Orbital to ECI
The relationship between the orbital and ECI coordinate systems is based on the
spacecraft's instantaneous ECI position and velocity vectors. The rotation matrix to
convert from orbital to ECI can be constructed by forming the orbital coordinate system
axes in ECI coordinates:
p
v
Teci/orb 
spacecraft position vector in ECI
spacecraft velocity vector in ECI
rotation matrix from orbital to ECI
b3 = p / p
b2 = (b3 x v) /  b3 x v
b1 = b2 x b3
(nadir vector direction)
(negative of angular momentum vector direction)
(circular velocity vector direction)
Teci/orb = [ b1 b2 b3 ]
7. ECI to ECR
The transformation from ECI to ECR coordinates is a timevarying rotation due primarily
to Earth rotation, but also containing more slowly varying terms for precession,
astronomic nutation, and polar wander. The ECI to ECR rotation matrix can be expressed
as a composite of these transformations:
Tecr/eci = A B C D
A = Polar Motion
B = Sidereal Time
C = Astronomic Nutation
D = Precession
Each of these transformation terms is defined in Reference 2 and is described in detail in
Reference 7.
Version 3.2
18
7/9/98
IAS Geometric ATBD Version 3
8. ECR to Geodetic
The relationship between ECR and geodetic coordinates can be expressed simply in its
direct form:
e2 = 1  b2 / a2
N = a / (1  e2 sin2(lat))1/2
X = (N + h) cos(lat) cos(lon)
Y = (N + h) cos(lat) sin(lon)
Z = (N (1e2) + h) sin(lat)
where:
X, Y, Z
lat, lon, h
N
e2
a, b

ECR coordinates
Geodetic coordinates
Ellipsoid radius of curvature in the prime vertical
Ellipsoid eccentricity squared
Ellipsoid semimajor and semiminor axes
The closedform solution for the general inverse problem (which is the problem of interest
here) involves the solution of a quadratic equation and is not typically used in practice.
Instead, an iterative solution is used for latitude and height for points that do not lie on
the ellipsoid surface. A procedure for performing this iterative calculation is described in
section 6.4.3.3 of Reference 4.
9. Geodetic to Map Projection
The transformation from geodetic coordinates to the output map projection depends on
the type of projection selected. The mathematics for the forward and inverse
transformations for the Universal Transverse Mercator (UTM), Lambert Conformal
Conic, Transverse Mercator, Oblique Mercator, Polyconic, Polar Stereo Graphic, and the
Space Oblique Mercator (SOM) are given in Reference 8. Further details of the SOM
mathematical development are presented in Reference 17.
3.1.4 Time Systems
Four time systems are of primary interest for the IAS geometric algorithms. These
include: International Atomic Time (Temps Atomique International or TAI), Universal
Time  Coordinated (UTC), Universal Time corrected for polar motion (UT1), and
Spacecraft Time (the readout of the spacecraft clock, nominally coincident with UTC).
Spacecraft Time is the time system applicable to the spacecraft time codes found in the
Level 0R PCD and MSCD. UTC (which the spacecraft clock aspires to) is the standard
reference for civil timekeeping. UTC is adjusted periodically by whole leap seconds to
keep it within 0.9 seconds of UT1. UT1 is based on the actual rotation of the Earth and
is needed to provide the transformation from stellarreferenced inertial coordinates (ECI)
Version 3.2
19
7/9/98
IAS Geometric ATBD Version 3
to terrestrialreferenced Earthfixed coordinates (ECR). TAI provides a uniform,
continuous time stream which is not interrupted by leap seconds or other periodic
adjustments. It provides a consistent reference for resolving ambiguities arising from the
insertion of leap seconds into UTC (which can lead to consecutive seconds with the same
UTC time). These and a variety of other time systems and their relationships are
described in Reference 4. The significance of each of these four time systems with
respect to the IAS geometric algorithms is described below.
1. Spacecraft Time
The Landsat 7 Mission Operations Center will attempt to maintain the onboard
spacecraft clock time to within 145 milliseconds of UTC (Reference 5, section 3.7.1.3.16)
by making clock updates nominally once per day during periods of ETM+ inactivity.
Additionally, the spacecraft PCD includes a quadratic clock correction model which can
be used to correct the spacecraft clock readout to within 15 milliseconds of UTC,
according to Reference 1 . The clock correction algorithm is:
dt = ts/c  t update
t UTC = t s/c + C0 + C1 dt + 0.5 C 2 dt 2
where:
dt
t UTC
t s/c
t update
C0
C1
C2
= spacecraft clock time since last update
= corrected spacecraft time (UTC +/ 15 milliseconds)
= spacecraft clock time
= spacecraft clock time of last ground commanded clock update
= clock correction bias term
= clock drift rate
= clock drift acceleration
The Landsat 7 onboard computer is using corrected spacecraft clock data to interpolate
the ephemeris data (which is inserted into the PCD) from the uplinked ephemeris (which
is referenced to UTC). The residual clock error of 15 milliseconds (3 sigma) leads to an
additional ephemeris position uncertainty of approximately 114 meters (3 sigma) which is
included in the overall geodetic accuracy error budget. Since the time code readouts in
both the PCD and MSCD are uncompensated, the clock correction algorithm must be
applied on the ground to get the Level 0R time reference to within 15 milliseconds of
UTC.
Note that if the onboard computer were not using the clock drift correction model for
ephemeris interpolation the ground processing problem would become more complicated.
Since the spacecraft time is used as a UTC index to interpolate the ephemeris points, the
times associated with these ephemeris points are true UTC  they just do not correspond
to the actual (corrected) UTC spacecraft time. In this case, the PCD time code errors of
Version 3.2
20
7/9/98
IAS Geometric ATBD Version 3
up to 145 milliseconds would cause a temporal misregistration of the interpolated
ephemeris and the associated orbital orientation reference coordinate system, relative to
the spacecraft clock. In making the spacecraft clock corrections, the times associated with
the ephemeris points would not be changed in recognition of the fact that they are
interpolation indices against a UTC time reference.
2. UTC
As mentioned above, UTC is maintained within 0.9 seconds of UT1 by the occasional
insertion of leap seconds. It is assumed that the insertion of leap seconds into the
spacecraft clock will be performed on the appropriate dates as part of the daily clock
update, and that this clock update will be coordinated with the ephemeris updates
provided by the FDF, so that the spacecraft clock UTC epoch is the same as the FDF
provided ephemeris UTC epoch. Although for any given scene or subinterval UTC
provides a uniform time reference, it will be beneficial for the IAS to have access to a table
of leap seconds relating UTC to TAI to support timerelated problem tracking and
resolution. This information will be obtained from the National Earth Orientation Service
(NEOS) web site at http://maia.usno.navy.mil.
3. UT1
UT1 represents time with respect to the actual rotation of the Earth, and is used by the
IAS algorithms which transform inertial ECI coordinates or lines of sight to Earthfixed
ECR coordinates. Failure to account for the difference between UT1 and UTC in these
algorithms can lead to ground position errors as large as 400 meters at the equator
(assuming the maximum 0.9 second UT1  UTC difference). The UT1  UTC correction
typically varies at the rate of approximately 2 milliseconds per day, corresponding to an
Earth rotation error of about 1 meter. Thus, UT1  UTC corrections should be
interpolated or predicted to the actual image acquisition time to avoid introducing errors
of this magnitude. This information will be obtained from the National Earth Orientation
Service (NEOS) web site at http://maia.usno.navy.mil.
4. TAI
Although none of the IAS algorithms are being designed to use TAI in normal operations
it is included here for completeness. As mentioned above, there is the potential for
system anomalies to be introduced as a result of timing errors, particularly in the vicinity
of UTC leap seconds. The ability to use TAI as a standard reference which can be related
to UTC and spacecraft time using a leap second file, will assist IAS operations staff in
anomaly resolution. For this reason, although there are no IAS software requirements for
a leap second reference it is strongly recommended that access be provided to the ECS
SDP Toolkit leap second file.
Version 3.2
21
7/9/98
IAS Geometric ATBD Version 3
3.1.5 Mathematical Description of Algorithms
3.1.5.1 ETM+ Instrument Model
3.1.5.1.1 Detector/Focal Plane Geometry
The location of each detector for each band is described in the ETM+ array reference
coordinate system by using a pair of angles, one in the alongscan, or coordinate Y,
direction and one in the acrossscan, or coordinate X, direction. In implementing the IAS
level 1 processing algorithm, each of these angles is separated into two parts: the ideal
band/detector locations used to model the relationship between the sensor image space
and ground object space, and the subpixel offsets unique to each detector which are
applied during the image resampling process. This approach makes it possible to
construct a single analytical model for each band from which nominal detector projections
can be interpolated and then refined during resampling, rather than performing the image
to ground projection computation for each detector. The capability to rigorously project
each detector, including all subpixel detector placement and timing offsets, has been
developed to support geometric characterization and calibration.
Alongscan angles are modeled as the sum of the band odd detector angle from the optical
axis, common to all detectors for one band, and the subpixel offset and time delay unique
to each detector within the band. Note that the IAS algorithms assume that the angular
displacement between even and odd detectors in each band is exactly compensated for by
the time delay used to nominally align the even and odd detectors. Band center location
angles, relative to the sensor optical axis, are stored in the Calibration Parameter File. The
odd detector offset from band center (1.25 IFOVs for the 30 meter bands) is subtracted
from the band center angle to yield the nominal detector locations. The alongscan angle is
computed as:
along_angle = bandoff_along[BAND]  odd_det_off[BAND]
where:
along_angle
bandoff_along
BAND
odd_det_off
= nominal alongscan location of detector
= alongscan band center offset (indexed by band)
= band number (1..8)
= odd detector offset from band center (indexed by band)
The individual detector offsets from their ideal locations are stored in the Calibration
Parameter File in the form of detector delays which incorporate both sample timing
offsets and detector alongscan positioning errors. These subpixel corrections are
applied during the image resampling process.
Version 3.2
22
7/9/98
IAS Geometric ATBD Version 3
Similarly, the acrossscan angles are separated into the nominal detector location
component, used for line of sight construction and projection, and the subpixel detector
placement corrections unique to each detector, applied during resampling. The nominal
detector location is computed as:
cross_angle = bandoff_cross[BAND] + ((ndets[BAND]+1)/2  N)*IFOV[BAND]
where:
cross_angle
bandoff_cross
BAND
ndets
N
IFOV
= nominal acrossscan location of detector N
= acrossscan band center offset (indexed by band)
= band number (1..8)
= number of detectors in band (indexed by band)
= detector number (1..ndets)
= detector angular field of view (indexed by band)
Using these nominal detector locations in the line of sight projection computations makes
it possible to rigorously project the first and last detectors in each band and to interpolate
the image space to ground space transformation for the interior detectors. This reduces
the computational load in level 1 processing where the unique subpixel detector
placement corrections are applied by the resampling algorithm. To precisely project an
individual detector to ground space, the subpixel alongscan and acrossscan corrections
are added to the nominal detector location.
3.1.5.1.2 Scan Mirror Model
The ETM+ scan mirror rotates about the sensor X axis providing nearly linear motion in
the alongscan direction in both the forward and reverse directions. The scan mirror’s
departure from linear motion is characterized prelaunch using a set of fifthorder
polynomials which model the repeatable acceleration and deceleration deviations from
linearity for forward and reverse scans. This fifthorder polynomial is then adjusted for
high frequency roll jitter components using ADS data from the PCD. Fifthorder mirror
polynomials are also used to model the acrossscan mirror deviation for forward and
reverse scans. Thus, a total of four fifthorder polynomials characterize the scan mirror’s
repeatable deviation from nominal linear scanning: 1) alongscan deviation for forward
scans, 2) acrossscan deviation for forward scans, 3) alongscan deviation for reverse
scans, and 4) acrossscan deviation for reverse scans. This set of four polynomials
constitutes the nonlinear portion of the scan mirror profile.
The ETM+ instrument contains two sets of Scan Mirror Electronics (SME) used to
control the mirror motion. Either SME can be used to drive the mirror in either of its
operational modes, Scan Angle Monitor (SAM) mode or Bumper mode. For each mode,
each set of SME exhibits a characteristic set of mirror profile polynomials as well as scan
Version 3.2
23
7/9/98
IAS Geometric ATBD Version 3
start and scan stop mirror angles. This leads to a total of four sets of mirror profile
polynomials: 1) SME1 SAM mode, 2) SME2 SAM mode, 3) SME1 Bumper mode, and
4) SME2 Bumper mode. The appropriate mirror profile coefficients and scan start/stop
angles must be selected for a particular data set based on the SME and mode settings in
the PCD.
Within each scan, the scan mirror assembly measures the deviation from the nominal time
of scan from scan start to midscan (first half scan error), and from midscan to scan stop
(second half scan error) and records this information in the scan line data minor frames for
the following scan. These scan time errors are measured as departures from nominal, in
counts, with each count being 0.18845 microseconds.
The first half and second half scan errors are used to correct the nominal scan profile by
adjusting the fifthorder alongscan polynomial to compensate for the actual midscan and
end of scan times. This correction is applied by rescaling the polynomial coefficients to
the actual scan time and by computing correction coefficients which are added to the
linear and quadratic terms in the fifthorder polynomial. The calculation procedure is
given in section 3.1.5.3.3.
3.1.5.1.3 Scan Line Corrector Mirror Model
The Scan Line Corrector (SLC) mirror rotates the ETM+ line of sight about the sensor Y
axis to compensate for the spacecraft’s along track motion during the scanning interval.
Although the IAS algorithms support fifthorder correction polynomials for the SLC
profile, in practice, no significant SLC nonlinearity has been observed for Landsats 4, 5,
and 6. Refer to section 3.1.5.3.4 for more information on the Scan Line Corrector Mirror
Model.
3.1.5.2
Landsat 7 Platform Model
3.1.5.2.1 Sensor Pointing (Attitude and Jitter) Model
Sensor line of sight vectors constructed using the ETM+ instrument model must be
related to object (ground) space by determining the orientation and position of the sensor
coordinate system with respect to object space (ECI coordinates) as a function of time.
The sensor orientation is determined using the sensor pointing model which includes the
effects of sensor to spacecraft (navigation reference base) alignment, spacecraft attitude as
determined by the onboard computer using the inertial measurement unit gyros and star
sightings, and the high frequency attitude variations (jitter) measured by the angular
displacement assembly. The inputs to this model are the attitude quaternions, gyro drift
estimates, raw gyro outputs, and ADS outputs contained in the Payload Correction Data,
and the sensor alignment matrix measured prelaunch and updated in flight by the Image
Version 3.2
24
7/9/98
IAS Geometric ATBD Version 3
Assessment System using the sensor alignment calibration algorithm described in section
3.1.5.4.1.
The Landsat 7 attitude control system outputs estimates of the spacecraft attitude every
PCD major frame (4.096 seconds). These estimates are provided in the form of attitude
quaternions which describe the orientation of the spacecraft navigation reference system
with respect to the ECI (J2000) coordinate system. These estimates are derived from the
IMU gyro data combined with periodic star transit observations from the Celestial Sensor
Assembly. As a byproduct of the integration of the gyro and star tracker observations
the ACS also computes estimates of the gyro drift rates. A gyro drift rate estimate is also
included in each PCD major frame, but only changes after a star sighting.
The raw outputs from each of the three gyro axes are recorded in the PCD every 64
milliseconds. In Landsats 4, 5, and 6 there was a 28 millisecond delay between the gyro
sampling time and the PCD major frame start time ( i.e., gyro sampling was indexed from
PCD start time  28 milliseconds). From Reference 1 , it appears that this offset does not
apply for Landsat 7. The 15.625 Hertz gyro sampling rate is not sufficient to capture the
high frequency spacecraft angular dynamic displacement (jitter). This sampling frequency
will capture attitude variations up to the Nyquist frequency of 7.8125 Hz but will alias
the higher frequencies to the extent the gyros respond to these higher frequencies. This
response is expressed in the IMU transfer functions. For Landsats 4 and 5, the transfer
functions for each IMU axis were taken to be the same.
The high frequency jitter information is measured by the angular displacement sensors
(ADS) which are sampled at 2 millisecond intervals and are sensitive to frequency
components between approximately 2 Hz and 125 Hz. Like the gyros, the ADS samples
are offset from the PCD major frame start time. Unlike the gyros, the ADS axes are
sampled at different times. The ADS samples begin at the PCD major frame time plus
375 microseconds for the X axis, plus 875 microseconds for the Y axis, and plus 1375
microseconds for the Z axis. The 500 Hz ADS sampling frequency could capture jitter
frequencies up to 250 Hz, but the output of the ADS is bandlimited to the maximum
expected frequency of 125 Hz by a presampling filter. Historically, the three ADS axes
have had different transfer functions. For Landsats 4 and 5, the ADS transfer functions
were not separated from the presampling filter but were provided as a net response. For
Landsat 6, the presampling filter transfer function was provided separately.
The gyro and ADS transfer functions are used to blend the gyro/attitude and ADS/jitter
data in the frequency region where their pass bands overlap from 2 to 7 Hz. The IAS
uses the classical approach defined by Sehn and Miller (see Reference 10 ) to construct
the crossover filters used to combine the gyro and ADS data. This is described in more
detail in section 3.1.5.3.3.1.
Version 3.2
25
7/9/98
IAS Geometric ATBD Version 3
3.1.5.2.2 Platform Position and Velocity
The spacecraft state vectors contained in the Payload Correction Data stream provide
position and velocity information, in the Earth Centered Inertial of epoch J2000
coordinate system, every PCD major frame (unlike earlier Landsats which only provided
new state vectors every other major frame). These state vectors are interpolated onboard
using the uplinked predicted ephemeris, and may contain significant errors. Gross
blunders (i.e., transmission errors) are detected and removed during PCD preprocessing.
The systematic errors which accumulate in the predicted ephemeris are estimated and
removed as part of the Precision Correction algorithm described in section 3.1.5.3.8.
3.1.5.3
Level 1 Processing Algorithms
The diagrams that follow describe the highlevel processing flows for the IAS level1
processing algorithms. Figure 3.1.51 describes the process involved in initialization and
creation of the Landsat 7 geometry model. Figure 3.1.52 shows the process of creating a
geometric correction grid and the application of that grid in the resampling process.
Figure 3.1.53 describes the process of refining the Landsat 7 geometry model with
ground control, resulting in a precision geometry model. Figure 3.1.54 again describes the
creation of a geometric correction grid (this time precision), and resampling with terrain
correction.
Detailed algorithms for each of main process boxes in these diagrams are given in the
sections that follow.
Version 3.2
26
7/9/98
IAS Geometric ATBD Version 3
FDF
Ephemeris
PCD Quality
Statistics
Cal
Parameter
File
Alignment
Mirror Profiles
Transfer Functions
Earth Constants
Instrument Constants
Process PCD
Compare I and Q
Detect Outliers
Compute Statistics
Correct S/C Clock
PCD
Validated PCD
Create Satellite Model
Process Attitude Data
Process Ephemeris Data
Construct Mirror Model
Compute and Validate
Scene Corners
Clock Correction
Coefficients
Process MSCD
Compare I and Q
Detect Outliers
Compute Statistics
Correct S/C Clock
MSCD
Satellite
Model
(S)
Validated MSCD
0R
Metadata
MSCD Quality
Statistics
Attitude Stats
Corner Fit
Model Initialization
9/19/96
Figure 3.1.5 1 Model Initialization
Processing
Parameters
Projection, Rotation,
Resampling Options,
Pixel Size, Framing Options Bands to Process
Generate Correction Grid
Compute Scene Framing
Grid Input Space
Relate Input Space to
Output Space
Processing
Parameters
1Rc
Image
(18)
Cal
Parameter
File
MTFC Coeffs
Det. Offsets
Resampling
(18)
Correction
Grid
(S)
Compute Scan Gap
Extend Image
Apply Detector Delays
Resample/Interpolate (2D)
1Gs
Image
(18)
Pixel/Line/Band
Lat/Long
Satellite
Model
(S)
Scan Gap
Statistics
Generate & Project LOS
(Call Model)
Compute Time and Angles
Construct LOS Vector
Intersect LOS with Earth
Rectification and Resampling
9/19/96
Figure 3.1.5 2 Rectification and Resampling
Version 3.2
27
7/9/98
IAS Geometric ATBD Version 3
GCP Library
GCP Correlation
(LAS Tools)
1Gs
Image
(8)
Cross, Phase,
Edge
Processing
Parameters
Ground
Control
Points
Inverse Mapping Lat,Lon,H
(18)
Correction
Grid
(S)
Fit Results
Apriori Weights
Parameterization
Options
Compute Predicted Time,Lat’,Lon’
Px’,Py’,Pz’
Lat/Long
Vx’,Vy’,Vz’
Compute Time
Compute State
Ephemeris
Precision Correction Updates
Update Model
Compute Ephemeris & Attitude
Attitude Updates
Updates
Detect GCP Outliers
Satellite
Model
(P)
Add Ephemeris
Updates to Model
Add Attitude
Updates to Model
Satellite
Model
(S)
Precision Correction Solution
9/19/96
Figure 3.1.5 3 Precision Correction Solution
Processing
Parameters
Projection, Rotation,
Resampling Options,
Pixel Size, Framing Options Bands to Process
Generate Correction Grid
Compute Scene Framing
Grid Input Space
Relate Input Space to
Output Space
Processing
Parameters
Pixel/Line/Band
1Rc
Image
(18)
Cal
Parameter
File
MTFC Coeffs
Det. Offsets
Resampling
(18)
Correction
Grid
(P)
Compute Scan Gap
Extend Image
Apply Detector Delays
Resample/Interpolate (2D)
1Gt
Image
(18)
Pixel/Line
Lat/Long
Pixel Offset
Terrain Correction
Generate & Project LOS
(Call Model)
Locate Nadir Per Scan
Compute Table of Offsets
Interpolate H from P/L
Compute Pixel OffNadir
Look Up Offset
Compute Time and Angles
Construct LOS Vector
Intersect LOS with Earth
Satellite
Model
(P)
DTM
Precision/Terrain Correct Image
9/19/96
Figure 3.1.5 4 Precision/Terrain Correction
Version 3.2
28
7/9/98
IAS Geometric ATBD Version 3
3.1.5.3.1 PCD Processing
The PCD contains information on the state of the satellite and sensor. The PCD
information is transmitted in both the InPhase (I) and Quadrature (Q) channels of the
wideband data. The two PCD data streams are converted to engineering units and stored
in the HDF format for the Level 0R product. The I channel is used as the default PCD
for processing. The PCD is recorded for the entire subinterval or full interval. The PCD
information required for geometric processing includes: the ephemeris, the spacecraft
attitude, the gyro data, the gyro drift rate data, the ADS data, the SME flag, the gyro unit
flag, and the spacecraft clock corrections. The ephemeris and spacecraft attitude are
updated every major frame, the gyro data is updated 64 times per axis in a major frame,
and the ADS is updated 2048 times per axis in a major frame. There are 4 PCD major
frames in a PCD cycle. Also included in the PCD are housekeeping values.
Validate Ephemeris Data
The PCD is transmitted in both the I and Q channels of the telemetry data. Under normal
conditions, the I and Q channels should be identical. Therefore, a comparison between
the ephemeris values contained in the I and Q channels is used as one of the validation
tests. Any differences detected are flagged and reported in the processing history file.
The I channel is assumed to be correct unless the other validation checks determine it to
be in error.
The Landsat 7 orbit's semimajor axis, inclination and angular momentum calculated from
the ephemeris data should not deviate substantially from their nominal values. A second
validation compares these values to their nominal values and any large deviations are
flagged. If a large deviation is detected, the Q channel's values are checked. If the Q
channel values are correct, they are used instead of the I channel values. The ephemeris
points that do not pass the I and Q channel tests are not used in the interpolation routine.
The average and the standard deviation from the nominal value for the semimajor axis,
inclination and angular momentum for the scene to be processed are saved for trending.
The equations to calculate the satellite's semimajor axis, inclination, and angular
momentum are as follows:
Angular momentum =  R x V 
where R is the satellite's position vector and V is the satellite's velocity vector
Inclination = acos(H dot k / H)
where H = R x V and k is the unit vector in the zdirection.
Semimajor axis = µ / (2.0 * E)
Version 3.2
29
7/9/98
IAS Geometric ATBD Version 3
where E = V2 / 2.0  µ / R and µ = G * M, where G is the Earth's gravity constant and
M is the mass of the Earth.
Validate Spacecraft Attitude Data
The onboard computer (OBC) calculates a flight segment attitude estimate every 512
milliseconds. The OBC outputs one of eight sets of spacecraft attitude data in the
telemetry every 4.096 seconds. The attitude information is output in the form of Euler
parameters that specify the vehicle attitude relative to the Earthcentered inertial frame
(J2000). The spacecraft's attitude is calculated using star crossings, gyro data and ADS
data. Care must be taken when using the spacecraft's attitude data, because star crossing
updates may cause discontinuities. The IAS converts the Euler parameters to roll, pitch,
and yaw to be used in the model.
The I and Q channel values are compared and any differences flagged. The spacecraft's
attitude data should not deviate greatly from a calculated value which was found using a
linear interpolation of the values around the value to be checked. The linear interpolated
value will be used to validate the spacecraft's attitude. Any large deviations are flagged.
If a large deviation is detected, the Q channel is validated. If the Q channel value passes
the validation test, it will be used instead of the I channel's value. If the attitude data in
both the I and Q channels have large deviations, the interpolated value is used. The Euler
parameters should satisfy the following equation.
EPA12 + EPA 22 + EPA 32 + EPA 42 = 1 +/ ε
where EPA1, EPA2, EPA3, and EPA4 are the Euler parameters which have been converted
by the LPS. The ε is some very small error term due to system noise.
The I channel values are validated using the above equation and any deviation is flagged.
If a deviation occurs, the Q channel values are checked. If both channels fail the validation
test, the interpolated value for each of the Euler parameters will be used.
The averages and standard deviations of the Euler parameters and the deviation from the
equation above are stored for trending.
Validate Gyro Data
The I and Q channel values are compared and any differences flagged. The gyro data
should not deviate greatly from a calculated value. This calculated value is found using a
forward and backward differencing of the values around the value in question. The steps
involved are:
Version 3.2
30
7/9/98
IAS Geometric ATBD Version 3
1) Check the current gyro value against a tolerance.
2) Do a forward prediction for each point using the two previous points. Compare this
predicted value against observed value.
3) Do a backward prediction for each point using the next two points. Compare this
predicted value against observed value.
Any point that passes (1) and either (2) or (3) is a good point. If a deviation is detected
in the I channel, the Q channel is checked. If the Q channel passes the validation test, it is
used in the processing. If both the I and Q channel values fail, the value in question is
flagged as an outlier. Once all values are checked the flagged values are replaced. Linear
interpolation is used to replace all outliers. The first valid points before and after the
outlier are found. These two values are then used to interpolate all outliers that lie in
between. Note: Care must be taken to account for counter resets in this check. The
register resets after (223  1) positive and 223 negative.
Validate Gyro Drift Data
The I and Q channel values are compared and any deviations are flagged. The values for
each major frame are compared. Any changes in the I channel are compared to the values
in the Q channel to determine if the changes have occurred in the Q channel. If both
channels display the same change and the change is determined to be within an acceptable
level, then it is assumed that a star sighting has occurred. The PCD major frame that
displays the change is flagged and the magnitude of the change is saved for trending.
Validate ADS Data
The I and Q channel values are compared and any differences flagged. The ADS data
should not deviate greatly from a calculated value which was found using a linear
interpolation of the values around the value in question. Any large deviations from the
calculated value are flagged. If a large deviation is detected, the Q channel value is checked
using the same interpolation method. If the Q channel's value passes the validation test,
the Q channel's value replaces the I channel's value. If both channels fail the test, the
value in question is corrected using the interpolated value.
The Linear interpolation test has been verified using a empirical test of Landsat 5 data.
The results showed that the difference between linear predicted and actual ADS data
points are typically within 10 counts. For a test scene over Iowa, once 5 actual outliers
were removed which had deviations of several thousand counts, the maximum deviations
in roll, pitch, and yaw were 12, 11, and 11 counts, respectively.
Validate Spacecraft Time Correction
Version 3.2
31
7/9/98
IAS Geometric ATBD Version 3
The I and Q channels values are compared and any differences flagged; the number of
occurrences will be saved for trending. The Mission Operations center will attempt to
maintain the onboard spacecraft clock to 145 milliseconds. A validation test will
calculate the correction of the clock and this value should not exceed the 145 millisecond
value. If the absolute value of the clock correction is much greater than 145 milliseconds,
the Q channel values are checked. If the Q channel clock correction values are less than
145 milliseconds, the Q channel’s values are used. If both the I and Q channel values are
greater than 145 milliseconds, then the clock coefficients are set to zero.
The clock correction is:
dt = spacecraft time  clock update time
ClockCorrectionI = TimeCoef1I + TimeCoef2I * dt + 0.5 * TimeCoef3I * dt 2
Validate and Save Instrument Time On
The I and Q channels values for the instrument ontime are compared. Any differences are
flagged and the values saved. The difference between the instrument on time and the start
of the first major frame is saved for trending and analysis.
Correct PCD Spacecraft Time
Add the clock correction to the major frame times.
Other Geometric PCD Parameter Validations
Other PCD geometric processing parameters, such as the SME mode and gyro unit flags
are validated by comparing the I and Q channel values. Any differences are flagged and
reported. The I channel is used as the default value.
3.1.5.3.2 MSCD Processing
For Landsat 7, the MSCD data is contained in the HDF formatted Level0R product.
The counted line length, scan direction, First Half Scan Error (FHSERR) and Second Half
Scan Error (SHSERR) are associated with the previous scan. The scan start time is,
however, for the current scan. Due to the counted line length, scan direction, FHSERR
and SHSERR being associated with the previous scan, the number of MSCD records must
be one more than the number of scans. The MSCD for the output product is subsetted to
match the imagery when the L0R product is generated.
The values from the MSCD required for ETM+ processing are:
FHSERR (First Half Scan Error)
SHSERR (Second Half Scan Error)
Scan Direction
Version 3.2
32
7/9/98
IAS Geometric ATBD Version 3
Scan Start Time
Counted line length
Validate Scan Direction
The I and Q channels are checked for consistency and any differences are flagged and
reported in the processing history file. The IAS assumes the I channel has the correct
value.
The flag for scan direction is a 1 or a 0. The scan direction is for the previous scan. The
validation test checks the first scan direction flag and the direction flags there after. Any
errors in the direction flags are corrected using the first valid direction flag as a reference.
The errors are flagged and reported in the processing history file and are also saved for
trending analysis.
Validate FHSERR and SHSERR
The MSCD is transmitted in both the I and Q channels of the telemetry data. Under
normal conditions the I and Q channels should be identical. Therefore, a comparison of
the FHSERR and SHSERR values contained in the I and Q channels will be used as one of
the validation tests. Any differences detected are flagged and reported in the processing
history file. The I channel is assumed to be correct unless the other validation check
determines it to be in error.
The FHSERR and SHSERR values should not deviate greatly from their nominal values.
Once these nominal values for the FHSERR and SHSERR have been characterized, the
FHSERR and SHSERR are checked for deviations from their nominal value. Large
deviations are flagged and the average difference and its standard deviation are saved for
trending. If a large deviation is detected, the Q channel is validated. If the Q channel
passes the validation check, the Q channel value is used in processing.
Validate Scan Start Time
The I and Q channels are checked for consistency and any differences are flagged and
reported in the processing history file. The IAS assumes the I channel has the correct
value.
The difference in time between start of scan times should not vary greatly from the
nominal value. The start of scan time is validated by comparing the difference in time
between scans. Any large variation from nominal is flagged and the Q channel is checked
for its difference between start times. If the Q channel values pass the validation check, it
will be used during processing. Otherwise, the scan start time is corrected using the
Version 3.2
33
7/9/98
IAS Geometric ATBD Version 3
nominal value. The average difference and the standard deviation between the scans times
for the scene being processed are saved for trending.
Validate Counted Line Length
The I and Q channel values are compared and any differences are flagged and reported in
the processing history file. The IAS will assume the I channel has the correct value.
The counted line length is the number of minor frames sampled from the start of scan to
the end of scan mark. The counted line length should be the same as the truncated
calculated line length which is found using the FHSERR, SHSERR, DWELLTIME (and
the conversion of counts to seconds). The counted line length should not deviate greatly
from its nominal value. The validation algorithm for the counted line length checks the
deviation of the counted line length from its nominal value and its deviation from the
truncated calculated line length. Large deviations from the nominal counted line length and
any deviation from the calculated line length are flagged. If a large deviation is detected,
the Q channel counted line length is compared to the calculated line length and if it passes,
the Q channel is used in the processing. The average counted line length and standard
deviation of the counted line length are saved for trending.
Calculate the first half scan time, the second half scan time, and the total scan time using:
For each forward scan:
first_half_time
= T fhf  FHSERR * Tunit
second_half_time = T shf  SHSERR * Tunit
total_scan_time = first_half_time + second_half_time
For each reverse scan:
first_half_time = T fhr  FHSERR * Tunit
second_half_time = Tshr  SHSERR * Tunit
total_scan_time = first_half_time + second_half_time
Calculate the line length using:
calculated_line_length = total_scan_time / DWELLTIME
where:
T fhf
T shf
T fhr
T shr
T unit
The forward nominal first half scan time
The forward nominal second half scan time
The reverse nominal first half scan time
The reverse nominal second half scan time
The conversion factor from counts to time
Correct Spacecraft Time in the MSCD
Version 3.2
34
7/9/98
IAS Geometric ATBD Version 3
The correction value for the spacecraft time is found in the process PCD module. This
value must be added to the MSCD start of scan times.
3.1.5.3.3 Platform/Sensor Model Creation
The Platform/Sensor model creation algorithm reads the telemetry data from the PCD and
the MSCD. These data are converted to a form which is usable by the geometric model.
The gyro and ADS data are combined and converted to the common Navigational
Reference Base system to create a time ordered table of attitude information. The
ephemeris data is converted to the ECEF coordinate system and fit with a polynomial
which is a function of time. The mirror scan 5th order polynomial is corrected for each
scan's variation in velocity.
Calculate Satellite State Vector
Transformation from J2000.0 system to the Earth Fixed system: For each valid
ephemeris data point, the following coordinate transformation is used to convert the
satellite position and velocity from the inertial system (J2000.0 system for Landsat 7; the
trueofdate system was used for Landsat 4, 5) to the Earth Fixed system:
Transform the J2000.0 system to the meanofdate system through the precession
matrix (Landsat 7 only).
Transform from the meanofdate system to the trueofdate system through
nutation matrix (Landsat 7 only).
Transform from the trueofdate system to the Earth Fixed system through Earth
rotation and polarmotion matrix. The Earth Orientation Parameters (UT1UTC,
xp, yp) are passed in from the Calibration Parameter File.
Generation of polynomial for orbit interpolation:
Two sets of polynomial coefficients are generated, one for the orbit position and
velocity in the Earth Fixed system and the other for the orbit position and
velocity in the Inertial system (J2000.0 for Landsat 7; trueofdate system for
Landsat 5). The methods for generating the coefficients for the two sets are the
same: solving the Vandermonde equation system using an algorithm given in
Reference 6. The order of the polynomial is automatically chosen as the number
of valid ephemeris data points in the scene, i.e., an evendetermined fit with no
redundancy.
Version 3.2
35
7/9/98
IAS Geometric ATBD Version 3
The polynomial in the Earth Fixed system is used to interpolate the orbit position and
velocity at each (grid point) sensor time, for the purpose of constructing a look vector and
intersecting the Earth model in the Earth Fixed system.
The polynomial in the Inertial system is used to interpolate the orbit position and
velocity at each gyro data point time, for the purpose of correcting the portion of gyro
observation caused by satellite motion.
Process Gyro Data
The Inertial Measurement Unit (IMU) which houses the gyroscopic units are sampled
every 64 msec for each of the three axes. The sampled values are referred to as the gyro
data. This results in 64 samples per axis per PCD major frame or 256 samples per axis
per PCD cycle. The gyro data registers reset if a positive count occurs when the register's
value is a positive 223  1, or if a negative count occurs when the register's value is a
negative 223. The time of the gyro data samples are referenced from the start of a PCD
cycle by the equation.
Gyro_data_sample_time = 64N milliseconds,
where N = 0 ... 255
The OBC calculates the gyro drift rates for each of the axes, using the results of star
sightings. The results of the OBC gyro drift calculations are in the ACS system. The
gyro drift parameters are updated asynchronously at a rate of up to once per minute.
When a gyro drift value changes, the value is updated at the PCD time code for the cycle
minus 8.192 seconds.
Processing Summary
The highlevel flow of gyro processing is as follows:
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
The register resets are accounted for.
For each axis, the first gyro value is subtracted from the rest.
The gyro data are rotated to the ACS system.
The orbital drift component are removed.
The gyro, gyro drift, and quaternion data are combined.
The data from step 5 are rotated to the ADS system.
The data from step 6 are filtered using the Gfilter.
The data from step 7 are synchronized with the ADS data.
Synchronize the ADS data.
The data from step 8 are combined with the ADS data.
The combined data are filtered using the Ffilter.
The combined data are rotated to the ACS system.
Version 3.2
36
7/9/98
IAS Geometric ATBD Version 3
For more detail on step 5 set section 3.1.5.3.3.1 on Smoothing Low Frequency
Information.
For more detail on steps 6, 8, 9, 10, 11 and 12, see section 3.1.5.3.3.2 on Combining Low
and High Frequency Attitude Information.
Processing Gyro Data to a Common (ADS) Reference System
Processing of the gyro data to a common ADS reference system is as follows:
1. Since the LPS converts the gyro data counts to arcseconds, the register reset
correction must be accomplished using arcseconds as the units. Thus, a positive reset
occurs at 511705.027 arcseconds and a negative reset occurs at 511705.088 arcseconds. Once a reset is detected, the values following the reset are offset by the
appropriate reset value.
2. Due to the gyro data measuring changes in attitude, the first value for each of the axes,
which corresponds with the spacecraft’s attitude values, is subtracted from the rest of
the values.
3. Each axis of the gyro data is rotated to the ACS reference system using a prelaunch
derived transformation matrix. This is accomplished using the following:
Define A as the rotation matrix due to roll. Define B as the rotation matrix due to
pitch and define C as the rotation matrix due to a yaw. Then the Perturbation
matrix P is defined as:
P = CBA
Use the Transformation matrix from the gyro axes to the ACS [IMU2ACS] and
its inverse to find the P matrix in the ACS reference system or P':
P' = [IMU2ACS] P [IMU2ACS]1
a b c
= d e f
g h i
Then the transformed roll (r') is:
r' = tan 1( h / i)
where h and i are elements of P'
Then the transformed pitch (p') is:
Version 3.2
37
7/9/98
IAS Geometric ATBD Version 3
p' = sin1(g)
where g is an element of P'
and the transformed yaw (y') is :
y' = tan1( d / a) where d and a are elements of P'
4. The drift of the gyro data due to the orbital motion is removed using the polynomial
for the ephemeris data in the inertial system derived in the Calculate Satellite State
Vector subroutine.
5. Combine the gyro, gyro drift and quaternion data into one data stream.
(Details of step 5 can be found in section 3.1.5.3.3.1, Smoothing Low frequency
Attitude Information)
6. Each axis of the gyro data is transformed to the ADS reference system. This is
accomplished by combining the ACS to ETM+ transformation matrix and the ETM+
to ADS transformation matrix. The ACS to ETM+ transformation matrix is the
inverse of the ETM+ to ACS matrix. The ETM+ to ACS matrix is first determined
prelaunch and is updated after launch. Its values will be contained in the Calibration
Parameter file. The ETM+ to ADS matrix is determined prelaunch and its members
are contained in the Calibration Parameter File.
The transformation from the ACS system to the ADS is accomplished using the
following:
Define A as the rotation matrix due to roll. Define B as the rotation matrix due to
pitch and define C as the rotation matrix due to a yaw. Then the Perturbation
matrix P is defined as:
P = CBA
First combine the Transformation matrix from the ACS to the ETM+ T and the
transformation matrix from the ETM+ to the ADS system R. Calculate the
combined matrix S inverse. Then using S and its inverse to find the P matrix in
the ADS reference system or P':
S = RT
Version 3.2
38
7/9/98
IAS Geometric ATBD Version 3
P' = SPS 1 =
a b c
d e f
g h i
Then the transformed roll (r') is:
r' = tan 1(h / i)
where h and i are elements of P’
Then the transformed pitch (p') is:
p' = sin1(g )
where g is an element of P’
Then the transformed yaw (y') is:
y' = tan1(d / a)
where d and a are elements of P’
7. Apply the Gfilter. (Refer to Section 3.1.5.3.3.1).
Steps 7,8, 9, 10, and 11 are shown in more detail in the section on attitude processing,
Section 3.1.5.3.3.2.
12. Calculate the inverse of step 7, P’’ = S1P’S
Spacecraft Attitude Processing
The OBC calculates a flight segment attitude estimate every 512 milliseconds. The OBC
outputs one of the eight sets of data in the telemetry every 4.096 seconds or once a PCD
major frame. The attitude is output in the form of Euler parameters that specify the
vehicle's attitude relative to the EarthCentered Inertial reference frame. The Euler
parameters (EPA 1, EPA2, EPA3, EPA4) are components of the reference quaternion, as
propagated from the gyro data, which defines the spacecraft's attitude. Components 1
through 3 define the Eigen axis of the rotation in ECI coordinates, and the fourth
component defines the rotation about the Eigen axis.
There are four attitude estimates in each PCD cycle. The time associated with the
attitude data contained within the PCD can be derived from the time code contained in
the words 96 through 102 of the first PCD major frame in the cycle. The derivation for
Landsat 7 (see Reference 1, page 3−29) is as follows:
Version 3.2
39
7/9/98
IAS Geometric ATBD Version 3
PCD MajorFrame
Time Computation
PCD time code − 8.192 seconds
PCD time code − 4.096 seconds
PCD time code
PCD time code + 4.096 seconds
0
1
2
3
Processing
The PCD attitude quaternions define the attitude of the spacecraft at a reference time.
This attitude measurement reflects the gyro information up to that point. The IAS
determines the spacecraft low frequency attitude by processing the gyro, gyro drift and
quaternions into one data stream by using a Kalman filter/smoother. This integrates
attitude measurements from all time into a single best estimate of the attitude over an
entire interval. The average associated with this data stream serves as the scene reference
attitude state. The average is removed from the data stream. This data stream is then
combined with the ADS information to get the change in state of the attitude from this
reference state. To accomplish this, the Euler parameters need to be converted to ACS
roll, pitch, and yaw values.
The direction cosines (transformation) matrix from the ACS reference axis to the ECI
reference system (ACS2ECI) is defined as:
ACS2ECI =
EPA 2  EPA 2  EPA 2 + EPA 2 2(EPA EPA  EPA EPA ) 2(EPA EPA + EPA EPA )
1
2
3
4
1
2
3
4
1
3
2
4
2
2
2
2
2(EPA EPA + EPA EPA )  EPA + EPA  EPA + EPA
2(EPA EPA  EPA EPA )
1
2
3
4
1
2
3
4
2
3
1
4
2
2
2
2
2(EPA1EPA 3  EPA 2 EPA 4 ) 2(EPA 2 EPA 3 + EPA1EPA 4 )  EPA1  EPA 2 + EPA 3 + EPA 4
The ACS2ECI transformation matrix can also be defined as the product of the inverse of
the spacecraft's attitude perturbation matrix P and the transformation matrix from
the orbital reference system to the ECI reference system (ORB2ECI).
ACS2ECI = [ORB2ECI][P 1]
The orbital reference system to ECI matrix must be defined at the same time as the
spacecraft's attitude matrix. The ORB2ECI matrix is constructed of the negative
of the spacecraft's position vector divided by its magnitude (zdirection vector),
the cross of the zdirection vector with the spacecraft's velocity vector
divided by its magnitude (ydirection vector) and the cross of the ydirection
vector with the zdirection vector divided by its magnitude (xdirection vector).
Version 3.2
40
7/9/98
IAS Geometric ATBD Version 3
The resulting ORB2ECI matrix is:
X x Yx Z x
ORB2ECI = X y Yy Z y
X z Yz Z z
Where Zx is the x component of the zdirection vector,
Z y is the y component of the zdirection vector,
Z z is the z component of the zdirection vector,
Xx is the x component of the xdirection vector,
Xy is the y component of the xdirection vector,
Xz is the z component of the xdirection vector,
Yx is the x component of the ydirection vector,
Yy is the y component of the ydirection vector,
Yz is the z component of the ydirection vector.
The roll, pitch and yaw values are contained in the P1 matrix, thus:
P1
a b c
= [ORB2ECI] 1[ACS2ECI] = d e f
g h i
Then the spacecraft roll (r) is:
r = tan 1(h / i)
where h and i are elements of P1
Then the spacecraft pitch (p) is:
p = sin1 (g)
where g is an element of P1
Then the spacecraft yaw (y) is:
y = tan1(d / a)
where d and a are elements of P1
Convert ADS to Common Reference System
The Attitude Displacement Sensor Assembly (ADSA) consists of three nominally
orthogonal ADS sensors. The ADSA is mounted on the ETM+ telescope. A digital
sample count of 0 is the maximum positive angular displacement and a digital count of
4095 is the maximum negative angular displacement. The LSB of each count is
Version 3.2
41
7/9/98
IAS Geometric ATBD Version 3
250 / 211 microradians. The nominal zero angular displacement output of the ADSA is
2048 counts. Each axis of the ADSA is sampled every 2 milliseconds, resulting 8192
samples of each axis in a PCD cycle. The axes are sampled sequentially as given below:
ADSA AXIS
roll
pitch
yaw
Sampling Times
PCD Time Code + (2N + 3/8) milliseconds
PCD Time Code + (2N + 7/8) milliseconds
PCD Time Code + (2N + 11/8) milliseconds
where N = 0,1,2,...,8191 for the given axis
Each axis of the ADSA has a nominal 2.0 to 125.0 Hz bandwidth. The transfer function
to rotational motion is measured prelaunch and is supplied in the Calibration Parameter
File. The nominal relative alignment matrix between the ADS and the ETM+ Sensor is
determined prelaunch and is also supplied in the Calibration Parameter File.
Processing
The first three steps are summarized here. Refer to Section 3.1.5.3.3.2, Combining Low
and High Frequency Attitude Information, for more information.
1. Each of the samples for each axis is transformed to rotational motion using the
transfer function.
2. Synchronize the ADS data
3. Combine the ADS data with the smoothed low frequency attitude data.
4. Rotate to the ACS system. Due to the ADSA being mounted on the ETM+
telescope, the transformation to the ACS reference system requires the ADSA
measurements to be transformed to the ETM+ reference system and then transformed
through the ETM+ sensor to ACS alignment matrix. Since the ETM+ to the ACS
reference system matrix will be determined after launch, this transformation must be a
two step process.
Define the transformation matrix from the ADSA to the ETM+ sensor as
ADS2ETM and the transformation matrix from the ETM+ to the ACS as
ETM2ACS. Then the transformation matrix between the ADSA and the ACS
(ADS2ACS) is:
ADS2ACS = [ETM2ACS][ADS2ETM]
Version 3.2
42
7/9/98
IAS Geometric ATBD Version 3
To transform the axis of the ADSA, the perturbation matrix P must first be
calculated. To accomplish this, we must define a matrix A as the rotation matrix due
to roll, a matrix B as the rotation matrix due to pitch and a matrix C as the rotation
matrix due to yaw. Thus, the perturbation matrix is:
P = CBA
Then using the transformation matrix ADS2ACS ands its inverse, the perturbation
matrix in the ACS reference P' is found:
a b c
P' = [ADS2ACS]P[ADS2ACS]1 = d e f
g h i
Then the transformed roll (r') is:
r' = tan 1(h / i)
where h and i are elements of P'
Then the transformed pitch (p') is:
p' = sin1(g )
where g is an element of P'
Then the transformed yaw (y') is:
y' = tan1(d / a)
where d and a are elements of P'
Calculate Scan Mirror Correction
The scan mirror (as introduced in section 3.1.5.1.2) is modeled using a linear motion with
a 5th order polynomial correction. A 5th order polynomial is determined during the
ground testing of the instrument, this along with an adjustment due to high frequency roll
jitter components combine to give an overall 5th order polynomial correction. Along with
the high frequency jitter component, in orbit perturbations within the satellite causes the
scan mirror velocity to change from scan to scan. The scan mirror has sensors attached to
it to measure the mirror position at the start of scan, midscan, and end of scan. The data
from these sensors are output in the form of counts, which are the deviation from nominal
scan mirror position in time. For example, the scan mirror crossed the midscan detector
so many counts ahead (or behind) of the nominal time. The unit of the counts are
16/84.9037e6 seconds or 0.18849 microseconds.
In the IAS model of the scan mirror motion, the mirror starts at a positive angle (A sm) and
becomes more negative as a function of time, until it reaches A me, where A sm is the angle
from the start of scan to the midscan position and A me is the angle of the mirror from
Version 3.2
43
7/9/98
IAS Geometric ATBD Version 3
midscan to end of scan. This sign convention is used to place the scan mirror motion in
the spacecraft coordinate system. Note: This convention is different than the convention
used in Reference 1 and 2.
Remaining in the true time domain, for the forward and reverse scan direction, a positive
result from the polynomial corresponds to the mirror being ahead of its linear location and
a negative result corresponds to the mirror being behind its linear location. However,
whether the results are added to the linear portion or subtracted depends on the slope of
the linear motion. For example, using the IAS sign conventions, adding a positive
polynomial value to the forward scan direction causes the scan mirror to be behind its
linear location. Thus, for the forward scan direction, a positive polynomial value must be
subtracted from the linear location or the sign of each of the coefficients must be changed.
The latter was chosen. For a reverse scan, a positive polynomial value added to the linear
location corresponds to the mirror being ahead of it linear location. Thus, adding the
polynomial value is correct. Note: The MSCD values of scan direction, counted line
length, first half scan error and second half scan error are for the previous scan.
Define:
bj = 0, 5
dj = 0,5
Asmf
Amef
Asmr
Amer
T fhf
T shf
T fhr
T shr
Tn
T unit
mprofilei,j
The forward scan polynomial coefficients for the SME mode
The reverse scan polynomial coefficients for the SME mode
The forward start to middle angle for the SME mode
The forward middle to end angle for the SME mode
The reverse start to middle angle for the SME mode
The reverse middle to end angle for the SME mode
The forward nominal first half scan time
The forward nominal second half scan time
The reverse nominal first half scan time
The reverse nominal second half scan time
The nominal scan time
The conversion factor from counts to time
The corrected mirror polynomial for each scan
Compute the actual scan timefind true first half time (T fh), true second half time (Tsh),
and total scan time (ts):
For a forward scan:
T fh = T fhf  FHSERRi * T unit
T sh = T shf  SHSERRi * T unit
t s = T fh + T sh
(The first half scan time)
(The second half scan time)
( The total scan time)
For a reverse scan:
T fh = T fhr  FHSERRi * T unit
T sh = T shr  SHSERRi * T unit
(The first half scan time)
(The second half scan time)
Version 3.2
44
7/9/98
IAS Geometric ATBD Version 3
t s = T fh + T sh
( The total scan time)
Calculate the correction to the linear model needed to make the midscan angle equal to the
observed value (zero) at time Tfh. (see Figure 3.1.55)
For a forward scan:
Af = −(Amef T fhf
+
AsmfT shf) / ts
For a reverse scan:
Ar = −(Amer T fhr
+
AsmrT shr) / ts
Figure 3.1.5 5 Scan Mirror Profile Deviations
Rescale the fifthorder polynomial coefficients to account for the difference between the
actual and nominal scan time:
For a forward scan:
bj = bj * (T n / ts) j where j = 0 to 5
For a reverse scan:
dj = dj * (T n / ts)j where j = 0 to 5
Version 3.2
45
7/9/98
IAS Geometric ATBD Version 3
Calculate the amplitude component of the correction quadratic (B f) due to the
polynomial location of the mirror to the linear location. This is done by evaluating the
rescaled fifthorder polynomial at time Tfh. ( f(b,Tfh) in figure 3.1.55 ).
For a forward scan:
Bf = b0 + (b1 * T fhf ) + (b2 * T fhf 2) + (b3 * T fhf 3) + (b4 * T fhf 4) + (b5 * T fhf 5)
For a reverse scan:
(Note: a () sign is applied to the d coefficients)
Br = d0 + (d1 * T fhr ) + (d2 * T fhr 2) + (d3 * T fhr 3) + (d4 * T fhr 4) + (d5 * T fhr 5)
The net correction to the mirror profile at the midscan time (T fh) is (Df and Dr) is:
For a forward scan:
Df = A f  Bf
For a reverse scan:
Dr = A r  Br
Using the constraints a(0) = a(ts) = 0.0,
the coefficients of the quadratic are:
and
a(tsh) = D,
For a forward scan:
a2 = D f / (Tfhf * T shf)
For a reverse scan:
a2 = D r / (Tshr * T fhr )
For both scan directions:
a1 = t s * a2
a0 = 0.0
Integration of Roll Jitter into the Scan Mirror Model
The high frequency roll attitude components (jitter) must be considered together with the
scan mirror profile since the scan mirror is not rigidly attached to the instrument structure
in the along scan direction but, rather, is driven by a flex pivot. This acts as a mechanical
buffer which does not transmit the high frequency along scan jitter directly to the mirror.
The roll jitter terms therefore have two effects: the ETM+ instrument rotates relative to
the orbital coordinate system and the scan mirror rotates relative to the ETM+ instrument
(in the opposite direction). Effectively, the scan mirror does not “see” the high frequency
roll jitter. The roll rotation of the ETM+ instrument (i.e., the focal planes) independent of
the scan mirror causes the effective line of sight to be deflected in the opposite direction
from that caused by a rigid rotation of the entire focal plane/scan mirror assembly.
[Reference 35] This is depicted in Figure 3.1.56.
Version 3.2
46
7/9/98
IAS Geometric ATBD Version 3
The scan mirror angle relative to the rest of the ETM+, and to the orbital coordinate
system, is known at the three points where it is measured by the scan angle monitor
(SAM): scan start, midscan, and scan end. The line of sight at these three times will be
based on the known mirror angles at beginning, middle, and end of scan modified by the
high frequency roll at the corresponding times, so the combined mirror/roll jitter model
must match these known angles at these three times. This is accomplished by modifying
the scan mirror profile for each scan to include the scan start, midscan, and scan end high
frequency roll effects. This is implemented as a roll jitter correction to the scan mirror
polynomial: θ$(t ) . The differential jitter relative to this profile is then: θt − 2θ$(t )
where θt is the measured roll jitter and θ$(t ) is doubled to account for the fact that a scan
mirror rotation of θ$(t ) leads to a line of sight deflection of 2θ$(t ) . As described above
and depicted in Figure 3.1.56, this differential jitter leads to a deflection of the line of
sight in the opposite direction, or: − [θt − 2θ$(t )] . So, the net deflection of the line of sight
relative to the orbital coordinate system is: − θt + 2θ$(t ) . We implement this deflection by
inverting the sign on the roll jitter component (θt ) and including the scan mirror profile
adjustment (θ$(t )) in the mirror model. All that remains is to determine what the scan
mirror profile adjustment should be.
As noted above, the net line of sight deflection must account for the roll angles at scan
start (θs ) , midscan (θm ) , and scan end (θe ) . If the mirror profile adjustment is
constructed so that it is equal to the measured roll jitter at scan start, midscan, and scan
end (i.e., it modifies the mirror profile to include roll jitter effects at those three points)
then the net line of sight deflection at T=0, T=tfh, and T=Ts will be:
θnet = −θt + 2θt = θt where t = s, m, or e for scan start, midscan, and scan end.
This is the desired result, so the mirror correction term, θ$(t ) , must be constructed so
that:
θ$(0) = θs
θ$(t fh ) = θm
θ$(Ts ) = θe
Version 3.2
47
7/9/98
IAS Geometric ATBD Version 3
Scan Mirror
Focal Plane
Nominal LOS
a) Nominal Focal Plane/Scan Mirror Geometry
Scan Mirror
Focal Plane
Rotated
Scan Mirror
θ
Rotated
Focal Plane
θ Rotated LOS
Nominal LOS
b) Rigid Rotation of Focal Plane and Scan Mirror
θ
Scan Mirror
Rotated
Scan Mirror
Rotated
Focal Plane
Deflected LOS θ θ Rotated LOS
Nominal LOS
c) Rotation of Focal Plane Independent of Scan Mirror
Figure 3.1.56 Effect of Roll Jitter on Line of Sight Direction
Version 3.2
48
7/9/98
IAS Geometric ATBD Version 3
This can be accomplished using a correction quadratic similar to that used to perform the
midscan correction. The coefficients of this polynomial can be shown to be:
c0 = θ s
θe − θs (θe − θm )t fh + (θs − θm )t sh
c1 =
−
Ts
t fh t sh
c2 =
(θe − θm )t fh + (θs − θm )t sh
Ts t fh t sh
This makes the net mirror profile:
b0′ = b0 + a 0 + c0
b1′ = b1 + a1 + c1
b2′ = b2 + a 2 + c2
b3′ = b3
b4′ = b4
b5′ = b5
where the bi are the nominal mirror profile coefficients normalized for actual scan time,
the ai are the midscan correction quadratic coefficients computed from the first half and
second half scan errors, and the ci are the roll jitter correction quadratic coefficients
defined above. Inverting the sign on the measured roll jitter data and using this modified
mirror profile creates the desired net line of sight deflection.
3.1.5.3.3.1 Smoothing Low Frequency Attitude Information
Purpose:
In the IAS release 1 model the quaternion and gyro drift observations from the
PCD major frame immediately preceding the scene was used as a reference point. All
subsequent attitude data for the scene was referenced from this point, therefore if a stellar
observation occurred later in the scene this change in reference frame was not be reflected
in the model. This change can be thought of as an update in the reference frame of the
system. The update could produce a discontinuity in the attitude data stream. To
incorporate all of the stellar observation information in the PCD quaternion samples, a
Kalman filter is used as a post processing filter of the attitude data stream. The Kalman
filter is implemented in a forwardpredictor, backwardsmoother configuration. Under
this scenario, the filtering serves as a forward recursive sweep while the smoothing serves
as a backward sweep. The backwardsmoother processing propagates the stellar
observations, represented by the quaternion data, backward in time to provide a best
estimate of the spacecraft attitude given data from all times in the subinterval. The
smoothed data stream will be a new attitude data stream representing the combined gyro,
quaternion, and drift data. This smoothed data set will represent the low frequency
Version 3.2
49
7/9/98
IAS Geometric ATBD Version 3
attitude information of the satellite. The data is then combined with the ADS information
to produce the overall attitude of the sensor/satellite.
The rest of this section contains three parts. The first is a brief description of the
filtering and smoothing equations. The reader is encouraged to pursue the references given
later in the text for further details on the subject. The second portion will involve
calculating the equations that are specific to this application. The final portion will
explain the details of running the process with actual data.
Algorithm:
The discrete Kalman filter assumes a random process modeled as follows:
[ X ]K +1 = [φ ]K [ X ]K + [q ]K
Where
[X]K = (n x 1) state vector at time t K
[φ]K = (n x n) state transition matrix relating XK to X K+1
[q]K = (n x 1) process noise, white, uncorrelated, with covariance Q
Measurements of the process are modeled as:
[Z ]K = [H ]K [ X ]K + [v ]K
Where
[Z]K = (m x 1) measurement vector at time tk
[H]K = (m x n) observation matrix relating state vector at tk to measurement vector
[v]K = (m x 1) measurement error, white, uncorrelated, with covariance R
The first step involves using the Kalman filter as a forward prediction process. The
Kalman filtering equations are (note all equations refer to the discrete domain):
(Note: superscript P refers to prediction)
Given an initial state vector and covariance matrix [X]P0 and [P]P0:
For k=1,…,N (number of data points)
Filter data
[K ]K
(
= [P ]K [H ]K [H ]K [P ]K [H ]K + [R ]K
P
T
P
T
)
−1
[ X ]K = [ X ]Kp + [K ]K ([Z ]K − [H ]K [ X ]Kp )
Version 3.2
50
7/9/98
IAS Geometric ATBD Version 3
[ P]K = ([I ] − [K ]K [H ]K )[ P]Kp
Predict
[ X ] Kp +1 = [φ ]K [X ]K
[P ]Kp +1 = [φ ]K [P ]K [φ ]TK
+ [Q ] K
Where:
[I] is the identity matrix
[K] is the Kalman Gain matrix
[P] is the state error covariance matrix.
[Q] = E[qq]
[R] = E[vv]
The second step involves using the filtered and predicted data to produce a smoothed data
set. Let,
[X]PK = estimate of [X] given measurements through tK1
[P]PK = error covariance associated with estimate [X]K
[X]K = filtered estimate through tK
[P]K = state error covariance associated with filtered estimate at tK
Using the definitions above a new notation is:
[X]KK1 = [X]PK
[P]KK1 = [P]PK
[X]KK = [X]K
[P]KK = [P]K
Smooth data
For K=N1,N2,..,0
[ X ]K  N = [ X ]K K + [ A]K ([ X ]K +1 N − [ X ]K +1K )
[ A]K = [ P]K K [φ ]K +1, K [ P]−K1+1K
T
Version 3.2
51
7/9/98
IAS Geometric ATBD Version 3
[ P]K  N = [ P]K K + [ A]K ([ P]K +1 N − [ P]K +1K )[ A]TK
For more information on the Kalman filter and smoothing data refer to references 23 and
34.
This is all the background on the Kalman filtering that will be given here. If more
information is desired the reader should refer to one of the references stated above. The
rest of this document will focus on the how the matrix equations specific to the attitude
processing were obtained and some of the processing flow issues.
The system will first be described for the continuous time case, then the corresponding
discrete time case will be solved for. The equations for the continuous case are:
•
X = FX + w
Y = CX
Where
X = state vector
F = state matrix
Y = output response vector
C = matrix relating state to measurement
w = vector of white noise processes
F,G,C may be time varying matrices.
The state matrix for the attitude processing contains three states. The states are attitude,
gyro rate, and gyro drift. The state transition matrix is:
0 1 0 X 1
[X ] = 0 0 0 X 2 +
0 0 0 X 3
•
0
w
2
w3
X1 = attitude
X2 = gyro rate
X3 = gyro drift
w = process noise (white)
The model above states that the gyro rate is the time derivative of the attitude and the
time derivatives of the gyro rate and drift are independent white noise processes.
The output response vector has the following form:
Version 3.2
52
7/9/98
IAS Geometric ATBD Version 3
1 0 0 X 1
[Y ] = 0 1 − 1 X 2 +
0 0 1 X 3
r1
r
2
r3
r = measurement noise (white)
The observation matrix above states that there are three kinds of measurements involved
in the process:
1) Quaternion observations which measure the attitude.
2) Gyro/gyro rate observations which measure the attitude rate biased by the gyro drift.
3) Gyro drift observations which measure the bias in the gyro measurements.
The gyro measurements have no absolute origin, therefore they have meaning only relative
to each other. By differencing the gyro samples and dividing by the sampling rate they
can be converted to attitude rate measurements, the quaternions are then used as the
reference. Once the gyro samples are converted into rate the drift can be directly applied
since it is measured in units of radians per unit of time. The measurements then used by
the process, as measured by the satellite, are quaternions, gyro rate, and gyro drift.
The matrix [F] does not contain any time varying elements and is small
enough that the direct method can be used to determine the discrete transition matrix.
The equation takes the form:
[φ ]k
[ {[[I ]s − [F ]] }]
−1
= L−1
t = dt
Where
L1 is the inverse Laplace transform.
[I] = identity matrix
[]1 = matrix inversion
dt = sampling interval
This gives the following
.
1
s
s − 1 0
= L−1 0 s 0 = L−1 0
0 0 s
0
−1
{
L−1 [[I ]s − [F ]]
−1
}
1
2
s
1
s
0
0
1 t 0
0 = 0 1 0
0
0
1
1
s
Substituting dt for t in the above matrix gives the discrete state transition matrix.
Version 3.2
53
7/9/98
IAS Geometric ATBD Version 3
The discrete process noise matrix can also be solved directly. The discrete process noise
is found from the equation:
[Q ]k = ∫0 ∫0 [φ (dt , v)][G (u )]E [w(u ) wT (v)][G (v)]T [φ (dt , v)]T dudv
dt
dt
Noting that w2 and w3 are uncorrelated and that σ is the variance of the associated white
noise process:
0 0
G (u ) E w(u ) w (v) G (v) = 0 σ 2
0 0
[
T
]
0
0
σ 3
Performing the remaining matrix operations:
1 u 0 0 0
0 1 0 0 σ
2
0 0 1 0 0
2
0 1 0 0 uvσ 2
0 v 1 0 = vσ 22
σ 3 0 0 1 0
uσ 2
σ 22
2
0
0
0
σ 32
Solving the double integral for the matrix above gives the discrete process noise matrix:
uvσ 22
2
∫0 ∫0 vσ 2
0
dt dt
uσ 22
σ 22
0
dt 4σ
0 34 2
dt σ 2
0=
2
σ 32 0
dt σ 2
2
3
2
dt 2σ 22
0
0
2 2
dt σ 3
0
The observation matrix relating the state vector to the output/measurement vector is the
same for the continuous and discrete cases.
[H] = [C]
Under this system model, each axis can be processed independently. The attitude of the
satellite will be contained in the first state of the state vector. The gyro and quaternion
data for Landsat 7 is not sampled at the same instant in time or at the same rate. This will
need to be accounted for prior to and during processing. The data will also have to be in a
Version 3.2
54
7/9/98
IAS Geometric ATBD Version 3
common coordinate system. The drift is measured in the attitude control system (ACS).
The quaternions are given as Euler parameters and can be converted to the ACS (ATBD
section 3.1.5.3.3 Satelite Attitude Processing). The gyro data is measured in the inertial
measurement unit (IMU) axis and can be converted to the ACS (ATBD section 3.1.5.3.3
Processing Gyro Data). The system therefore works in the ACS. The output will also be
in the ACS. The time offset is accounted for by looking at the sampling times associated
with each data source.
Gyro samples
PCD time code + 0.064N seconds N=0,1,2,…,255
Attitude
PCD
Major Frame
0
PCD time code – 8.192 seconds
1
PCD time code – 4.096 seconds
2
PCD time code + 0.0 seconds
3
PCD time code + 4.096 seconds
Drift
PCD
Major Frame
0
PCD time code – 8.192 seconds
At the first PCD cycle, PCD major frame 0, the gyro data will correspond to the
quaternion value of the first PCD cycle, major frame 2. This will be the initial/start
conditions for the filtering. The gyro drift value will be zero until the data stream reaches
a point with valid drift data. The initial covariance error matrix, [P], was set equal to the
process noise matrix. This value may need to be changed after some test runs with
different data sets. Due to the fact that the gyro, attitude and drift data have different
sampling rates, there will be three different observation matrices. One observation matrix
when only the gyro is present, a second when both the gyro and quaternions are present,
and a third when all three measurements are available. When only the gyro is available,
the observation matrix will be:
X1
[Y ] = [0 1 − 1] X 2 + [r ]
X 3
When both the gyro and quaternions are available the observation matrix will be
Version 3.2
55
7/9/98
IAS Geometric ATBD Version 3
Y1 1 0 0 X 1
=
+
Y2 0 1 − 1 X 2
r1
r2
When the gyro, quaternion, and gyro drift measurements are available the observation
matrix will be:
Y1 1 0 0 X 1
Y = 0 1 − 1 X +
2
2
Y3 0 0 1 X 3
r1
r
2
r3
Note: Although it is stated that “gyro measurements are available” it is actually the gyro
rate that the system is working with at this point.
The state transition matrix will be the same for all these cases. The final step is to
determine the values associated with the measurement noise matrices and the standard
deviations associated with the process noise matrix. The measurement error for the gyro
corresponds to the increment associated with one gyro count or 0.061 arcseconds. The
measurement error for the gyro rate is found by dividing 0.061 arcseconds by the
sampling rate of the gyro, 0.064 seconds. The measurement error used for the quaternion
corresponds to the error associated with the Celestial Sensor Assembly and is 4.2 arcseconds [Reference 35]. The measurement error for the drift and the standard deviations
for the process noise were determined by trial and error. The criteria for choosing values
for trial however fell under the following premise.
The drift values were considered much less reliable then either the gyro or
quaternion values.
The measurement error for the quaternion was chosen to be larger than the gyro
(rate) measurement error. This allowed the quaternions to act as a very low
frequency reference and the gyro (rate) to represent the higher frequency
reference. This will keep the scantoscan consistency needed for processing.
The attitude processing should not effect the scantoscan accuracy to any
measurable degree.
The measurements for the gyro (rate) were considered to carry more weight then
the model while the opposite was considered for the drift. The values were
chosen such that more consideration was given to the model for the drift while for
the gyro (rate) more consideration was given toward the measurements.
Through trial and error this led to a values of 50.0 arcseconds for the drift measurement,
a standard deviation of 0.1 for the process noise associated with the gyro rate, and a
Version 3.2
56
7/9/98
IAS Geometric ATBD Version 3
standard deviation for 0.00001 for the process noise associated with the drift. It should
be noted that these values were determined from data associated with the Landsat 4
satellite, the values may need to change slightly to reflect Landsat 7. These values were
used to process one hundred consecutive major frames of Landsat 4 attitude data.
Processing the quaternion, gyro and drift produced a change of less than one meter in the
difference between adjacent gyro samples. This provided proof that the effect on
between scan alignment was negligible.
Once all the data is processed, filtered and smoothed, the new attitude data can be used
by the G and F filters, along with the ADS information, to get the overall satellite
attitude. To be consistent with the old way of processing however, the output from the
filtering and smoothing should have the first value kept as a frame of reference and
subtracted from all other values.
3.1.5.3.3.2 Combining Low and High Frequency Attitude Information
Much of the premise for the design of the Attitude Data processing algorithm comes from
the document Landsat Thematic Mapper Attitude Data Processing by G. J. Sehn and S.
F. Miller [Reference 10]. References 38, 39, and 40 are also good sources on the subject.
The following symbols will be used throughout this section:
FFT
Hg
P
Ha
PHa
G(f)
F(f)
function
j
Fast Fourier Transform with a complex component of 1 on the exponent
associated with the forward transform (1 for the reverse).
Frequency response of the gyro sensor
Frequency response of the ADS prefilter
Frequency response of the ADS sensor
Frequency response of the ADS sensor combined with prefilter
Frequency response of G filter
Frequency response of F filter
Magnitude response of function
Indicates the imaginary part of a complex number
The gyro has the following transfer function for Landsat 4 and 5:
Hg = ak * (1.0 + K * s) / (1.0 + tau * s) * [(s / wn)2 + (2.0 * zeta * s) / wn + 1.0]
where ak, K, tau, wn, and zeta are measured constants and
s = 2j * PI * freq.
The ADS and prefilter combined have the following transfer function for
Landsat 4 and 5:
Version 3.2
57
7/9/98
IAS Geometric ATBD Version 3
Hads = (A5 * s5 + A4 * s 4 + A3 * s 3) /
(s6 + B5 * s 5 + B4 * s 4 + B3 * s 3 + B2 * s 2 + B1*s + B0)
where A5, A4, A3, B5, B4, B3, B2, B1, and B0 are measured constants and
s = 2j * PI * freq.
The Landsat 4 and 5 satellites have two sensors to measure satellite attitude
perturbations. One of the sensors measures the low end frequencies while the other
measures high end frequencies. These are called the gyro and attitude displacement sensor
(ADS) respectively. The gyro is sampled every 64 milliseconds and the ADS is sampled
every 2 milliseconds. The samples taken from these two systems cannot be simply added
together to get the overall satellite attitude perturbations due to a discontinuity between
the transition regions of the responses of the two systems and due to the sensor phase
shift properties. A plot of the frequency responses of the two systems is shown in
Figure 3.1.57.
Figure 3.1.5 7 Magnitude Response of Gyro and ADS
The magnitude of the sum Hg + PHa is shown in 8 Figure 3.1.58.
Version 3.2
58
7/9/98
IAS Geometric ATBD Version 3
Figure 3.1.5 8 Magnitude Response of Gyro Plus ADS
The magnitude of the sum Hg + PHa is shown in Figure 3.1.59.
Figure 3.1.5 9 Magnitude Response of Gyro + Magnitude Response of ADS
These figures help to explain the discontinuity between the system responses.
Version 3.2
59
7/9/98
IAS Geometric ATBD Version 3
To combine these two signals, along with the quaternions and gyro drift values, a network
is used that combines the traditional G.J. Sehn and S.F. Miller approach along with the
Kalman filtering that was explained in section 3.1.5.3.3.1
Sync
Compensate
for Gyro
Reset
Rotate
to ACS
Apply
F
filter
Subtract
Orbital
Motion
Kalman
Filter
Figure 90
Rotate
to
ADS
Apply
G
Filter
S
y
n
c
Attitude Processing Network
This algorithm combined with the necessary steps involved in implementing a workable
system is listed below. A diagram of the algorithm steps is given in Figure 3.1.510.
1. The register resets are accounted for.
2. For each axis, the first gyro value is subtracted from the rest.
3. The gyro data is rotated to the ACS.
Version 3.2
60
7/9/98
Rotate
to
ACS
IAS Geometric ATBD Version 3
4. The orbital drift component is removed.
(Details of steps 1 through 4 can be found in section 3.1.5.3.3, Process Gyro Data)
5. Smooth gyro, gyro drift and quaternions.
(Details of step 5 can be found in section 3.1.5.3.3.1)
6. Rotate the data from step 5 to the ADS system.
(Details of step 6 can be found in section 3.1.5.3.3, Process Gyro Data)
7. Filter the data from step 6. This is done with a filter traditionally called the presumming G filter or just G filter. This filter will help to eliminate any frequencies
above a predetermined cutoff while also providing an adequate transition for
combining the gyro and ADS. The gyro data is not prefiltered as the ADS data is, the
G filter will serve this purpose. The requirements for the G filter are to provide unity
gain at low frequencies and to allow the ADS system to provide the gain at high
frequencies.
The G filter takes on the following values in the frequency domain:
Define:F1
= 2.0
F2 = 7.8125
freq = current frequency
Note: F1 was chosen as the frequency where Hads reaches 1.0 dB while F2 was
chosen to be the Nyquist frequency associated with the sampling of the gyro data.
if freq is less than F1 then Vk = 1
if freq is greater than F1 but less than F2 then
a = (F2  freq) * (F2  freq);
b = (freq  F1) * (freq  F1);
Vk = (a * Hpdf + b * PHa) / (a + b);
(Note: This is just an interpolation over the crossover region of the two systems.)
if freq is greater than F2 then
Vk = PHa
G(freq) = 0.0
Version 3.2
61
7/9/98
IAS Geometric ATBD Version 3
The G filter can then be solved for by looking at the response of the attitude
processing system before the F filter is to be applied [Reference 10].
GHg + PHa = Vk
The G filter can be solved for by writing this equation as:
(GHg + PHa) * conjugate{(GHg + PHa)} = Vk * Vk
Multiplying this out and noting that G is real so that G = conjugate{G} the equation
becomes:
G2Hg2 + 2GRe{ Hg conjugate{PHa} } + PHa2  Vk2 = 0
The quadratic equation can then be used to solve for G, the solution is chosen for the
minimum magnitude. This solution produces the following equations for the G filter.
Hc = conjugate{PHa}
Gpa = Hg * Hc
a = Hg * Hg;
b = 2.0 * Gpa.re;
c = PHa * PHa  Vk * Vk;
Re{G} = (b  (b * b  4.0 * a * c)1/2) / (2.0 * a);
Im{G} = 0.0;
Due to the fact that the gyro data is not prefiltered there is aliasing associated with
the measured signal. The G filter is applied in the spatial domain to try to remove
some of the aliasing at frequencies above 2 Hz and also to provide for a smooth
transition between the gyro and ADS transition region. Using the definition above of
the G filter in the frequency domain a spatial domain filter is achieved as follows.
Define:
Filter size = N
Padded filter size = M
Note: This explanation assumes that values for the filter are stored in an array. If
speaking in the frequency domain then the buffer location 0 in the array is the DC term
and location N/2 is the Nyquist frequency.
7.1 Determine an appropriate size for the G filter. Sehn and Miller state that a G
filter of length 31 guarantees less than 1% filter approximation error. It should
Version 3.2
62
7/9/98
IAS Geometric ATBD Version 3
also be noted that since the G filter is determined in the frequency domain it
makes sense to choose a value as close to a power of 2 as possible since this is
what the FFT routine needs. If a routine was used that did not need powers of 2
then this would not be an issue. The IAS uses a value of 63.
7.2 Find the nearest power of 2 larger then the associated filter size N. Zero pad the
data to this length (this will make the data of size M).
7.3 Step through the number of values chosen (filter length padded to a power of 2)
calculating the frequency response of the filter, the equations were given above.
This will be calculated at frequency increments that are dependent on the gyro
sampling rate and the number of samples from step 2. In practice, only one half
plus one of the filter size needs to be calculated due to step 7.4.
7.4 The G filter needs to be real in the spatial domain therefore the Fourier Transform
must be the symmetrical complex conjugate with respect to the origin, or for in
this case the Nyquist frequency. Reference 6 gives a good explanation on the Fast
Fourier Transform and some of its properties which should further explain this
process.
7.5 Calculate the inverse FFT of the results from step 7.4.
7.6 Adjust the data so that it is of the original desired filter length N. It should be
symmetrical about the center of the filter.
7.7 Apply the filter coefficients calculated in step 7.6 to the gyro data. When the
filter coefficients are applied to the ends of the gyro data the data is mirrored
about its end points.
8. Synchronize (up sample) the data from step 7 so that agrees with the 2 msec sampling
rate of the ADS data. The prototype uses a series of filters that allows better control
on the response of the overall interpolation (final sampling rate). A further
understanding and explanation of this process can be found in Reference 33. The
filtering coefficients used in the IAS are:
Num
Filter 1
Filter 2
Filter 3
Filter 4
Filter 5
Filter 6
Version 3.2
scale
1.0
2.0
16.0
32.0
256.0
346.0
h0
1.0
1.0
9.0
19.0
150.0
208.0
h1
0.0
0.0
1.0
3.0
25.0
44.0
63
h2
0.0
0.0
0.0
0.0
3.0
9.0
h3
0.0
0.0
0.0
0.0
0.0
0.0
h4
0.0
0.0
0.0
0.0
0.0
0.0
7/9/98
IAS Geometric ATBD Version 3
Filter 7
Filter 8
Filter 9
512.0
802.0
8192.0
302.0
490.0
5042.0
53.0
116.0
1277.0
7.0
33.0
429.0
0.0
6.0
116.0
0.0
0.0
18.0
These numbers are reflected around the point of interest (the point to be
interpolated). The filter coefficients are multiplied by the corresponding pixel it is
associated with and the sum is then divided by the scale factor. This will give an
interpolated point. For example, if interpolating a point between n and n+1 with filter
4 (pix(n) defines gray value at point n).
sum =
h1 * pix(n1) +
h0 * pix(n) +
h0 * pix(n+1) +
h1 * pix(n+2)
new pix = sum / scale[4]
The sampling start time and the time between samples for the smoothed low
frequency data is:
x,y,x axis
PCD + 64N N = 0,…255
while the ADS is sampled at increments of:
x axis PCD + 2N + 3/8 msec
y axis
PCD + 2N + 7/8 msec
z axis
PCD + 2N + 11/8 msec
where N = 0,…,8191
To get the smoothed low frequency data to have a time difference of 0.002 secs
between samples the data is interpolated 5 times with the coefficients listed above.
To account for the offsets in the ADS data, interpolation with the coefficients above
is performed 4 times. This will allow a time difference between samples to be chosen
such that every 16th sample will fall on an even increment of 0.002N secs from a PCD
cycle. Since there is no valid data at N=0 the first value for each axis of the time
synchronized ADS will be zero. The offset is accounted for by stepping ahead by a
set number of data points. The first valid synchronized ADS will occur at 0.002
secs.
(Note: Interpolating 4 times gives the ADS a time difference of 0.000125 secs)
xaxis (0.002  0.000375) / 0.000125 = 13 samples
yaxis (0.002  0.000875) / 0.000125 = 9 samples
xaxis (0.002 – 0.001375) / 0.000125 = 5 samples
Version 3.2
64
7/9/98
IAS Geometric ATBD Version 3
Every 16th sample is then chosen after taking into account the offset. This will return
the time difference between samples to 0.002 secs.
9. Sum the synchronized data from step 8.
10. Apply the F filter to the data obtained in step 9. Since the ADS data is prefiltered
and due to the fact that the gyro data has been G filtered, the F filter can be applied in
the frequency domain (aliased frequencies have been suppressed). Both Sehn and
Miller and EOSAT determine a spatial domain filter for the F filter in much the same
way that the G filter was determined. This approach would apply the F filter in the
spatial domain and would eliminate the need for the data in step 9 to be a power of 2.
The F filter takes on the following values in the frequency domain:
F(freq) = 1.0 / [PHa + G(freq)Hg]
The following steps were taken to apply the F filter.
10.1 Calculate the FFT of the data from step 9. It should be noted that if the FFT is
used then the number of data elements, when summed, must be a power of 2.
10.2 Step through the data from step 10.1 calculating the response of the F filter
given the definition above. This will be calculated at frequency increments that
are dependent on the ADS sampling rate and the number of samples from step
10.1.
10.3 Reflect complex conjugate of the data from step 10.2 around the Nyquist data
point (this is analogous to what was done for the G filter). This will assure that
the inverse FFT of this data set will be real.
10.4 Multiply the data from steps 10.1 and 10.3.
10.5 Calculate the inverse Fourier transform of the data
10.6 This data will be the roll, pitch, and yaw due to attitude perturbations as
measured from the gyro, gyro drift, quaternions and ADS of the satellite.
11. The combined data is rotated to the ACS.
(Details of steps 6 and 11 can be found in section 3.1.5.3.3, Convert ADS to a Common
Reference System)
Attitude Processing Sign Conventions
Version 3.2
65
7/9/98
IAS Geometric ATBD Version 3
The PCD quaternions are used to construct an orientation matrix which rotates the ACS
axes to ECI as described on page 39 of the ATBD. This is premultiplied by a matrix
which rotates ECI to Orbital coordinates, constructed from the ephemeris data, to create
an ACS to Orbital rotation matrix (body to orbit). This matrix is mapped to the roll,
pitch, yaw perturbation matrix:
cos( p) cos( y) sin( r ) sin( p) cos( y) + cos( r ) sin( y) sin( r ) sin( y) − cos( r ) sin( p) cos( y)
− cos( p) sin( y) cos( r ) cos( y) − sin( r ) sin( p) sin( y) cos( r ) sin( p) sin( y) + sin( r ) cos( y)
sin( p)
− sin( r ) cos( p)
cos( r ) cos( p)
From the last row and first column of this matrix:
roll = tan1(m32/m33) = tan 1(m32/m33)
pitch = sin1(m31)
yaw = tan1(m21/m11) = tan 1(m21/m11)
These are the body to orbit rotation angles.
Orbital motion is removed from the gyro data as follows. The ephemeris data is used to
construct the ECI to orbital rotation matrix at the first PCD major frame time
(ECI2ORB0). At each subsequent gyro point, the ephemeris is interpolated and used to
construct the instantaneous ECI to orbital rotation matrix which is inverted to form the
instantaneous orbital to ECI rotation matrix ORBt2ECI. These arrays are multiplied
(ECI2ORB0 ORBt2ECI) to form the rotation from the instantaneous orbital coordinate
system to the orbital coordinate system at the beginning of the PCD interval,
ORBt2ORB0. The same logic used above for the quaternions is used to extract roll, pitch,
and yaw values from this composite matrix yielding body to orbit rotation angles where
the “body” coordinate system is taken to be the ideal ACS which is perfectly aligned with
the instantaneous orbital coordinate system, relative to the reference orbital coordinate
system at PCD start. These “ideal” angles are subtracted from the gyro samples so that
the gyro data reflects deviations from the ideal orbital coordinate system, or body to orbit
rotation angles. These gyro body to orbit angles are combined with the body to orbit
angles derived from the quaternions, in the Kalman filtering procedure.
After the quaternion and gyro data are combined the first filtered attitude sample value is
subtracted from all of the samples and saved as a zero frequency (DC) reference. This DC
reference is sign inverted to form an orbit to body rotation reference. The residual gyro
samples are combined with the ADS data and, in the course of this filtering process, are
sign inverted so that the blended high frequency attitude data represent orbit to body
Version 3.2
66
7/9/98
IAS Geometric ATBD Version 3
perturbations. The orbit to body DC reference and the orbit to body high frequency
perturbations are added to form the net orbit to body rotation which is used in the
satellite model.
3.1.5.3.4
Line of Sight Generation and Projection
The algorithms that follow are used by the Geometric Correction Grid algorithm. The
inputs to this algorithm are line, sample, and band. This algorithm also uses the
information which was compiled and converted into a usable form in the Create Model
algorithm. The result of this algorithm is the geodetic latitude and longitude for a point on
the Earth, or a height h above the surface, which corresponds to the entered line, sample,
and band of the input image.
Calculate Sample Time
The scan line number is found by subtracting 1 from the line number and then dividing by
the number of lines per scan for the given band and then truncating.
The detector number is found by taking the decimal value remaining of the scan line
number, before truncating, and multiplying by the number of lines per scan.
N = number of lines per scan
LineNumber − 1
S = remainder
N
detector = N – floor(S*N)
The time into the current scan is found by:
For forward scans:
time_in_scan = (sample  1.0 + bandoffset) * DWELLTIME
For reverse scans
time_in_scan = (scan[*scan_line].counted_line_length  sample + bandoffset) *
DWELLTIME
For Landsat 7, the bands are aligned using left justified padding and so the padding must
be removed from the calculation as a function of band. The L0R scan line offset (SLO)
date are used to accomplish this.
Version 3.2
67
7/9/98
IAS Geometric ATBD Version 3
The time from start of image is found using the MSCD start of scan value for the current
scan minus the MSCD start of scan value for the first scan plus the time in scan.
Calculate Viewing Angles
Find the acrossscan angle of a detector due to its location within the band. Note: This
example is for the multispectral bands.
detector_angle = (LINES_SCAN / 2.0  detector + 0.5) * IFOVMS
Find the alongscan 5th order polynomial correction. For forward scans use the timeinscan value to determine the correction; for the reverse scans use the total scan time minus
the timeinscan value to determine the correction value. Do the same for the acrossscan
5th order correction.
The polynomial coefficients and scan start and stop angles are provided as measures of
mirror rotation. Since the instrument optical axis is being reflected off of the scan mirror
the effect on the line of sight in object space is twice the mirror angle. The scan mirror
model and the focal plane geometry must be combined to calculate object space alongscan
angles in the sensor coordinate system. The IAS sign convention is defined so that
positive scan angles are in the direction of the sensor Y axis (as described in section 3.1.2.
The alongscan angle for band “BAND” is calculated as a function of time “TIME” as
follows:
For forward scans:
s2m = A fsm = start to midscan angle
s2e = A fsm + A fme = start to end angle
along_angle = 2.0 * (s2m  s2e * TIME/Ts + f( b, TIME ))
+ bandoff_along[BAND]  odd_det_off[BAND]
where f( b, TIME ) is the corrected fifthorder polynomial model.
For reverse scans:
s2m = A rsm
s2e = A rsm + A rme
along_angle = 2.0 * (s2m + s 2e * TIME/Ts + r( b, T s  TIME ))
+ bandoff_along[BAND]  odd_det_off[BAND]
where r( b, TIME ) is the corrected fifthorder polynomial model.
The acrossscan angle is found using a linear interpolation of the Scan Line Corrector
Mirror's position, which is a function of time in scan. Also accounted for is the Scan
Mirror's acrossscan 5th order correction, the detector angle, and the acrossscan focalplane offset.
Version 3.2
68
7/9/98
IAS Geometric ATBD Version 3
Thus, the model for the SLC angle as a function of time in scan is:
slc_angle = slc_start[SCANDIR]  slc_max[SCANDIR] * TIME/Ts
where:
slc_angle = acrossscan pointing angle (positive in the direction of flight)
slc_start = SLC angle at the beginning of scan (indexed by scan direction)
SCANDIR = index for scan direction (forward or reverse)
slc_max = total SLC angle (indexed by scan direction)
TIME = time in scan
Ts = total scan time
The same equation applies for forward and reverse scans, with the appropriate selection
of the SLC start and SLC total angles. The sign convention and direction of SLC motion
is the same in either case: positive angles are in the direction of flight (toward the sensor
X axis) with the SLC rotating from a positive to a negative angle in each scan.
The SLC angle must be combined with the acrossscan detector angle and the scan mirror
acrossscan correction polynomial to compute the total acrossscan look angle for
detector N in band BAND:
For forward scans:
cross_angle = slc_start[FWD]  slc_max[FWD]*TIME/Ts
+ 2.0 * f( b, TIME ) + bandoff_cross[BAND]
+ ((ndets[BAND]+1)/2  N)*IFOV[BAND]
where f( b, TIME ) is the scan mirror acrossscan correction polynomial
evaluated at the time of interest and all other components are as defined
previously.
For reverse scans:
cross_angle = slc_start[REV]  slc_max[REV]*TIME/Ts
+ 2.0 * r( b, TIME ) + bandoff_cross[BAND]
+ ((ndets[BAND]+1)/2  N)*IFOV[BAND]
where r( b, TIME ) is the scan mirror acrossscan correction polynomial
evaluated at the time of interest. Note that the only difference between
forward and reverse scans is the selection of the appropriate SLC
parameters, as mentioned above, and the application of the scan mirror
acrossscan correction polynomial.
Version 3.2
69
7/9/98
IAS Geometric ATBD Version 3
Calculate LOS Vector
This algorithm uses the time in which a pixel was imaged to find the attitude, position,
and velocity of the spacecraft. The along and acrossscan angles are used to construct a
ideal line of sight (LOS) vector in sensor space. This LOS vector is transformed to the
spacecraft navigation reference base and is then transformed using the attitude information
to the Orbital Coordinate System (OCS). The OCS LOS vector is then transformed to
the Earth centered coordinate system.
1. To find the attitude, this algorithm uses the time from first ADS value to the time of
pixel imaging to look up the ADS and gyro combined attitude in the combined attitude
table. This table was constructed in the Process Attitude module. Added to the roll,
pitch, and yaw values from the table are the spacecraft's roll, pitch, and yaw values
(the DC reference values).
2. The “move satellite” algorithm uses the time from the reference ephemeris point to
the time of pixel imaging to interpolate a new, ECEF, satellite position and velocity.
The polynomials for this interpolation were defined in the Create Model module in
the Calculate Satellite State Vector algorithm.
3. The Find LOS algorithm uses the along and acrossscan angles to construct a LOS
vector in the sensor's coordinate system.
losx = sin(across_angle) * cos(along_angle)
losy = sin(along_angle)
losz = cos(across_angle) * cos(along_angle)
4. The Attitude algorithm transforms the LOS vector to the spacecraft's ACS coordinate
system. This vector is then transformed using the spacecraft's attitude to the OCS,
and then the vector is transformed to the Earth Centered Earth Fixed (ECEF)
coordinate system.
Transform the LOS from the sensor coordinate system to the ACS coordinate system:
3 x3
ACS
Sensor
= Transforma tion
LOS Matrix
LOS
The 3 x 3 transformation matrix from the sensor to the ACS is determined prelaunch.
As part of the IAS calibration activities, this matrix, from time to time will be checked
Version 3.2
70
7/9/98
IAS Geometric ATBD Version 3
and, if required, updated. The members of this matrix will be contained in the
Calibration Parameter File.
Transform the ACS LOS to the OCS using the attitude information by:
Orbital
−1
−1
−1 ACS
= [YAW ] [PITCH ] [ROLL ]
LOS
LOS
Transform the Orbital LOS to the ECEF coordinate system.
ECEF
LOS =
3 x3
Orbital
Transformation
LOS
Matrix
The 3X3 OCS to ECEF transformation matrix is found using the unit vectors of the
axes for the Orbital Coordinate System (OCS) determined in the ECEF coordinate
system. The Z axis of the OCS is defined as the negative of the spacecraft's position
vector. The Y axis of the OCS is the negative of the spacecraft's angular momentum
vector or (V X R). The X axis is the cross product of the Y axis with the Z axis (Y X
Z).
Intersect LOS Vector With Earth
This algorithm finds the geocentric location on the ellipsoid of the intersection of the LOS
vector to the ellipsoid or at a constant elevation for the entire scene above the ellipsoid.
The satellite's position vector, the LOS vector, and an elevation are the input to this
algorithm and the geocentric location of the intersection is output in both the Cartesian
and spherical coordinates.
This algorithm uses an iterative method to locate the intersection. The satellite's
geocentric latitude is used as the first guess in this process. This latitude is used to find
the geocentric radius of the ellipsoid plus the elevation at this latitude. Using the angle
between the satellite's position vector and the LOS vector, the radius of the Earth, and the
plane trigonometry's law of cosines, the magnitude of the LOS vector is found. The
intersection is the addition of the satellite's position vector and the scaled LOS vector.
However, this location was determined using a slightly incorrect latitude. Thus, the
calculated intersection’s geocentric latitude is used as the second guess. This process is
continued until the LOS magnitude difference between iterations is small (< 0.01 meter).
Another approach is to use the reduced latitude in the calculation. This allows for a direct
solution for a point on the surface of the ellipsoid. The solution is as follows:
Version 3.2
71
7/9/98
IAS Geometric ATBD Version 3
Define Rs as the position vector of the satellite, LOS as the line of sight vector and Rt as
the position vector of the target in question. The position vector of the target using the
reduce Latitude psi, its longitude Lt and the semimajor a and the semiminor b axes is:
Rt = X ti + Ytj + Z tk
Xt = a * cos (psi) * cos(Lt)
Yt = a * cos (psi) * sin(Lt)
Z t = b * sin (psi)
The solution requires the following equation to be true.
Rt = Rs + m * LOS
Where m is the magnitude of the final LOS vector. The current LOS vector has
a magnitude of 1.
or as separate components:
Xt = a * cos (psi) * cos(Lt) = Xs + m * Xlos
Yt = a * cos (psi) * sin(Lt) = Y s + m * Ylos
Z t = b * sin (psi) = Zs + m * Z los
A vector Rt' with a magnitude of 1 can be constructed by dividing the X t and Y t
components by a and dividing the Z t component by b.
Rt' = X t'i + Yt'j + Z t'k
Xt' = X t / a = cos (psi) * cos (Lt) = Xs / a + m * Xlos / a
Yt' = Yt / a = cos (psi) * sin (Lt) = Y s / a + m * Ylos / a
Z t' = Z t / b = sin (psi) = Zs / b + m * Zlos / b
Or
Rt' = Rs' + m * LOS'
Using the law of cosines for a plane:
Rt'2 = m2 * LOS'2 + Rs'2  2 * m * LOS' * Rs' * cos (w) = 1
Version 3.2
72
7/9/98
IAS Geometric ATBD Version 3
where w is the angle between the satellite's position vector and the line of sight vector.
Due to Rt' being defined to have a magnitude of 1, the only unknown left in the equation
is m.
Rearranging the equation to use the quadratic formula to solve for m results in:
0 = m2 * LOS'2  2 * m * LOS' * Rs' + Rs'2  1
and thus, the solution for Rt is:
Rt = Rs + m * LOS
The geocentric location is then converted to a geodetic location. The theory behind this
routine is found in Reference 20, page 398.
Speed of Light Correction
Due to the noninfinite velocity of light, the velocity of a point on the surface of the Earth
and the velocity of the satellite, cause pixel location errors if not accounted for. The speed
of light correction due to the Earth's rotation is submeter and is currently not accounted
for. The speed of light correction for due to the satellite's velocity is approximately a
constant for Landsat 7. The magnitude is about 15 meters and is in the along track
direction and is applied in the model.
3.1.5.3.5 Correction Grid Generation
The ETM+ Geometric Correction Grid defines the piecewise (gridded) geometric
transformation of imagery in satellite perspective to corrected imagery in a user specified
output projection of userspecified geographic extent. The grid is used by the geometric
transformation function (the Resampler) to rectify the image data to the corrected image
output space.
This algorithm assumes that the Landsat satellite model has been initialized using the
input image’s associated PCD and MSCD files. It is also assumes that the output
projection system is defined and read from a projection definition file. [Reference 19].
This algorithm also requires the frame (geographic extent) of the resulting output space
Version 3.2
73
7/9/98
IAS Geometric ATBD Version 3
image, and the dimension of each pixel in the output coordinate space. Note: The PCD
and MSCD files are used by the Landsat satellite model to define the input image's (0R or
1R) geometric space.
Note: Refer to the LAS Geometric Manipulation Overview Document [Reference 18] for
more information on specific projection parameters for a given projection, and for
information on the Projection Transformation Package, a coordinate transformation
subroutine library based on the U.S. Geological Survey's General Cartographic
Transformation Package (GCTP).
Defining the Frame
The first step in the gridding process is to determine the geographic extent of the output
image to be generated by the resampler. This geographic extent of the output image space
is referred to as the output space “frame”, and is specified in output image projection
coordinates. There are four different methods that are used to determine the output frame.
Method 1:
The user defines the upperleft and lowerright corner coordinates of the area of
interest in geographic (latitude/longitude) coordinates. These coordinates are then
projected to the output projection coordinate system using the Projection
Transformation Package. This usually results in a nonrectangular area so a
minimum bounding rectangle is found (in terms of minimum and maximum X and
Y projection coordinates) in the resulting output space. This minimum bounding
rectangle defines the output space frame. The output image pixel size is then
applied to the projection space to determine the number of lines and samples in
the output space.
Method 2:
The user defines a rectangular area in the output space by specifying upperleft
and lower right output space projection coordinates. The output image pixel size
is then applied to the projection space to determine the number of lines and
samples in the output space.
Method 3:
The user defines an upperleft output space projection coordinate, the output
image pixel size, and the number of lines and samples in the output image. The
coordinate pair must be the minimum X and maximum Y projection coordinate.
The maximum X and minimum Y coordinates (the other end of the bounding
rectangle) are then calculated.
Method 4:
Version 3.2
74
7/9/98
IAS Geometric ATBD Version 3
The user allows the framing software to project the 4 corners of the input L0R
data set to the Earth ellipsoid. This will produce the minimum bounding rectangle
which contains all of the input image data.
Method 5:
The user specifies a path oriented Landsat product in either the Space Oblique
Mercator (SOM) or Universal Transverse Mercator (UTM) projection. In this
case, the framing coordinates are not userspecified. The frame is a preset number
of lines and samples based on the Landsat Worldwide Reference System (WRS)
scene size and the maximum rotation needed to create a path oriented product.
For a pixel size of 30 meters, a rectified image in the SOM or UTM projection is
approximately 6440 lines by 6850 samples in size.
For diagrams and more information on the framing process, refer to the LAS Geometric
Manipulation Package Overview Document [Reference 18].
Path Oriented Framing and Rotation
A path oriented projection is basically a rotation from a typical "mapnorth" up
projection that better represents the near polar orbit of the Landsat satellite (nominal
inclination angle of 98.2 degrees). The path oriented projection also does a better job of
maintaining the scan line to the output space line.
The first step in generating a path oriented projection is to calculate the center location in
the output space that the frame is to be rotated about. This is done by computing the
nominal WRS scene center for the current path and row. The frame is centered in the
output projection space at this projection coordinate and will be rotated about this point.
The procedure for computing the WRS scene center and the associated rotation angle is:
Given:
Parameter
Earth Parameters
SemiMajor Axis
SemiMinor Axis
L7 Orbit
Parameters
Inclination
L7 WRS
Parameters
Version 3.2
Value
Source
6378137.0 m
6356752.314 m
CPF:Earth_Constants:Semi_Major_Axis
CPF:Earth_Constants:Semi_Minor_Axis
98.2 degrees
CPF:Orbit_Parameters:Inclination_Angle
75
7/9/98
IAS Geometric ATBD Version 3
Number of Paths
Number of Rows
WRS Repeat Cycle
Desc Node Row
Long of 001/060
Input Parameters
Path
Row
233
248
16 days
60
64.6 degrees
CPF:Orbit_Parameters:WRS_Cycle_Orbits
CPF:Orbit_Parameters:Scenes_Per_Orbit
CPF:Orbit_Parameters:WRS_Cycle_Days
CPF:Orbit_Parameters:Descending_Node_Row
CPF:Orbit_Parameters:Long_Path1_Row60
Integer 1233
Float 0 < R <
249
User
User
Find:
WRS Scene Center Latitude
WRS Scene Center Longitude
Scene Heading
Method:
Convert input angles to radians:
Inclination_Angle = Pi / 180 * Inclination_Angle
Long_Path1_Row60 = Pi / 180 * Long_Path1_Row60
Compute the Earth’s angular rotation rate:
earth_spin_rate = 2 * Pi / (24 * 3600)
Note: We use the solar rotation rate rather than the sidereal rate as called for in the
General Electric Landsat D PIR which describes the Worldwide Reference System, in
order to account for the orbital precession which is designed to make the orbit sun
synchronous. Thus, the apparent Earth angular velocity is the inertial (sidereal) angular
velocity plus the precession rate which, by design, is equal to the solar angular rate.
Compute the spacecraft’s angular rotation rate:
SC_Ang_Rate = 2 * Pi * WRS_Cycle_Orbits / (WRS_Cycle_Days * 24 * 3600)
Compute the central travel angle from the descending node:
Central_Angle = (Row  Descending_Node_Row)/Scenes_Per_Orbit * 2 * Pi
Compute the WRS geocentric latitude:
WRS_GCLat = asin( sin(Central_Angle) * sin(Inclination_Angle) )
Compute the longitude of Row 60 for this Path:
Long_Origin = Long_Path1_Row60  (Path  1) * 2 * Pi / WRS_Cycle_Orbits
Compute the WRS longitude:
Version 3.2
76
7/9/98
IAS Geometric ATBD Version 3
Delta_Long = atan2( tan(WRS_GCLat)/tan(Inclination_Angle),
cos(Central_Angle)/cos(WRS_GCLat) )
WRS_Long = Long_Origin  Delta_Long  Central_Angle * Earth_Spin_Rate /
SC_Ang_Rate
Make sure the longitude is in the range +/ Pi:
While ( WRS_Long > Pi )
WRS_Long = WRS_Long  2*Pi
While ( WRS_Long < Pi )
WRS_Long = WRS_Long + 2*Pi
Compute the scene heading:
Heading_Angle = atan2( cos(Inclination_Angle)/cos(WRS_GCLat),
cos(Delta_Long)*sin(Inclination_Angle) )
Convert the WRS geocentric latitude to geodetic latitude:
WRS_Lat = atan( tan(WRS_GCLat) * (Semi_Major_Axis/Semi_Minor_Axis)
* (Semi_Major_Axis/Semi_Minor_Axis) )
Convert angles to degrees:
WRS_Lat = WRS_Lat * 180 / Pi
WRS_Long = WRS_Long * 180 / Pi
Heading_Angle = Heading_Angle * 180 / Pi
Round WRS lat/long off to the nearest whole arc minute:
WRS_Lat = round( WRS_Lat*60 ) / 60
WRS_Long = round( WRS_Long*60 ) / 60
Return the answers:
WRS_Lat
WRS_Long
Heading_Angle
The WRS heading angle, from geodetic north, will need to be converted from the WRS
scene center to a frame orientation angle in map coordinates. The can be done by first
converting the WRS scene center lat/long to map projection X1, Y1 and moving to a point
slightly north of the WRS scene center. This is done by adding 1 microradian
(0.2seconds) to the WRS scene center latitude and projecting this point to X2, Y2. Next
the azimuth of this line is computed in grid space as the arctangent of the X difference
(X2X1) over the Y difference (Y2Y1). This is the grid azimuth of geodetic north at the
WRS scene center. Adding this angle to the geodetic heading to give the grid heading. A
standard framed scene puts the satellite direction of flight at the bottom of the scene so
the scene orientation angle is the grid heading plus or minus 180 degrees. If the grid
Version 3.2
77
7/9/98
IAS Geometric ATBD Version 3
heading is less than zero add 180 degrees. If the grid heading is greater than zero subtract
180 degrees. This is the scene orientation angle to be used with the WRS scene center and
the fixed standard scene size to calculate the frame corners.
It should be noted that the preceding approach does not take into account the skewing of
the image data. Skewing of the image is due to the Earth’s’ rotation relative to the
spacecraft’s absolute velocity.
A delta X and delta Y are calculated from the center of the projection space to the four
corners of the output frame area using the output pixel size and the number of lines and
samples in the output space.
Example: Calculation of the upper right delta X and delta Y
delta_x = (number_samples  center_sample) * pixel_size
delta_y = (center_line  1.0) * pixel_size
The deltas and the scene orientation angle are then used to calculate the rotated projection
coordinate of each corner.
ur_proj_x = delta_x * cos(angle) + delta_y * sin(angle) + center_proj_x
ur_proj_y = delta_x * sin(angle) + delta_y * cos(angle) + center_proj_y
A least squares routine is used to generate a first order polynomial which maps output
image (line, sample) coordinates to rotated output projection coordinates.
The output pixel size and projection coordinates are saved in the mapping grid file.
Gridding The Input Space
For the process of rectifying the ETM+ satellite imagery to a given projection, a regular
grid is defined over the input space, one set per scan. This separate, forward mapping is
necessary to properly correct for gaps between scans and because an inverse satellite
model cannot be defined. The first and last lines from each scan are represented in the
gridding of the input. The difference between the last line of a scan (line 16 for TM multispectral bands, line 8 for the thermal band, and line 32 for the Pan band) and the first line
of the next scan, is used to determine the scan gap. These scan gaps are corrected in the
ETM+ geometric resampling process. Refer to the Landsat 7 Wide Band Data Format
Book (Reference 1) for more information about the Landsat 7 instrument and scan gap.
The grid spacing in the alongscan (sample) direction 96 pixelsa factor that divides into
the total number of input image samples and that coincides with the ADS sampling rate.
Version 3.2
78
7/9/98
IAS Geometric ATBD Version 3
(The typical value used in for Landsat 4 and 5 was 128 pixels.) Figure 3.1.511 illustrates
the gridding of an input image:
Figure 3.1.5 11 Input Image Gridding Pattern
Relating Input to Output Space
The Landsat 7 mapping grids are defined in image (line, sample) coordinates which map
input (level0) line, samples to output, projected space line, samples. At each grid
intersection (input line and sample) the initialized Landsat satellite model is called to
calculate the geometric latitude and longitude at that point. Using the projection definition
information as input to the Projection Transformation Package, the output projection
coordinate is then calculated for each latitude and longitude point. From the output
projection coordinates an output line and sample location, in the defined output frame, is
calculated for each point. For a path oriented product, the polynomial is evaluated at each
output projection coordinate. This mapping takes care of both the rotation of the
projection space and the scaling from projection coordinates to image (line, sample)
coordinates. The polynomial (called once for x coordinates and once for y coordinates) is
of the form:
x' = a0 + a1x + a2y
where: x' is the output line or sample and
x and y are output space projection coordinates
For the other output projections (when there is no rotation between the image coordinate
system and the projection coordinate system) the output line and samples are calculated
directly using the following equations:
out_line = (ul_proj_y  proj_y) / pixel_size + 1
out_samp = (proj_x  ul_proj_x) / pixel_size + 1
Version 3.2
79
7/9/98
IAS Geometric ATBD Version 3
The input line and sample, output line and sample, and the projection information are
written to a grid file. This geometric mapping grid file defines the transformation between
the satellite image and the specified output projection. The grid file is used by the
resampler to geometrically rectify the input image to the userspecified output image
frame, projection, and pixel size.
The above process is done to create a grid for each band that is to be processed. The net
effect is that there is a grid for each band in the grid file. This is necessary because of
timing offset differences between bands.
3.1.5.3.6 Resampling
This section describes the IAS image resampling algorithm.
Calculate Scan Gap
This algorithm will calculate the scan gap and the sample shift between each grid cell of
two consecutive scans. This information will be used to create a synthesized or extended
image that will fill any gaps or "missing data" between two scans. A scan gap will refer to
a line offset between consecutive scans that is larger than the nominal scan line spacing
while a misalignment will refer to sample offset between consecutive scans in the image.
A least squares routine is used in this algorithm to calculate mapping coefficients. The
algorithm follows:
1. Loop through each grid cell calculating the scan gap and misalignment. This is done
using the mapping equations for each grid cell:
Forward mapping equations:
So = a0 + a1 * Si + a2 * Li + a3 * Si * Li
Lo = b0 + b1 * Si + b2 * Li + b3 * Si * Li
Where
Si is the input sample
So is the output sample
Li is the input line
Lo is the output line
a0,a1,a2,a3,b0,b1,b2,b3 are the bilinear mapping coefficients.
Inverse mapping equations:
Version 3.2
80
7/9/98
IAS Geometric ATBD Version 3
Si = c0 + c1 * So + c2 * Lo + c3 * So * Lo
Li = d0 + d1 * So + d2 * Lo + d3 * So * Lo
Where
Si is the input sample
So is the output sample
Li is the input line
Lo is the output line
c0,c1,c2,c3,d0,d1,d2,d3 are the bilinear mapping coefficients.
The four corners of each grid cell are used to solve for the mapping coefficients.
By using the forward and inverse mapping equations for each grid cell of a scan along
with the mapping equations for the grid cell of the scan directly below, the maximum
gap is obtained. (If scan N cell M is being processed then it is compared with scan
N+1 cell M). This is shown in figure 3.1.5 12.
Also during this time the mapping coefficients are calculated (forward and inverse) for
each grid cell and the inverse coefficients are stored in the grid structure.
These values are obtained from the following steps. (Scan N cell M is used to
designate one grid cell):
1. Loop through each grid cell of each scan, except for the last scan since this scan
cannot have a gap. Repeat steps 28 for number of grid cells.
2. Create a fictitious "extended" pixel 17E located directly (one line length) below
scan N cell M.
3. Find the output to input (inverse coefficients) and input to output (forward
coefficients) mapping coefficients for scan N cell M.
4. Store the inverse mapping coefficients into the grid structure.
5. Map pixel 17E to the output space using the forward mapping equations for
scan N cell M.
6. Find the output to input mapping coefficients for scan N+1 cell M.
7. Map the output pixel created in step 5 back to the input space using the
inverse mapping equations for scan N+1 cell M (call this point xg,yg).
Version 3.2
81
7/9/98
IAS Geometric ATBD Version 3
8. The difference between (xg,yg) and sample one line one of scan N+1 cell M
represents the misalignment and the scan gap for scan N cell M. See figure
3.1.512.
9. Find inverse mapping coefficients for grid cells of last scan and store these
coefficients into the grid structure.
2. The maximum scan gap in the image is found. The maximum scan gap and the
resampling kernel size will determine the size of the intermediate extended image.
These two parameters will determine the number of gap, or extended, lines needed
(gap_lines).
A gap of greater than 1/32 of a pixel is considered large enough to add one extra line to
each scan. This number represents the pixel increments of the resampling kernel.
The image must also be extended according to the resampling kernel size. A value of
kernel size minus one will make it possible to only need one scan of data at a time
during resampling. If this was not done then at the bottom of each scan the
resampling kernel would step into the scan directly below. Due to scan gap and
misalignment this data would need to be resampled to "fit", evenly spaced pixels, the
current scan being processed. To avoid a redundancy of this step it is only done once
up front during extended image generation.
See Reference 26 and Reference 27 for more information about the characteristics of scan
gap. Refer to Reference 28, Reference 29, and Reference 18 for more information about
geometric corrections and transformations.
Creating the Extended Image
References 26 and 27 describe scan gap correction. It should be noted that these
documents refer to LANDSAT 4 and 5. These documents also follow the traditional
three pass approach for extending the gap and creating an output image.
Each scan is extended by both the maximum number of lines missing due to the scan gap
and due to the resampling kernel size.
Algorithm:
1. Create weight tables. A set of cubic convolution weights are set up with 1/32 pixel
increments. This table will be used to resample the data in the sample direction. The
weights are set up for a shift of 0.0 to 1.0 pixels at 1/32 increments. A set of spline
weights are set up with 1/32 pixel increments. This table will be used to resample in
the line direction. The weights are set up for a shift of 0.0 to the maximum gap
Version 3.2
82
7/9/98
IAS Geometric ATBD Version 3
possible from the system. A value of 3.0 covers the maximum gap possible in the
LANDSAT 5 satellite. A value of 6.0 will be needed for the LANDSAT 7 satellite,
to cover the PAN band.
2. Loop through the image creating extra image data after each grid cell of every scan,
except for last scan which does not need to be extended. (For the processing of scan
N cell M); Repeat steps 36 for appropriate number of extended lines needed.
3. Get misalignment and scan gap for the current grid cell.
4. Determine data needed.
a) If the gap is zero then the lines are extracted only from the scan below and
resampled to adjust for scan misalignment and detector delays.
b) If the scan gap is greater than zero then the first extended line is created from
the last two lines from scan N and first two lines from scan N+1. These
lines are used as long as an extended line falls in between their line locations.
When the extended line starts to fall into scan N+1, going beyond line one of
scan N+1, then the lines used must be incremented appropriately. This can
be explained by looking at figure 3.1.512. If x1=1.0, x2=2.0, x3=4.5, and
x4=5.5 corresponds to the line locations of the last two lines in scan N and
the first two lines in scan N+1 respectively, then extended lines will be
created at 3.0 and 4.0 using these four lines. In this case both extended lines
fall between these four lines. However, when an extended line is created at
5.0 a new set of lines need to chosen such that x1=3.0, x2=4.0, x3=5.5, and
x4=6.5. At this point the first two lines are the previously created extended
lines while the other two lines used will be from scan N+1 and will
correspond to lines two and three. This pattern will then follow with two
extended lines being used along with two lines from scan N+1 for all other
lines created.
c) If the scan gap is less than zero then the first extended line is created with
two lines from the scan N and scan N+1. The lines taken correspond to the
last two from scan N and two from scan N+1 incremented by the
appropriate offset. If the scan gap is −0.5 then lines 15 and 16 from scan N
would be used along with lines 2 and 3 from scan N+1. After the first
extended line is created it is used for further processing and only one new line
from scan N+1 is needed. For the example just given the second iteration
would include line 16 from scan N, the extended line created in the previous
iteration, along with lines 3 and 4 from scan N+1.
Version 3.2
83
7/9/98
IAS Geometric ATBD Version 3
5. The lines taken from scan N are resampled to adjust for detector delays. The data
will be resampled using cubic convolution for a constant shift across the whole line,
where the shift is equal to the detector delay of that line.
6. The lines taken from scan N+1 are resampled to adjust for detector delay and
misalignment. The data will be resampled using cubic convolution so that it is scaled
and shifted to account for geometric misalignment and scaling differences between the
scan N cell M and scan N+1 cell M. It will also take into account detector delays.
The transformation is treated as a linear relationship, mapping the lines from scan
N+1 so that the samples line up in the sample direction.
This is done by the following relationship
S1N1 = seg * ns + 1
S2N1 = (seg+1)*ns + 1
S1N2 = S1N1  misalignment[seg];
S2N2 = S2N1  misalignment[seg+1];
m = (S2N2  S1N2) / (S2N1  S1N1);
b = S1N2  S1N1 * a1 + delay;
new_sample = m * old_sample + b;
where
S1N1 = first sample from scan N
S2N1 = second sample from scan N
S1N2 = first sample from scan N + 1
S2N2 = second sample from scan N + 1
m and b = mapping coefficients
out_sample = newly mapped sample
in_sample = sample from scan N+1
delay = detector delay for line from scan N+1
seg = current scan cell
misalignment = scan misalignment between two cells
ns = length of a cell
7. A new line is produced using the four lines created from steps 4 and 5. This is done
using a cubic spline which will take into account the unequal spacing between pixels
due to scan gap. The scan gap is treated as a linear function over the scan cell.
This is done by the following relationship
m = (gap_at_end  gap_at_start) / ns;
Version 3.2
84
7/9/98
IAS Geometric ATBD Version 3
b = gap_at_start;
new_gap = m * sample + b
where
gap_at_start = size of scan gap at the start of grid cell
gap_at_end = size of scan gap at the end of grid cell
m and b = mapping coefficients
ns = length of a cell
sample = sample position of current resampling step. These
values are from 0 to SAMP_RATE  1.
new_gap = gap associated with sample
See figures 3.1.514 and 3.1.515.
Notes on the Cubic Spline
The cubic spline weights used were developed by Fischel [Reference 30] and take on the
form of (See figure 3.1.513):
w1 = (U2/(3*h1*h2*h2))*(2*FG);
w2 = (x3x)/h2(F/(3*h1*h2*h2*h2))*(2*U2*(h1+h2)+h1*V2)
+(G/(3*h1*h2*h2*h2))*(U2*(h1+h2)+2.0*V2*h1);
w3 = ((xx2)/h2)+(F/(3*h2*h2*h2*h3))*(2*U2*h3+V2*(h2+h3))
(G/(3*h2*h2*h2*h3))*(U2*h3+2*V2*(h2+h3));
w4 = (V2/(3.0*h2*h2*h3))*(2.0*GF);
where
h1 = x2  x1;
h2 = x3  x2;
h3 = x4  x3;
F = (x3x)*((x3x)*(x3x)h2*h2);
G = (xx2)*((xx2)*(xx2)h2*h2);
U2 = h2/(h1+h2);
V2 = h2/(h2+h3);
Figure 3.1.516 shows the spline weights calculated for different size gaps. It should be
noted that at larger gaps the weights show greater variations. When the PAN band is used
for LANDSAT 7 this could cause some anomalies in cases of large varying grey values
between pixels that are located in an area with a large scan gap.
The PseudoInverse Model
Version 3.2
85
7/9/98
IAS Geometric ATBD Version 3
Due to the fact that there is no inverse model a method is needed to find which output
pixel correspond to which input pixel. This is done using the inverse mapping
coefficients for each grid cell. To do this it must be determined which grid cell goes with
what output pixel. A "rough" polynomial is set up which gives an approximate input
pixel for a given output pixel. An approximate grid cell location can then be determined
from this pixel. After this approximation is performed the correct grid cell can be found
using a trial and error process with the inverse mapping coefficients for the grid cells.
Calculating the Rough Polynomial
A first order polynomial was used for the "rough" mapping polynomial. The data used to
produce this polynomial was taken from the grid input and output pixel locations. The
reasoning behind using the first order was mostly from a speed and computational stand
point. The algorithm is as follows:
1. Read input and output lines/samples (not all are needed).
2. Use least squares routine to calculate a first order polynomial.
This creates a "rough" polynomial.
Si = a0 + a1 * So + a2 * Lo + a3 * So * Lo
Li = b0 + b1 * So + b2 * Lo + b3 * So * Lo
Where
Si is the input sample
So is the output sample
Li is the input line
Lo is the output line
a0,a1,a2,a3 are the polynomial mapping coefficients.
b0,b1,b2,b3 are the polynomial mapping coefficients.
Finding the Corresponding Grid Cell
The algorithm to determine the corresponding input grid cell for a given output grid cell is
as follows: (Note: If the inverse mapping routine is caught or jumping back and forth
between two scans then it is assumed that the output pixel falls in the gap between the
two scans. In this case the scan above the gap is passed back as the appropriate grid cell.)
Map output pixel to current_input_space_point using “rough” polynomial.
Set current_grid_cell from the current_input_space_point.
Loop
Version 3.2
86
7/9/98
IAS Geometric ATBD Version 3
Calculate the new_input_space_point using the inverse polynomial for the
current_grid_cell.
Calculate the new_grid_cell from the new_input_space_point.
If
current_grid_cell and new_grid_cell are the same correct grid cell has been found,
exit loop.
Else
Copy new_grid_cell to current_grid_cell.
Save current_grid_cell to history list.
If
Current_grid_cell is the same as the grid cell from previous iteration.
If
Matching grid cell is from the two iterations ago and the intervening grid cell
was from the same row, use last grid cell as the correct grid cell.
Else
Set flag to indicate the output pixel maps to the gap, use the grid cell in the
cycle from the upper most row as the correct cell.
Exit loop.
See figures 3.1.517, 3.1.518, and 3.1.519.
The Resampling Process
The output image is created by using the extended image and the image grid. Currently a
two dimensional cubic convolution kernel is used for resampling.
This processing flow may need to be revised if a new weight table is needed due to
detector delays or due to the MTF. This process does not contain any precision or
terrain processing concerns.
The algorithm is as follows:
1. Create resampling table weights. (Refer to the section on detector delays [section
3.1.5.3.6.1] for more information).
2. Step through the output space (for each pixel); Repeat steps 3  8 for all output
pixels
3. Use rough polynomial to map an output pixel to input space.
4. Use the “Find Corresponding Grid Cell” algorithm to determine which input image
grid cell the output pixel falls in along with the input pixel location for the
corresponding output pixel.
5. Determine (and get) data needed in extended image.
6. Determine correct resampling weights.
Version 3.2
87
7/9/98
IAS Geometric ATBD Version 3
•
If cubic convolution or MTF compensation is being used, find the fractional
portion of the input line/sample to look up the correct set of resampling
weights.
• If nearest neighbor resampling is being used then find the
nearest input pixel.
7. Perform resampling; Multiply each input pixel by the corresponding weight and
sum the results.
8. Store output pixel.
See References 12, 27, and 28 for more information about the resampling process.
Figure 3.1.5 12 Extended Pixels and Scan Alignment
Version 3.2
88
7/9/98
IAS Geometric ATBD Version 3
Figure 3.1.5 13 Calculation of Scan Gap
Version 3.2
89
7/9/98
IAS Geometric ATBD Version 3
Figure 3.1.5 14 Scan Misalignment and Gap
Figure 3.1.5 15 Extended Scan Lines
Version 3.2
90
7/9/98
IAS Geometric ATBD Version 3
Figure 3.1.5 16 Cubic Spline Weights
Figure 3.1.5 17 Inverse Mapping with “Rough” Polynomial
Version 3.2
91
7/9/98
IAS Geometric ATBD Version 3
Figure 3.1.5 18 “Rough” Polynomial – First Iteration
Figure 3.1.5 19 Results Mapped back to Input Space
Figure 3.1.5 20 Nearest Neighbor Resampling
Version 3.2
92
7/9/98
IAS Geometric ATBD Version 3
Nearest Neighbor Resampling
To perform nearest neighbor resampling, changes to the algorithm described above must
be made. To get the true nearest neighbor the scan gap is not extended. This presents
two scenarios; one case where the pixel falls within a grid cell and the resampling follows
the traditional nearest neighbor routine and the second case where the pixel falls in the
scan gap area. In the second scenario there is extra work that must be performed to
determine which pixel is really the true nearest neighbor. The confusion in determining
the true nearest neighbor lies in the fact that each grid cell defines its own coordinate
system. When a pixel falls within the gap the transformations between both grid cells
must be used to determine which pixel is truly closest for nearest neighbor resampling.
This is demonstrated in figure 3.1.520.
Linear Interpolation Resampling
Due to the effects of detector delays the linear interpolation employed is slightly different
then the classical bilinear approach. Detector delay are shifts in the sample direction with
respect to the raw imagery. When these offsets are taken into account during linear
interpolation, the sample actually represents the centroid of the polygon that is made up
of the four samples that surround the output sample location. These values are solved
directly, there is not a table of values representing a set of offsets as there is for cubic
convolution resampling.
Notes on the Correction Grid
Input lines in the grid must contain the top and bottom lines of a scan (for scan N line 1
and 16 must be the contained in the input lines of the grid). This is made necessary by
the way the scan gap is calculated. This algorithm expects the sample spacing between
grid points to evenly divide into the total number of samples in the image. It should also
be noted that the pixel locations given in the grid refer to the pixel center location. See
Section 3.1.5.3.5, Correction Grid Generation, for more information on the Gridding
process.
Notes on the Least Squares Routine
This algorithm uses the QR decomposition functions which are based on information in
References 6 and 31.
3.1.5.3.6.1 Detector Delays
A detector delay represents the virtual distance of each detector from its respective band
center. It contains the net effect of detector misplacement on the focal plane and the
Version 3.2
93
7/9/98
IAS Geometric ATBD Version 3
effective sampling delays due to electronics and scan angle monitor delays. The manner in
which they are defined is illustrated in Figure 3.1.521.
Since the delays are different for each detector of each band, they cannot be centered in
the sparse geometric correction grid. Therefore, the detector delay information has been
incorporated into the weight tables for resampling. The way in which they have been
applied assumes that the nominal position of the detectors has been used in determining
the alongscan angle in the forward satellite model ( Figure 3.1.521). It is the difference
between the real and ideal detector positions that is applied to the input sample when
determining the resampling weight table. Following this method, the difference is added to
the input samples of forward scans and subtracted from reverse scans. In order to
simplify the application of detector delays in calculating the resampling weights, the
minus sign has been included in the definition of the delay correction term for the reverse
scan so that it can be added regardless of scan direction. So:
new sample = old sample + delay correction
where:
 (delay  nominal) : Forward
delay correction = 
 (delay  nominal) : Reverse
Figure 3.1.5 21 Detector Delay Definition
Detector delays are defined as displacements relative to band center.
Version 3.2
94
7/9/98
IAS Geometric ATBD Version 3
Prior to resampling, a set of weights are determined which will incorporate the nonideal
locations of the detectors. In Figure 3.1.522, the ideal line and sample locations are
shown as dotted lines. When an output (line, sample) pair is mapped to input space, it
generally falls on a noninteger (line, sample). Typically, the surrounding integer pixels are
used to generate weights for resampling (i.e., cubic convolution, MTF, bilinear, etc.). To
incorporate detector delays, it is necessary to determine the weights based on the distance
from the noninteger (line, sample) in input space to the actual detector location (see
Figure 3.1.522). Weights for the synthetic “detectors” in the extended lines are
calculated under the assumption that they act as ideally placed detectors.
The standard cubic convolution uses a foursegment, piecewise third order polynomial
over the range [2, +2] [see Reference 12]. It is defined so that beyond [2, +2] the
weights are zero. From Fig. 3.1.522, note that if the detectors were ideally located at S1...S+2, then it would be necessary to calculate weights at these four samples only (for
each of the four lines L1...L+2). However, since the detectors can be displaced from their
nominal positions in either the positive or negative direction, it is necessary to calculate
weights over the range S2...S+3. For every line, two of these weights will be zero, but
which two weights are zero will depend on the magnitude and sign of the displacement
and the location of the noninteger (line, sample).
Version 3.2
95
7/9/98
IAS Geometric ATBD Version 3
Figure 3.1.5 22 Resampling Weight Determination
Resampling weights are determined based on the x and y distances between the
noninteger (line, sample) and the real detector locations
3.1.5.3.7 Precision Correction
1. Defining the Precision Correction Problem
The geometric measurement in ETM sensor system can be regarded as the look vector,
l sc, in the spacecraft bodyfixed system. This vector is transformed into the
Version 3.2
96
7/9/98
IAS Geometric ATBD Version 3
Figure 3.1.5 23 Definition of Orbit Reference System
Orbit Reference Frame [Reference 2] (OB) system, [see Figure 3.1.523] through the
spacecraft attitude parameters:
l ob = TT(ξr, ξp, ξy) l sc.
(1.1)
where ξr, ξp, and ξy are roll, pitch, and yaw angles, T is the transformation matrix, and can
be expressed as
T (ξr, ξp, ξy) = R3(ξy)R2(ξp)R1(ξr)
cosξ cos ξ
cosξ r sinξ y + sinξ r sinξ p cosξ y sinξ r sinξ y  cosξ r sinξ p cosξ y
p
y
=  cosξ p sinξ y cosξ r cosξ y  sinξ r sinξ p sinξ y sinξ r cosξ y + cosξ r sinξ p cosξ y
sinξ p
 sinξ r cosξ p
cosξ r cosξ p
the approximation for small angles ξr, ξp , ξy
1
ξy  ξp
=  ξ y 1 ξ r
ξ p  ξ r 1
Version 3.2
(1.2)
97
7/9/98
IAS Geometric ATBD Version 3
where R1, R2, and R3 are the coordinate system transformation matrix for rotation around
x, y and z axis respectively.
The vector l ob is further transformed into the Earthfixed (EF system) system
l ef = Tf(ref , vef ) l ob
(1.3)
where Tf is the forward transformation for vectors from the OB system to the EF
system, as a function of the satellite position ref and velocity vef vectors in the EF
system. Note that vef should be the "inertial" velocity expressed in the EF system as
described by Reference 13. This vector, l ef, together with the satellite position vector, ref ,
is then used to intersect the ellipsoid Earth surface to define a point position, Ref , as the
target point on the Earth. This is the common forward image pixel geolocation calculation.
Mathematically, Ref is a function of l sc, ξr, ξp, ξy, ref , and vef .
Ref
=
F(l sc, ξr, ξp, ξy, ref , vef )
(1.4)
Because of errors in the satellite orbit ephemeris and attitude data, this calculated Ref is
different than the true location of the image pixel. If we know the true location of a
landmark pixel, Rcp , from other sources (i.e., base map, survey etc.), this point can be
taken as a Ground Control Point (GCP) to check the accuracy of the computed image
pixel location. The precision correction process uses the Ground Control Point
coordinates to estimate the correction to the satellite ephemeris and attitude data, so that
with the corrected parameters in equation (1.4) the calculated image pixel location, Ref ,
will be close to its true location, Rcp , (depending on the GCP positional accuracy).
To calculate the precision correction, the difference between Ref and Rcp is taken as the
observable, and the observation equation becomes:
dR = Rcp

F(l sc, ξr, ξp, ξy, ref , vef )
(1.5)
according to equation (1.4). However, the actual calculation of Ref is usually not an
explicit function of the orbit and attitude parameters, especially for the intersecting
procedure. Therefore it is inconvenient to linearize equation (1.5) with standard
estimation techniques. Instead, the calculation of look vector l cp corresponding to Rcp , in
the OB system, is much more explicit:
l cp = Ti (ref , v ef )
Version 3.2
(R
cp
− ref
)
(1.6)
R cp − ref
98
7/9/98
IAS Geometric ATBD Version 3
where (Rcp  ref ) is the lineofsight vector in the EF system corresponding to Rcp , and
Ti(ref , vef ) is the inverse transformation for the look vector from the EF system to the OB
system. If all the satellite attitude and ephemeris parameters are accurate, the l cp from
equation (1.6) and l ob from equation (1.1) should be equal. Since the measurement l sc is
accurate compared to the attitude and ephemeris information, any systematic difference
between l cp and l ob can be attributed to the attitude and orbit errors. Thus we can use the
difference between l cp and l ob as the observable
dl = l cp  l ob = Ti(ref , vef )
(R cp − ref )
R cp − ref
 T(ξr, ξp, ξy) l sc
(1.7)
The task of precision modeling is then to calculate the correction to those satellite
ephemeris and attitude parameters (i.e., ref , vef and ξ's) so that the residuals of dl after
correction are minimized for all selected GCP's. The orbit correction is modeled as a
linear function of time for each component in the OB system. Referred to as the short arc
method, this purely geometric method shifts and rotates the short arc of orbit defined by
the original ephemeris points to fit the GCP measurements.
2. Linearizing the Observation Equations
In equation (1.7), the calculation of l ob can also be carried out through
l ob = T i(ref , vef )(Ref  ref ) /  Ref  ref 
(2.1)
if Ref is more conveniently accessible. Since equation (2.1) is simply the inverse of
equation (1.4) and equation (1.3), the l ob calculated from equation (2.1) is the same as the
one in equation (1.1) except for the possible inclusion of numerical errors. However, it
should be mentioned that the true relationship between l ob and the parameters is always
equation (1.1). Equation (2.1) should not be confused with this because Ref in equation
(2.1) is not an independent variable but a function of equation (1.4). So, in observation
(1.7) information about the attitude parameters is contained in l ob, and the information
about orbit parameters comes from l cp .
Since the measurement of l sc is 2dimensional in nature, only 2dimensional information is
contained in equation (1.7) though there are 3 components involved. If a look vector
(either l cp or l ob) has the three components in the OB system
l = {xl, yl, zl}
(2.2)
the real information in these three components can be summarized in two variables like
the original look angle measurements. We chose the following two variables
Version 3.2
99
7/9/98
IAS Geometric ATBD Version 3
δ = atan (yl / zl )
(2.3)
ψ = atan (xl / zl )
(2.4)
So that the three components of equation (1.7) can be reduced to the two equations:
α = δ cp − δob
(2.5)
β = ψ cp − ψob
(2.6)
Note that in equation (2.3) and (2.4) the components of xl, yl, zl can be that of lineofsight vector instead of unit look vector, so that δcp and ψcp are explicit functions of orbit
position. In that case zl is approximately the height of the satellite.
If we define
true value = approximate value + correction
and differentiate equations (2.3), (2.4) with respect to the orbit position (for δcp and ψcp ),
differentiate equation (1.1) with respect to the satellite attitude (for δob and ψob) at their
corresponding approximate values, then equations (2.5) and (2.6) can be linearized as the
function of correction parameters
α ≅ (cos2δcp / h) dy − (cosδcp sinδcp / h) dz + dξr
β ≅ (1.0 / h) dx − dξp + tanδcp dξy
(2.7)
(2.8)
where dx, dy and dz are the correction to satellite position vector rob in the OB system,
and dξ's are the corrections to the satellite attitude angle ξ's. Other quantities are
functions evaluated at the approximate values of ref , vef and ξ's.
The linearization above is done by directly differentiating equation (2.3) and (2.4), with
transformation Ti regarded unaffected by the error in ref and vef . This, however, ignores
the curvature of the satellite orbit and the Earth, resulting in about 10% of error in the
coefficients of dx, dy and dz. A more accurate way to evaluate these coefficients is to
examine the sensitivity terms dψcp /dx, dδcp /dy, and dδcp /dz through the geometry of the
look vector (see Figure 3.1.524).
Version 3.2
100
7/9/98
IAS Geometric ATBD Version 3
Figure 3.1.5 24 Look Vector Geometry
Denote
R – the Radius of the Earth;
r – the radius of the satellite position;
h – the altitude of the satellite;
d – the magnitude of the look vector (from satellite to target);
δ – the alongscan angle of the look vector;
φ − the Earth centered angle between the satellite and the target;
γ – the zenith angle of the look vector at the target;
x,y,z – the coordinates of the satellite position in the OB system
then we have
R sin(δ + φ) = r sin δ
(2.9)
Differentiating the equation (holding R and r constant) yield
R cos(δ + φ)(dδ + dφ) = r cosδ dδ
(2.10)
Note that δ + φ = γ, and dφ = −dy / r, we have
α = dδ = (−b / (r d)) dy
(2.11)
Similarly for the crossscan direction, we have
Version 3.2
101
7/9/98
IAS Geometric ATBD Version 3
β = dψ = ((r  d cosδ) / (r d cosδ)) dx
(2.12)
For the effect of altitude error, differentiate equation (2.9) with respect to δ and r (holding
φ constant now) and noting dr = dz, we have
α = dδ = (sin δ / d) dz
(2.13)
Note that the dx, dy and dz in equations (2.11) through (2.13) are error terms, which are
opposite in sign to the correction terms. With this in mind, we can replace the correction
terms in equation (2.7) and (2.8) and rewrite the linearized observation equation as
α = (b / (r d)) dy − (sin δ / d) dz + dξr
(2.14)
β = ((r  d cosδ) / (rd cosδ)) dx  dξp + tan δ dξy
( 2.15)
where
b = R cosγ = sqrt(R2 − (r2sin2δ))
(2.16)
d = r cos δ  b
(2.17)
3. Weighted Least Square Solution of the Parameters
The correction parameters in equations (2.14) and (2.15) can be expanded to include the
correction to the change rates of the satellite attitude and position by defining
dx = dx0 + dxdotdt
and
dξr = dξr0 + dξrdotdt
(3.1)
If the major random error source in computing α and β is the error in the coordinates of
the Ground Control Point, then the covariance matrix for the observation equations (2.14)
and (2.15) should be the covariance matrix of Rcp in equation (1.6) mapped through
equations (1.6), (2.3) and (2.4).
Note that in the observation equations (2.14) and (2.15), α is only related to parameters
dy, dz and dξr, and β is only related to dx, dξp and dξy. The parameters are uncoupled in
the two observations. In the simplified case where observational error of α and β are
uncorrelated, the observation equations can be separated into two independent equations
and solved individually.
Define the parameter vectors as
Version 3.2
102
7/9/98
IAS Geometric ATBD Version 3
X1’ = {dξr0, dy 0, dz 0, dξrdot, dy dot, dz dot}
(3.2)
X2’ = {dx0, dξp0, dξy0, dxdot, dξpdot, dξydot}
(3.3)
where ' means transpose of a vector or matrix. Then, the two observation equations can
be written as
α = h1 X1 + εa
(3.4)
β = h2 X2 + εb
(3.5)
where
h1 = {1.0,b/(d r), sin δ/d, dt, b dt/(d r), sin δ dt/d}
(3.6)
h2 = {(r  d cosδ)/(r d cosδ), 1.0, tanδ, (rd cosδ) dt/(r d cosδ), dt, tanδ dt}
(3.7)
εa and εb are the random error of α and β respectively. With all GCPs included, the
observation equation can be written as
A = H1X1 + εA;
(3.8)
B = H2X2 +εB;
(3.9)
and the parameters can be solved by Weighted Least Square estimation as
X1 = (H1’WaH1)1 (H1’WaA)
(3.10)
X2 = (H2’WbH2)1 (H2’WbB)
(3.11)
where A and B are the observation vectors, composed of α and β for all the GCPs
respectively, H1 and H2 are corresponding coefficient matrix with h1 and h2 as rows
corresponding to each α and β, Wa and Wb are the diagonal weight matrix for A and B
respectively, composed of inverse of the variance of each individual εa and εb.
One problem in this solution is the nearly linear correlation between parameter dx and
dξ p in the observation equation (3.7). The alongtrack orbit error and the pitch angle error
have the very similar effect on β and the two parameters can not be well separated in the
solution without additional information. Including both parameters in the observation
equations results in a nearsingular normal equation and therefore unstable solution of the
Version 3.2
103
7/9/98
IAS Geometric ATBD Version 3
parameters. Similarly, high correlation exists between the crosstrack position and the roll
attitude errors in equation (3.6) and an illconditioned normal equation would result.
For the purpose of correcting the image, we do not have to distinguish between orbit
position correction and attitude correction parameters. Letting either the orbit or attitude
correction parameters absorb the existing errors will correct the image in a similar manner.
Therefore, we can chose to estimate either dx and dy or dξp and dξr. This can be done by
setting those coefficients in h1 and h2 corresponding to the unwanted parameters to zero.
In practice the allocation of correction offsets to ephemeris and attitude terms is
controlled by the selection of apriori weights for the correction parameters.
4. Separating Orbit Position and Attitude Errors
One of challenging tasks in the Landsat 7 Image Assessment System is to distinguish
satellite attitude error from the orbit positional error. The purpose of precision correction
estimation is not only to correct the image but also to extract information about the sensor
alignment which is reflected in the attitude correction parameters. In order to separate the
ephemeris error from the attitude error as much as possible, we should first use the most
precise ephemeris data available and correct systematic errors with available models.
Second we should use available a priori information in addition to the observation to cure
the illcondition of the normal equation in statistical estimation.
Let the observation equation be
Y = HX + ε;
Ε[ε] = 0, Cov[ε] = s 2C
(4.1)
and the a priori information of the parameters be
X = X + εx;
Ε[εx] = 0, Cov[εx] = q2Cx
(4.2)
then the normal equation for the Best Linear Unbiased Estimate (BLUE) [Reference 14] ,
X^ of X is
((l/s2)H’WH + (l/q2)Wx )X^ = (l/s2)H’WY + (l/q2)Wx X
(4.3)
where W and Wx are weight matrices
W = C1;
Wx = Cx 1
Version 3.2
(4.4)
104
7/9/98
IAS Geometric ATBD Version 3
The covariance matrix of X^ is
Cov[X^] = ((l/s2)H’WH + (l/q2)Wx )1
(4.5)
Usually, however, the Cov[ε] and Cov[εx] can not be exactly known. In the case of
ground control points, for example, the position error involves many factors like base
map error, and human marking error, etc. If there are unknown scale factors s2 and q2, we
can still obtain the Weighted Least Square (WLS) estimate from the normal equation
(H’WH + W x )X^ = H’WY + Wx X
(4.6)
In such case, the inverse of the normal matrix can not be taken directly as the Cov[X^].
Factor s2 and q2 should be estimated with appropriate variance component estimation
from the residual of the solution of equation (4.6). The weighted residual square
summation can be calculated as
V’WV = Y’WY  2 X^’M + X^’ NX^
(4.7)
Vx’Wx Vx = X’Wx X  2 X^’Wx X + X^’Wx X^
(4.8)
where
V = Y  HX^
(4.9)
Vx = X  X^
(4.10)
N = H’WH
(4.11)
M = H’WY
(4.12)
5. Estimating the Weight Factor for Observation and a priori Information
When the factors s2 and q2 are appropriately estimated, the weight matrix W and Wx
should be correspondingly corrected by factors 1/s2 and 1/q2 respectively, and equation
(4.6) should be resolved with the new weight matrices. In the new solution, information
from the observation and the a priori information are appropriately combined and the
(H' WH + Wx )1 is the Cov[X^].
One of the estimate of s2 and q2 is the Helmert type estimate. For the problem here, the
equation for the estimate can be derived following Helmert's variance component analysis
[Reference 15],
E s2 + D q2 = V’WV
Version 3.2
(5.1)
105
7/9/98
IAS Geometric ATBD Version 3
D s 2 + G q2 = Vx’Wx Vx
(5.2)
where
E = n  2 tr{Q N} + tr{Q N Q N}
(5.3)
G = m  2 tr{Q Wx } + tr{Q W x Q Wx }
(5.4)
D = tr{Q N Q Wx}
(5.5)
Q = (H’ W H + W x }
(5.6)
n = number of observations
m = number of parameters
Equation (5.1) and (5.2) do not guarantee positive solution of s2 and q2. In some cases,
especially for small s2 and q2, noise can drive the solution negative. Another type of
estimate, the iterative Maximum Likelihood Estimate (MLH) guarantees positive solution
[Reference 16], though the estimate s2 and q2 may not be statistically unbiased. The MLH
solution is obtained by iteratively solving equation (4.6) and
s2 = V’WV / n
(5.7)
q2 = Vx’Wx Vx / m
(5.8)
W = W / s2
(5.9)
Wx = Wx / q2
(5.10)
until s2 and q2 converge.
Note: The weight factor estimates procedure is very sensitive to the presence of blunders
in the input data and therefore should not be used until the outlier rejection procedure,
section 3.1.5.3.7.1, has been performed.
6. Propagating the Parameter Errors into Look Vector Error
The solution above provides the estimate of the corrections to the ephemeris and attitude
data as well as to their covariance matrix. The covariance information can be used as a
measure of precision for assessing the alignment errors of the sensor system. It can also
Version 3.2
106
7/9/98
IAS Geometric ATBD Version 3
be propagated to any pixel in the scene to evaluate the pixel location error after the
precision correction.
Given the sample time and alongscan look angle of a pixel, the coefficients h1 and h2 can
be calculated for α and β according to equation (3.6) and (3.7). The variance of α and β
are then calculated as
σα2 = h1 Cov[X^]h1’
(6.1)
σβ2 = h2 Cov[X^]h2’
(6.2)
respectively. These are the variance of the pixel location in sample and line directions due
to the uncertainty of the estimated precision correction parameters, they are in angles but
can be easily converted into IFOV according to the sensor system specifications.
3.1.5.3.7.1 Precsion Correction Outlier Detection Algorithm
Background
Outlier detection for the IAS precision correction solutions seeks to identify GCPs which
are likely to be in error due to miscorrelation. This is done by analyzing the GCP
residuals, taking into account the relative importance of the GCP as reflected in the
precision solution normal equation matrix. [References 36, 37]
Definitions:
A = matrix of coefficients (partial derivatives) relating parameters to observations
θ = parameter vector
X = observation vector
V = residual vector
C = observation covariance matrix
n = the number of observations
p = the number of parameters
A is n x p, θ is p x 1, X and V are n x 1, and C is n x n
Observation Equation:
Aθ = X  V
X = X true + E
Version 3.2
where E = error vector ~ G(0,C)
107
7/9/98
IAS Geometric ATBD Version 3
Aθtrue = X true
where θtrue is the “true” parameter vector
Aθ = X true + E  V
so V = E if θ = θtrue
Minimum Variance Parameter Estimate:
θ’ = [ATC1A]1ATC1X
Estimated Residual Error:
V’ = X  A[ATC1A]1ATC1X
Define Projection matrix P:
P = A[ATC1A]1ATC1
This matrix projects the observation vector into the parameter subspace (the column
space of A).
This projection is only orthogonal if C has the special structure described below.
Substituting:
V’ = X  PX = [I  P]X
[I  P] projects X into the parameter null space.
Looking at the Error Estimate V’:
V’ = [I  P]X = [I  P][Xtrue + E] = [I  P]Xtrue + [I  P]E
but [I  P]Xtrue = 0 since Xtrue lies entirely within the parameter subspace.
so V’ = [I  P]E = E  PE
Some comments about V’ and E:
1. For a given precision solution the elements of E are not random variables, they are
realizations of random variables.
2. V’ is an estimate of the actual (realized) error E which includes an estimation error
equal to PE.
Version 3.2
108
7/9/98
IAS Geometric ATBD Version 3
3. We cannot exactly recover E from [I  P]1V’ because [I  P] is singular (it is an n x
n matrix of rank np).
4. We can attempt to predict how accurate our estimate (V’) of E is likely to be by
looking at the estimation error R = PE.
5. Since we want the predicted accuracy to apply in general, we treat R as a random
vector which is a function of another random vector E.
Expected Value:
E[ R ] = E[ P E ] = P E[ E ] = P 0 = 0
T
Variance:
E[ R R ] = E[ P E ET PT ] = P E[ E ET ] PT = P C PT
Special Structure of Observation Covariance Matrix for IAS Precision Correction:
C = σ2I
since the observation errors are realizations of independent and identically distributed
zero mean Gaussian random variables with variance σ2.
Substituting into the equation for P yields:
P = A[(1/σ2)ATIA]1ATI(1/σ2) = Aσ2[ATA]1AT(1/σ2) = A[ATA]1AT
And the equation for the variance of R:
E[ R RT ] = σ2 P I PT = σ2 P
noting that PT = P and P P = P
so R ~ G( 0, σ2 P )
For a particular component of R r i:
E[ ri ] = 0
E[ ri2 ] = σ2 p ii
where p ii is the ith diagonal component of P
Looking at the equation for P we see that:
p ii = A iT [ATA]1 Ai
where AiT is the ith row of A
Considering a particular component of the Residual Error Vector V’:
vi = ei  ri
where ei is the corresponding component of the observation error vector
so vi is an unbiased estimate of ei with variance σ2 p ii
Version 3.2
109
7/9/98
IAS Geometric ATBD Version 3
Outlier Detection
If we knew what e i was we could test it against a probability threshold derived from its
standard deviation, σ, to determine if it is likely to be an outlier. Instead of ei we have vi
which includes the additional error term ri. Including the additional estimation error in the
threshold computation leads to:
σv2 = σ2 + σ2 p ii
where σ2 is the term due to the actual error variance and σ2 p ii is the term due to the
estimation error variance.
This may seem like cheating since ei and ri are not independent for a given realization. In
fact:
E[ vi2 ] = E[ (ei  ri)2 ] = E[ ei2  2eiri + ri2 ]
and ri = Σj p ij ej
2
2
E[ vi ] = σ (1  pii)
It is tempting to use vi / (1  p ii)1/2 for ei in the outlier test (or, equivalently, to test vi
against a threshold based on σ2 (1  pii)) but this becomes dangerous as pii approaches 1.
The factor pii can be interpreted as a measure of the uniqueness of or as the information
content of the i th observation. As pii approaches 1 the ith observation lies almost entirely
within the parameter subspace which implies that it is providing information to the
solution which the other observations do not. Note that such “influential” observations
can be identified from the structure of the coefficient matrix, A, without reference to the
observation residuals. Attempting to use 1/(1  pii)1/2 to rescale the residual vi to better
approximate ei will, in a sense, punish this observation for being important. Instead, we
view p ii as a measure of how poor an estimate of the actual error, ei, the residual, vi, is and
ignore the fact that vi will tend to be an underestimate of ei. We therefore use σv2 = σ2 (1 +
p ii) as shown above) to construct the outlier detection threshold.
One remaining problem is that we do not know exactly what σ2 is and must estimate it
from the observation residuals. This is done by scaling the apriori observation variance
using the variance of unit weight computed in the precision solution. The fact that we are
using an estimated variance to establish our outlier detection threshold modifies the
algorithm in two ways: 1) we compensate for the fact that removing a point as an outlier
will alter the computation of the variance of unit weight by removing one residual and
reducing the number of degrees of freedom; and 2) we base the detection threshold
computation on Student’s tdistribution rather than the Gaussian distribution.
The variance of unit weight is computed as:
var0 = VTC1V/(n  p) = VTV/σ02(n  p) = Σvj2/σ02(n  p)
where: n = number of observations,
p = number of parameters, and
Version 3.2
110
7/9/98
IAS Geometric ATBD Version 3
σ02 is the apriori variance.
The estimated variance is:
var = var0 σ02 = Σvj2/(n  p)
Removing the kth observation makes this:
vark = (Σvj2  vk2)/(n  1  p) = (n  p)/(n  p  1) * (Σvj2  vk2)/(n  p)
vark = (n  p)/(n  p  1) * var  vk2/(n  1  p)
To normalize the kth residual we divide it by the estimated standard deviation σ’ =
(var)1/2:
wk = vk / σ’
We can rescale this normalized residual to reflect the removal of this observation from the
variance estimate without having to actually compute a new variance:
wk’ = vk / σk’ = w k σ’/σk’ = w k (var/vark)1/2
var/vark = 1 / [(n  p)/(n  p  1)  vk2/var (n  p  1)] = (n  p  1)/(n  p  vk2/var)
var/vark = (n  p  1)/(n  p  wk2)
noting that vk2/var = wk2
wk’ = w k [(n  p  1)/(n  p  wk2)]1/2
Finally, we include the (1 + pkk) factor discussed above and our normalized and
reweighted residual becomes:
wk’ = w k [(n  p  1)/(1 + pkk)(n  p  wk2)]1/2
where: wk = vk / σ’
This normalized and reweighted residual is compared against a probability threshold
computed using Student’s tdistribution with (n  p) degrees of freedom.
3.1.5.3.7.2 Inverse Mapping
For the precision correction process the information needed from the systematic ground
control point includes the satellite position and velocity vectors in earth centered inertial
coordinates along with the latitude and longitude associated with the point. From the
reference ground control point the precision process needs the latitude and longitude along
with the height of the control point above the WGS84 ellipsoid. The steps involved in the
process include:
Version 3.2
111
7/9/98
IAS Geometric ATBD Version 3
1. Use Geometric Correction Grid to find an input line and sample for an associated
output line and sample of the systematic Ground Control Point (GCP) (see the
Resampling algorithm description).
2. Determine the scan time associated with the systematic GCP. Since the model
truncates to determine the scan number from the input line, the input line has a value
of 0.5 added to it so that it can be rounded to the nearest integer line.
3. Use an interpolating polynomial for satellite position and velocity to get the
position and velocity vectors for the time associated with the input line and sample of
systematic GCP.
4. Call the Projection Transformation Package to convert the projection coordinate of
the systematic image GCP to latitude and longitude.
5. Call the Projection Transformation Package to convert the projection coordinate of
the base image GCP to latitude and longitude.
3.1.5.3.7.3
Precision Ephemeris Update
The purpose of this algorithm is to update the systematic orbit model with the
information from the precision correction solution parameters, so that a precision image
can be generated.
Given:
tref
xob
y ob
z ob
vxob
vy ob
vz ob
reference time for the precision correction;
orbit position correction in the Orbit System x (alongtrack) direction
orbit position correction in the Orbit System y (crosstrack) direction
orbit position correction in the Orbit System z (altitude) direction
orbit velocity correction in the Orbit System x (alongtrack) direction
orbit velocity correction in the Orbit System y (crosstrack) direction
orbit velocity correction in the Orbit System z (altitude) direction
For each ephemeris data point, the algorithm takes the following ephemeris information:
t eph
x
y
z
vx
vy
Version 3.2
time for the ephemeris data point
orbit position in the Earth Fixed system x direction
orbit position in the Earth Fixed system y direction
orbit position in the Earth Fixed system z direction
orbit velocity in the Earth Fixed system x direction
orbit velocity in the Earth Fixed system y direction
112
7/9/98
IAS Geometric ATBD Version 3
vz
orbit velocity in the Earth Fixed system z direction
and performs the following steps:
1. Calculate the orbit position correction in Orbit System at the time
dxob = xob + vxob * (teph  tref)
dyob = yob + vyob * (teph  tref)
dzob = zob + vzob * (teph  tref)
2. Calculate the transformation matrix from the Orbit System to the Earth Fixed
system at the time
2.1. First, compose the transformation matrix, Teo, from the Earth Fixed system to
the Orbit System:
Wz, unit vector in negative radial direction, pointing to the Earth center:
(denoting X = {x, y, z} and V = {vx, vy, vz})
Wz =  X / X
Wy, unit vector in crosstrack direction ( angular momentum);
Wy = Wz cross V / Wz cross V
Wx, unit vector in alongtrack direction;
Wx = Wy cross Wz
W x (1) W x (2) W x (3)
Teo = W y (1) W y (2) W y (3)
W z (1) W z (2) W z (3)
2.2 Transpose the Teo matrix to get the transformation from the Orbit System to
the Earth Fixed system:
Toe
Version 3.2
W x (1) W y (1) W z ( 1)
= W x (2) W y (2) W z (2)
W x (3) W y (3) W z (3)
113
7/9/98
IAS Geometric ATBD Version 3
3. Calculating the orbit position correction in the Earth Fixed system:
dx W x (1) W y (1) W z (1) dx ob
dy = W (2) W (2) W (2) * dy
y
z
ob
x
dz W x ( 3) W y (3) W z ( 3) dz ob
4. Calculating the orbit velocity correction in the Earth Fixed system:
dvx W x (1) W y (1) W z (1) vxob
dvy = W (2) W (2) W (2) * vy
y
z
ob
x
dvz W x (3) W y ( 3) W z (3) vz ob
5. Update the orbit position and velocity in the Earth Fixed system:
x = x +
y = y +
z = z +
vx = vx
vy = vy
vz = vz
dx;
dy;
dz;
+ dvx;
+ dvy;
+ dvz;
This algorithm computes the orbit correction for each ephemeris data point in the Earth
Fixed system and adds to yield the new ephemeris data point information.
3.1.5.3.8 Terrain Correction
The terrain correction algorithm uses elevation data from a digital elevation model (DEM)
to generate offsets which will correct for the effect of terrain parallax. These corrections
are only required in the alongscan (sample) direction.
To avoid having to call the Projection Transformation Package for every point in the
output space, the DEM is processed to the same pixel size, projection, and frame as the
precision image and its elevations are relative to the WGS84 ellipsoid rather than mean sea
level or some other vertical datum. It is also expected that the precision corrections have
already been applied to the forward model. The process is as follows:
1. Read the DEM and get minimum and maximum elevations present in the scene.
2. Determine the location of nadir in terms of raw space input sample. The secant
method of root finding is used to determine the sample location at which the ycomponent of the line of sight vector (LOS) in orbital coordinates goes to zero for
Version 3.2
114
7/9/98
IAS Geometric ATBD Version 3
each scan. It stores these in an array. This segment also stores the latitude and
magnitude of the satellite position vector (in meters) for the center of the scene.
3. Using the minimum and maximum elevations found in Step 1, determine the
sample displacement due to terrain for every (offnadirsample, elevation) pair.
Offnadirsample should range from zero to 3600 (or some number which
guarantees that all possible input samples will be covered). Elevation should be in
meters and range from the minimum to the maximum in the scene. The algorithm
steps elevation in increments of 1 meter.
4. For every line, sample pair in the output image space, the terrain correction
algorithm needs the corresponding line, sample in input space. (See the inverse
model, section 3.1.5.3.6).
5. Given the current output line and sample, the elevation is read from the DEM
buffer from Step 1.
6. Given the corresponding input sample, it's distance from nadir in pixels is
determined from the calculation of nadir pixel in Step 2.
off_nadir_pixel = in_sample  nadir_sample
7. From the (off_nadir_pixel, elevation) pair determined in Steps 6 & 7, the value of
the terrain offset is looked up from the table generated by Step 3.
8. The next step in terrain correction, (applied in the Resampler) is to adjust the
input sample locations by the terrain offset and resample.
Note that the algorithm assumes a DEM is available in the same projection, datum and
framing as the output space as defined by the precision corrected grid; and that it has been
resampled to the same pixel size. It also assumes that for every line, sample in output
space, the corresponding line, sample can be found in input space. The offsets due to
terrain are applied in the resampler.
This algorithm was investigated using the 1°x 1° (3arcsec) DTED1 (Digital Terrain
Elevation Data Level1) available via anonymous FTP from the EROS Data Center.
Slightly better accuracy could be achieved using higher spatial resolution 7.5 x 7.5arcmin
(30m) USGS Digital Elevation Model (DEM) data; vertical accuracy is better than 15m
RMSE.
Terrain Correction Algorithm
Version 3.2
115
7/9/98
IAS Geometric ATBD Version 3
The terrain correction algorithm will generate a table of terrain induced sample offsets,
corresponding to (offnadirsample, elevation) pairs, as generated in Step 3 above. The
method used for determining these terrain corrections comes from Reference 22. The
following summarizes that document. A schematic of the geometry is shown in Figure
3.1.525.
Given distance in pixels of a sample in input space from the nadir pixel location, the
ground range, S, is given by:
S = pixel_size * off_nadir_pixel
Where pixel_size is the nominal pixel size in meters for the given band.
This is associated with an Earth centered angle, s:
s = S/Re
Where Re is the radius of the Earth reference ellipsoid at the scene center latitude.
The magnitude and direction of the line of sight vector relative to nadir can then be found
as:
LOS = sqrt(Re2 + (Re + Alt) 2  2*Re(Re + Alt)*cos(s))
d = asin(Re*sin(s)/LOS)
Where Alt is the altitude of the satellite above the ellipsoid.
Looking at Figure 3.1.525, elevation, h, at location S, will cause the information at the
top of the hypothetical mountain to be mapped to the location S + dS on the ellipsoid
used by the forward sensor model. What the terrain correction will do is determine the
offset, dS, so that when the time comes to resample, we can grab the information
projected on the ellipsoid at location S + dS and move it back to the mountain top at
location S. This is done by first determining the angular difference between the LOS to a
spot on the ellipsoid below the terrain feature and the LOS' to the projected location of
that feature on the ellipsoid.
dd = atan{(Re+Alt) * sin(d) * (1(Re+h)/Re) / [(Re+h) * sqrt(1(Re+Alt)2 *
sin2(d)/Re2)(Re+Alt)*cos(d)]}
z'' = asin[(Re+Alt) * sin(d+dd)/Re]
ds = z''  s  (d+dd)
dS = Re * ds
Version 3.2
116
7/9/98
IAS Geometric ATBD Version 3
terrain_offset = dS / pixel_size
Figure 3.1.5 25 Terrain Correction Geometry
Version 3.2
117
7/9/98
IAS Geometric ATBD Version 3
Then for resampling, rather than getting data from input(sample), the data for the location
in output space should be taken from input (sample + terrain_offset * (sign of
off_nadir_pixel)).
The assumption of a constant satellite altitude above the geoid for one Landsat scene was
tested with an absolute worst case scenario. A terrain feature with an elevation of 8,850m
(the height of Mt. Everest) was placed near the equator at the end of a scan line. The
orbit was modeled using the highest eccentricity allowed within specifications and in the
vicinity of the greatest rate of change of orbit radius. In this case, the maximum error
introduced was only ~0.03 panchromatic pixels. This is approximately the same size as
the resampling increment, 0.03125 pixels. In a test case using Landsat 5 data from an
actual test site, (path 28, row 30), the error from the altitude assumption for an elevation
of 500m at the end of scan was ~0.001 pixels. For these reasons, the assumption of a
constant satellite altitude over a scene was considered valid.
3.1.5.4 Geometric Calibration Algorithms
3.1.5.4.1 Sensor Alignment Calibration
The purpose of this algorithm is to use a sequence of precision correction solutions,
including both correction parameter estimates and estimated covariance information, to
estimate the underlying ETM+ instrument to Landsat 7 navigation reference base
alignment. A Kalman filter implementation is used to isolate the systematic alignment
trend from the attitude and ephemeris knowledge error which varies from scene to scene.
This algorithm is one step of the Kalman filter [Reference 23] for a linear dynamic model
(constant rate model). It takes one precision correction solution file (generated from one
scene) and updates the estimate of the alignment state parameters. It appends the newest
state parameter estimate to the orbit parameter sequence file, attitude parameter sequence
file, and also appends an innovation sequence file. This algorithm does not make a
decision on whether the alignment is anomalous or whether it should be corrected. This
algorithm provides the information (the sequence files) for an analyst to make that
decision.
This algorithm would be applied each time a calibration scene is processed through the
precision correction software, to update the alignment estimate. Operationally, the
precision correction processing will include the use of definitive ephemeris data, obtained
from the GSFC Flight Dynamics Facility (FDF), to help separate ephemeris position
errors from angular alignment and attitude errors. Error propagation analyses and
simulations run with the algorithm prototype software suggest that it will require 4 or 5
precision corrected calibration scenes to drive the alignment knowledge from its pre
Version 3.2
118
7/9/98
IAS Geometric ATBD Version 3
launch accuracy of approximately 240 arc seconds per axis (1 sigma) to the IAS calibrated
specification value of 24 arc seconds per axis.
Algorithm Description
Defining the Filter:
The state vector of this filter is a 12x1 vector:
X=
 roll 
 pitch 
 yaw 
 x

 y

 z

 r_dot 
 p_dot 
yaw_dot
 vx 
 vy 
 vz 
roll attitude angle, in microradians (urad);
pitch attitude angle, in urad;
yaw attitude angle, in urad;
alongtrack orbit position error, in meters (m);
crosstrack orbit position error, in m;
radial orbit position error, in m;
roll angle rate, in urad/day;
pitch angle rate, in urad/day;
yaw angle rate, in urad/day;
velocity in x, in m/day;
velocity in y, in m/day;
velocity in z, in m/day;
See the precision correction algorithm description for the definition of above parameters,
but note the units for the rate and velocity terms are different here. In fact, the rate and
velocity state vector components represent long term linear trends in the bias data rather
than the short term rate of change errors estimated by the precision correction process.
The covariance matrix for X is a 12x12 matrix P.
The observation vector for this filter is a 6x1 vector:
Y=
 roll
 pitch
 yaw
 x
 y
 z






roll attitude angle, in urad;
pitch attitude angle, in urad;
yaw attitude angle, in urad;
alongtrack orbit position error, in m;
crosstrack orbit position error, in m;
radial orbit position error, in m;
The covariance matrix for Y is a 6x6 matrix R. The observation is from the precision
correction solution, but with the velocity and attitude rate parameters discarded. The
velocity and rate estimated from precision correction are short term variations in a scene,
what is to be monitored in the alignment calibration is the long term drift.
Version 3.2
119
7/9/98
IAS Geometric ATBD Version 3
The precision correction estimates for the attitude bias terms: roll, pitch, and yaw; are
actually estimates of the combination of long term alignment error (what we are looking
for) plus short term alignment and attitude errors. These short term errors appear to be
constant, or nearly constant, biases over time intervals of minutes to hours, but will
decorrelate over the time intervals between calibration scenes (days). These short term
variations thus appear to be random noise superimposed on the alignment bias
observations taken from each calibration scene. This will manifest itself in the calibration
scene precision correction results as sequences of attitude bias estimates whose
dispersion is greater than would be expected based on the aposteriori covariance estimates
computed in the precision correction process. Any systematic linear trend present in the
attitude error will be indistinguishable from alignment and will be incorporated into the
alignment bias estimated by the Kalman filter.
The random attitude and dynamic alignment errors are, in effect, an additional
measurement error vector, added to the estimation error from the precision correction
solution, in the Kalman filter formulation. The observation covariance matrix, R, is
therefor the sum of the aposteriori covariance matrix from the precision correction
solution and the expected covariance of the “random” attitude and dynamic alignment
errors. According to the Landsat 7 System Specification [Reference 5], the expected
variance of the additional random attitude and dynamic alignment errors is less than or
equal to the sum of the squares of 45 arc seconds of attitude knowledge error (Reference
5, section 3.7.1.3.9), 8 arc seconds of dynamic ETM+ to Navigation Base Reference
alignment error (Reference 5, section 3.7.1.3.10), and 15 arc seconds of line of sight to
ETM+ reference dynamic alignment error (Reference 5, section 3.7.8.1.16.7). This yields
a net random component of 48 arc seconds. Since these errors are assumed to be
uncorrelated with each other, they are added to the diagonal elements of the R matrix
corresponding to the alignment biases.
Although we do not expect to detect any systematic offset in the position bias terms: x,
y, and z; they are included because of their high correlation with the attitude biases. This
is reflected in the observation covariance matrix where significant offdiagonal terms will
exist for xpitch and yroll. Any particular precision correction solution will resolve the
correlation between the parameters by allocating the along track and across track errors
between the ephemeris and attitude parameters based on their a priori weights. Thus,
some of the systematic alignment bias could end up allocated to the ephemeris correction
terms. Over multiple precision correction solutions the net ephemeris bias should be very
close to zero, so we establish small a priori state variances for the ephemeris bias terms
and allow the Kalman filter to reallocate systematic trends in the ephemeris corrections
back to the alignment parameters. In practice, since accurate postpass ephemeris will be
available for the calibration test scenes, most of the correction will be allocated to the
attitude terms in the precision correction solutions anyway.
Version 3.2
120
7/9/98
IAS Geometric ATBD Version 3
The dynamic model used for this filter is the constant rate (velocity) model:
X k = F( t k , t k1 )X k1
where the state transition matrix is:
I dt I
F(t k , t k1 ) =
0 I
I is a 6x6 identity matrix and dt = t k  t
t k.
k1
is the time difference between epoch t k1 and
The observation equation is
Y = H*X + v;
where H = [ I 0 ] and v is the net measurement error with covariance R.
Process noise can be added, with constant covariance matrix Q. Q is diagonal with 12
components. The process noise will initially be assumed to be zero. After several
precision correction solution sets are available for analysis, the innovation sequence and a
posteriori state covariance will be assessed to decide whether process noise should be
added to better model the dynamic behavior observed in the data.
Procedure:
For each run this algorithm performs the steps described in the following subparagraphs.
1. Setup option
Read in the input option file to define the input and output file names, options and
process noise variance.
2. Initialization
Read in the previous epoch time t p, state vector Xp, and covariance matrix Pp. For a regular
run this is read from the output state file of the previous run. For the first run, this is read
from an initial observation file, which contains the a priori knowledge of the state vector
Version 3.2
121
7/9/98
IAS Geometric ATBD Version 3
(alignment terms are zero plus or minus 240 arc seconds, ephemeris terms are zero plus or
minus one meter).
3. Read and preprocess an observation as follows:
3.1 Extract the observation data from the precision correction solution
Read in the current epoch time tc , precision correction observation Z, and covariance
matrix R. This is read from a solution file output by the precision correction algorithm.
The precision correction solution file contains both bias and rate correction terms whereas
the Kalman filter observation only uses the bias components. The terms corresponding to
roll, pitch, yaw, x, y, and z bias are extracted from the precision correction solution vector
and covariance matrix to assemble Y and R, respectively. The additional observation noise
due to scene to scene attitude and dynamic alignment error is accounted for by adding the
expected error variance (48 arc seconds)2 to the diagonal elements of R corresponding to
the roll, pitch, and yaw bias correction parameters.
3.2 Convert the attitude observation to an absolute pointing error observation
The precision correction solution contains updates to the alignment knowledge from the
Calibration Parameter File (CPF) used to process the test scene. We want the Kalman
filter to model the absolute alignment angles so that observations processed with different
versions of the CPF can be combined. To make this possible we must construct the
composite CPF alignment/precision correction angles for input to the Kalman filter.
The precision correction roll (r), pitch (p), and yaw (y) update values can be composed
into an alignment correction rotation matrix as follows:
∆A =
cos(p) cos(y) sin(r) sin(p) cos(y) + cos(r) sin(y)  cos(r) sin(p) cos(y) + sin(r) sin(y)
 cos(p) sin(y)  sin(r) sin(p) sin(y) + cos(r) cos(y) cos(r) sin(p) sin(y) + sin(r) cos(y)
sin(p)
 sin(r) cos(p)
cos(r) cos(p)
The Calibration Parameter File (CPF) alignment matrix is defined as the rotation from the
Navigation Base Reference (NBR) coordinate system to the ETM+ sensor coordinate
system. To convert a look vector constructed in ETM+ sensor coordinates (XETM) to the
orbital coordinate system (Xorb), where the precision correction solution is performed, the
following rotations must be performed:
Xorb = TTatt AT XETM
where: Tatt is the attitude rotation matrix, and A is the alignment matrix.
Version 3.2
122
7/9/98
IAS Geometric ATBD Version 3
Note that the transpose of A is used assuming that A is an orthogonal rotation matrix so
that its transpose is the same as its inverse.
The precision correction solution computes updates to the composite Tatt AT rotation to
best fit the observed ground control. This update takes the following form:
Xorb = ∆A TTatt AT XETM
where: ∆A is the alignment correction matrix defined above.
For small rotations, the ∆A and Tatt commute to a good approximation. Using the
maximum offnadir attitude pointing of 1080 arc seconds (900 arc seconds maximum offnadir commanded pointing per the Landsat 7 System Specification section 3.7.1.3.8.2
plus 180 arc seconds of pointing error per section 3.7.1.3.8) and the maximum expected
alignment error of 720 arc seconds (ibid. section 3.7.1.3.11), the error in commuting these
rotations is less than 0.01 arc seconds. This allows the equation above to be rewritten as:
Xorb = TTatt ∆A AT XETM
The composite rotation A’ = ∆A AT captures the net measured alignment offset. We
extract roll, pitch, and yaw angles from this composite matrix using the following
equations, where a’ij indicates the element from the ith row and jth column of A’.
roll’ = tan1(a’32/a’33)
pitch’ = sin1(a’31)
yaw’ = tan1(a’21/a’11)
These composite angles are inserted into the observation vector, Y, replacing the precision
correction output values.
4. Propagate Xp and Pp to the current epoch time tc .
Calculate:
dt = tc  t p;
Xc = F(t c ,t p) * Xp;
Pc = F(t c ,t p) * Pp * FT(t c ,t p);
Because of the special structure of F, the computation can be simplified if we divide the
X vector into 2 blocks:
X =
Version 3.2
 X1 


123
7/9/98
IAS Geometric ATBD Version 3
 X2 
and P into 4 submatrices:
P =
 P11  P12 
  +  
 P21  P22 
Then:
Xc1 = Xp1 + dt*Xp2
Xc2 = Xp2
Pc11 = Pp11 + dt*Pp12 + dt*Pp21 + dt*dt*Pp22
Pc12 = Pp12 + dt*Pp22
Pc21 = Pp21 + dt*Pp22
Pc22 = Pp22
And add the process noise if that option is selected:
Pc = Pc + Q
5. Calculate the Kalman gain matrix K:
K = Pc *HT*[H*Pc *HT + R]1
Because of the structure of H:
H*Pc *HT = Pc11
and:
Pc *HT =
 Pc11 


 Pc21 
which gives the Kalman gain matrix the following structure:
K=
 K1 


 K2 
 Pc11 
= 
 * [Pc11 + R]1
 Pc21 
6. Update the current state with the filtered innovation to create a new estimate and
update the new filtered state covariance:
Xn = Xc + K*(Y  H*Xc ) = Xc + K*(Y  Xc1 )
Pn = Pc  K*H*Pc
Version 3.2
124
7/9/98
IAS Geometric ATBD Version 3
which reduces to:
Pn =
Pc 
 K1*Pc11

 K2*Pc11
K1*Pc12 

K2*Pc12 
7. Append the orbit elements in Xn to the orbit sequence file.
8. Append the attitude elements in Xn to the attitude sequence file.
9. Append the innovation Y  H*Xc to the innovation sequence file.
10. Write out Xn and Pn to the new state file.
11. Update the Calibration Parameter File alignment matrix
Periodically, the analyst will decide that the alignment knowledge has improved
sufficiently to justify modifying the Calibration Parameter File. This would typically
occur after the first 5 or 6 calibration scenes have been processed postlaunch and then at
least once every quarter after the initial on orbit calibration. The update would be
performed using the most current roll, pitch, and yaw correction estimates in the attitude
sequence file.
The current filter estimates of roll (r), pitch (p), and yaw (y) are composed into a net
alignment correction rotation matrix, as was done for the differential corrections above, as
follows:
cos(p) cos(y) sin(r) sin(p) cos(y) + cos(r) sin(y)  cos(r) sin(p) cos(y) + sin(r) sin(y)
Anet =  cos(p) sin(y)  sin(r) sin(p) sin(y) + cos(r) cos(y) cos(r) sin(p) sin(y) + sin(r) cos(y)
sin(p)
 sin(r) cos(p)
cos(r) cos(p)
The net alignment rotation Anet is the transpose of the new alignment matrix, so the
updated alignment matrix for the CPF is computed as follows:
Anew = AnetT
Once this update is incorporated into the CPF, subsequent calibration scenes will be
processed with the improved alignment knowledge. Since this alignment knowledge is
combined with the precision correction observation, as described in section 2.2.3.2 above,
for input to the Kalman filter, the filter need not be reset after CPF updates. The state
covariance would also not be reset since it accurately reflects the confidence of the current
alignment knowledge.
Version 3.2
125
7/9/98
IAS Geometric ATBD Version 3
After a sufficient number of precision correction results have been collected, an analyst
may come up with some idea about the process noise by looking at the innovation
sequence and the output filter covariance. The process noise sigma can be fed back to the
filter through the input option file. Without the additional process noise, the filter may
begin to diverge after a certain time (the filtered state covariance will begin to grow).
Tuning the process noise is a difficult part of Kalman filtering practice, and some
experimentation will be needed to optimize the algorithm.
3.1.5.4.2 Focal Plane Calibration
Focal Plane  Band Placement Calibration Algorithm Description
The objective of this algorithm is to estimate updates to the focal plane locations of the
eight ETM+ band centers to improve band to band registration as specified in Reference
3.
Input to this algorithm includes the band to band comparison point correlation results
from the band to band characterization algorithm, with outliers detected and removed as
described in section 3.1.5.5.1. These are used in conjunction with the resampling grid file
(used to produce the level 1Gs test image) as input to the band to band calibration
process.
Algorithm:
The band location parameter vector is of the form:
∆b1x
∆b
1y
X = ...
∆b
8x
∆b 8y
where:
∆bix = band offset in the focal plane x (along scan) direction for band i, in
microradians
∆biy = band offset in the focal plane y (across scan) direction for band i, in
microradians
The band to band characterization test point correlation output data are assembled, point
by point, to construct a band to band calibration measurement vector:
Version 3.2
126
7/9/98
IAS Geometric ATBD Version 3
m 21P1
m 21L1
...
m 21PN
m 21LN
m 31P1
Y=
...
m 81LN
m 32P1
...
m 87LN
where: m jiPk = observed offset for point k, in output space pixels, between band j
and band i (reference band) = pixel jk  pixel ik
m jiLk = observed offset for point k, in output space lines, between band j
and band i (reference band) = line jk  line ik
The parameter vector and measurement vector are related by the following observation
equation:
dp
dp
dp
dp
−
... 0 ...
... 0
0 ... − dx
dy ik
dx jk dy jk
ik
X
A jik X =
dl
dl
dl
dl
0 ... −
−
... 0 ...
... 0
dx ik
dy ik
dx jk dy jk
m jiPk
=
+ V jik = Yjik + V jik
m jiLk
where: Vjik is the 2x1 vector of measurement residuals for point k, assumed to be
zero mean Gaussian errors with covariance matrix Cjik. Typically, the Cjik
would be assumed to be the same for all k, so Cjik = Cji.
σ 2p
C ji =
0
0
σ p ≈ σ l ≈ 0.25
σ 2l
The partial derivatives in the observation coefficient matrix Ajik, are computed
numerically using the resampling grids for bands j and i, from the input grid file, as
follows:
Version 3.2
127
7/9/98
IAS Geometric ATBD Version 3
For band i partials:
If band i was resolution adjusted for correlation with a lower resolution band,
adjust the measured pixel, line to match the grid:
If i=8 (i.e., this band is panchromatic)
p ik = 2*pik  1
lik = 2*lik  1
If j=6 (i.e., the other band is thermal)
p ik = 2*pik  1
lik = 2*lik  1
Use GRIDi to inverse map the test point pixel line coordinates: p ik, l ik to input
image space pixel, line coordinates s ik, d ik
Add a small epsilon value, e s, to s
ik
Use GRIDi to map ( sik + es , dik ) to output space ( pix , lix )
Add a small epsilon value, e d, to d
ik
Use GRIDi to map ( sik , dik + ed ) to output space ( piy , liy )
Compute the partials numerically, dividing by the band IFOV in microradians to
rescale the partial derivatives from output pixels per input pixel to output pixels
per microradian:
dp/dxik = (pix  pi)/es/IFOVi
dl/dxik = (lix  li)/es/IFOVi
dp/dyik = (piy  pi)/ed/IFOVi
dl/dyik = (liy  li)/ed/IFOVi
If band i was resolution adjusted for correlation with a lower resolution band,
adjust the partial derivative to match the actual measurements:
If i=8 (i.e., this band is panchromatic)
dp/dxik = (dp/dxik)/2
dl/dxik = (dl/dxik)/2
dp/dyik = (dp/dyik)/2
dl/dyik = (dl/dy ik)/2
If j=6 (i.e., the other band is thermal)
dp/dxik = (dp/dxik)/2
Version 3.2
128
7/9/98
IAS Geometric ATBD Version 3
dl/dxik = (dl/dxik)/2
dp/dyik = (dp/dyik)/2
dl/dyik = (dl/dy ik)/2
The same procedure is used with the grid file for band j, GRID j, to yield the band
j partials: dp/dxjk, dl/dxjk, dp/dyjk, and dl/dyjk. In this case, i and j are interchanged
in the procedure above.
The normal equations can be constructed directly from a sequence of observation
equations:
N = Σ AjikT Cji1 Ajik
summed over all k and all i,j combinations
T
1
L = Σ Ajik Cji Yjik
summed over all k and all i,j combinations
With only the band difference observations this system of equations will be indeterminate
since the solution is insensitive to the addition of constant offsets to the ∆b terms. To
stabilize the solution we add an additional observation requiring the band offset for band 8
(panchromatic) to be zero in each direction: Band 8 is constrained because it is used for
scan mirror calibration and sensor alignment calibration, making it the geometric
calibration reference band. This constraint is treated as another observation of the form:
0 0 0 0
A 00 =
0 0 0 0
0
Y00 =
C 00 =
0
0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 1
σ 20
0
0
σ 20
σ 0 ≈ 0.01 urad
After adding this observation to the normal equations, the minimum variance estimate for
X is:
Xest = N1 L
The measurement residuals for each input point can be calculated as:
Vjik = Ajik Xest  Yjik
These residuals would be assessed by band pair to determine whether there were
additional uncompensated biases in particular band combination measurements and to
evaluate the point matching accuracy of particular band combinations.
Output: The resulting parameter vector, Xest, contains the estimates of the 8 band offset
pairs (x and y directions) in microradians. These values would be used to update the Band
Offsets (expressed in microradians) in the Focal Plane Parameters parameter group in the
Calibration Parameter File, as defined in Reference 24.
Version 3.2
129
7/9/98
IAS Geometric ATBD Version 3
3.1.5.4.3 Scan Mirror Calibration
The purpose of the Scan Mirror Calibration segment of the Landsat 7 Image Assessment
System is to check and, if necessary, correct the along scan and across scan mirror
profiles.
As a preliminary approach to scan mirror calibration, a precisionterrain corrected
Landsat panchromatic band is matched to a reference image with correlation. Currently
the best candidate for the reference is the USGS digital orthophoto quadrangle (DOQ)
imagery. Using the precision mapping grid, correlation windows can be selected which lie
entirely within one scan thus allowing forward and reverse mirror scans to be separated.
For each successful correlation between a point in the reference image and a point in the
Landsat scene, the respective image coordinates (lines, samples) in output space are
mapped back through the inverse grid into input space (Resampling, section 3.1.5.3.6).
After being mapped to input space, the effects of terrain are removed from the two
samples to locate the (line, sample) input space coordinate associated with the measured
output space (line, sample) (Terrain Correction, section 3.1.5.3.8). The remaining
difference between the samples/lines is then taken to be the error in the scan mirror
profile. Using the forward model, the time in scan of each sample is then calculated from
the input sample (LOS Generation, section 3.1.5.3.4).
This correlation and inverse mapping procedure is used to generate a number of mirror
error data points for several forward and reverse scans, each with an associated time in
scan. Outlier rejection is performed by grouping the data points by scan direction and
scan angle. Points are removed that deviate from the mean..
Once the data has been filtered for outliers, the data from forward and reverse are fit with
fifth order Legendre polynomials as developed below. If the error in the mirror profile is
beyond some predetermined value, the existing mirror coefficients can be updated in the
calibration parameter file. The error polynomial coefficients are saved for trending.
Algorithm Development
The prelaunch mirror profiles are characterized as deviating from linear by a fifth order
polynomial
b(t) = b 0 + b1t + b 2t 2 + b3t 3 + b4t 4 + b5t 5
0 < t < Tnom
(1)
This can also be written as an orthogonal fifth order Legendre polynomial
p(x) = a0 + a1x + a2(3/2x2  1/2) + a3(5/2x3  3/2x) + a4(35/8x4  30/8x2 + 3/8)
Version 3.2
130
7/9/98
IAS Geometric ATBD Version 3
+ a5(63/8x5  70/8x3 + 15/8x)
1 < x < 1
(2)
where x = 2t/Tnom  1
In order to convert back and forth between the mirror coefficients, bi, and the Legendre
coefficients, ai, we equate the two polynomials above and gather terms of the same order.
This leaves us with the transformations
b0 = a0  a1 + a2  a3 + a4  a5
b1 = (2a1  6a2 + 12a3  20a4 + 30a5)/Tnom
b2 = (6a2  30a3 + 90a4  210a5)/Tnom2
b3 = (20a3  140a4 + 560a5)/Tnom3
b4 = (70a4  630a5)/Tnom4
b5 = (252a5)/Tnom5
a5 = (b5T nom5)/252
a4 = (b4T nom4)/70 + (b5T nom5)/28
a3 = (b3T nom3)/20 + (b4T nom4)/10 + (b5T nom5)5/36
a2 = (b2T nom2)/6 +(b3T nom3)/4 + (b4T nom4)2/7 + (b5T nom5)25/84
a1 = (b1T nom)/2 + (b2T nom2)/2 + (b3T nom3)9/20 + (b4T nom4)2/5 + (b5T nom5)5/14
a0 = b0 + (b1T nom)/2 + (b2T nom2)/3 + (b3T nom3)/4 + (b4T nom4)/5 + (b5T nom5)/6
(3)
(4)
As was shown in the section on scan mirror correction (Create Model, section 3.1.5.3.3),
the coefficients of the mirror scan profile are scaled by the total time of each scan which is
calculated from the first half and second half scan errors
bi’ = bi(T nom/Ts)i
bi’(T s)i = bi(T nom)i
(5)
(6)
From the transformations in Eqs. 4, the relationship in Eq. 6 leads directly to
aI` = ai
(7)
where aI’ are the coefficients of a polynomial in the form of Eq. 2 but where the auxiliary
variable x has been replaced by
x' = 2t/Ts 1
(8)
Therefore, using the actual scan time, Ts, instead of the nominal scan time, Tnom, to
compute the auxiliary variable, x, takes care of the scan time scaling without changing the
Legendre polynomial coefficients, a i .
Version 3.2
131
7/9/98
IAS Geometric ATBD Version 3
As discussed in section 3.1.5.3.3 (scan mirror correction), in addition to the fifth order
deviation from a linear profile, there is also an adjustment to the quadratic terms used to
fit the midscan error.
The linear portion of the mirror model is
L(t) = (amet + a sm(T s  t))/Ts
l(x) = (ame + asm)/2 + x(ameasm)/2
(9)
(10)
where asm and ame are the mirror angles at start and end of scan respectively. The
predicted mirror angle at the time of first half scan, tfh , is
L(tfh) = (amet fh + asmt sh)/Ts = l(ε)
(11)
t sh = T s  t fh
ε = 2t fh/Ts  1 = (t fht sh)/Ts
(12)
where
Thus as shown in the section on scan mirror correction, the midscan correction due to
deviation from linear (the difference between what the linear model predicts and what was
observed, i.e. zero) is
Af = (amet fh + asmt sh)/Ts = l(ε)
(13)
So the total midscan correction, including the mirror profile, is then
Df = A f  p(ε) = l(ε)  p(ε)
(14)
Again referring to the section on scan mirror correction, the coefficients for the quadratic
adjustment to the fifth order polynomial are
∆b0' = 0
∆b1' = (D fT s)/(tfht sh)
∆b2' = D f/(tfht sh)
(15)
Using the transformations in Eqs 4, these corrections can be written in terms of Legendre
coefficients as
∆a0 = ∆b0' + (∆b1'T s)/2 + (∆b2'T s2)/3 = (DfT s2)/(6tfht sh)
∆a1 = (∆b1'T s)/2 + (∆b2'T s2)/2 = 0
∆a2 = (Db 2'T s2)/6 = (DfT s2)/(6tfht sh)
Version 3.2
132
(16)
7/9/98
IAS Geometric ATBD Version 3
These can be rewritten as
∆a0 = (A f  p(ε))Ws = (l(ε) + p(ε))Ws
∆a1 = 0
∆a2 = (A f  p(ε))Ws = (l(ε) + p(ε))Ws
(17)
where Ws = T s2/(6tfht sh)
Adding these to the coefficients in Eq. 2 we get
p'(x') = a0 + a1x' + a2(3/2x'2  1/2) + a3(5/2x'3  3/2x') (18)
+ a4(35/8x'4  30/8x'2 + 3/8) + a5(63/8x'5  70/8x'3 + 15/8x')
+ p(ε)Ws(3/2x'2  3/2)  AfWs(3/2x’2  3/2)
since tfh ≈ t sh ≈ T s/2 and ε ≈ 0 looking at three special cases gives:
p’(ε) = p(ε) + 2/3 p(ε)(3/2ε2  3/2) + 2/3 l(ε)(3/2 ε2  3/2) ≈ l(ε)
p’(1) = p(1)
p’(1) = p(1)
(19)
From this we see that the quadratic correction does not change the scan mirror profile at
the beginning and end of scan where it is already constrained to be zero, and that any
mirror profile deviation at midscan is suppressed by the correction quadratic. Therefore,
any update to the mirror profile should obey the beginning and end of scan constraints:
∆p(1) = 0
∆p(1) = 0
(20)
Also, any update to the mirror profile at midscan will be suppressed by the quadratic
correction which implies that any errors in the mirror profile at midscan can be neither
measured nor corrected.
To account for these three special cases in our mirror correction model, we apply three
constraints to the correction profile, ∆p(x’), and treat them as additional observations:
∆p(1) = a 0 + a1 + a2 + a3 + a4 + a5 = 0
∆p(1) = a 0  a1 + a2  a3 + a4  a5 = 0
∆p(ε) = a0 + a1ε + a2(3/2ε2  1/2) + a3(5/2ε3  3/2ε)
+ a4(35/8ε4  30/8ε2 + 3/8) + a5(63/8ε5  70/8ε3 + 15/8ε) = 0
Version 3.2
133
(21)
7/9/98
IAS Geometric ATBD Version 3
where the ai's are now the corrections to the nominal profile. These constraints have the
practical effect of removing the low order terms, a0 , a1, and a2 from the correction model.
The first two constraints can be reformulated and simplified
a0 + a2 + a4 = 0
a1 + a3 + a5 = 0
(22)
Since the last constraint applies at ε which is very close to zero, it can be approximated
by:
a0  1/2a2 + 3/8a4 = 0
(23)
These constraints can be rearranged to determine the low order coefficients, a0 , a1 and a2
from the other three coefficients:
a0 = 7/12a4
a1 = a3  a5
a2 = 5/12a4
(24)
A fifth order polynomial which obeys these constraints can be used to model the higher
order mirror nonlinearity which cannot be removed by the midscan correction quadratic.
The measured deviations, converted to mirror angle deviations as a function of time,
provide observations of the ∆p(x’) mirror correction function. A weighted least squares
technique is used to fit a Legendre polynomial like Eq. 2 to these
The constraints are then applied to the coefficients, computed by the weighted least
squares procedure.
The user selects the weight for the observations based on the expected reference image
and correlation accuracy.
The existing mirror polynomial coefficients are converted to Legendre space according to
the transformation in Eqs. 4. If the higher order terms of the determined error polynomial
are deemed by the operator to be significant and reproducible, then they are added to the
existing coefficients and transformed back into mirror space according to Eqs. 3.
Inputs
As input, the mirror calibration algorithm requires a precision corrected grid, and a terrain
corrected image. It also requires the same digital elevation model which was used in the
terrain correction algorithm and the terrain correction table it generates. In order to
calculate the time in scan for each sample, it also needs to have the forward sensor model.
Finally it requires a reference image in the same framing as the terrain corrected image.
The current plan is to use DOQ imagery as the reference. Also, the PCD flags which
Version 3.2
134
7/9/98
IAS Geometric ATBD Version 3
indicate the scan mirror electronics (SME) number and scanning mode (Scan Angle
Monitor or Bumper) to know which mirror profile is applicable.
Outputs
The scan mirror calibration algorithm generates a set of fifth order Legendre polynomial
coefficients for corrections to the along and across scan mirror profiles for both forward
and reverse scans. If the deviation from the existing profile is deemed to be significant,
then the existing higher order mirror coefficients will be transformed into Legendre space,
have the deviations added to them and be converted back into mirror coefficients. The
lower order terms can be kept for characterization, though as shown above, they cannot
be uniquely determined and therefore should not be used for calibration. In addition, all of
the data generated by the correlation routines will also be stored.
Correlation Testing
One of the key components of the mirror calibration procedure is selecting an accurate
reference image source. Digital orthophoto quad data is 1m resolution, panchromatic,
aerial imagery which has been terrain corrected and rectified to a standard UTM
projection with an accuracy specification of 6m, 1 σ. They are available in 7.5 arcmin
quadrangles or 3.75 arcmin quarterquads. In order to use DOQs as a reference they need
to be downsampled to the same resolution as the ETM+ imagery and mosaicked to cover
several scans. In the case of Landsat 7 this means downsampling to the level of the
panchromatic band.
Ideally, we would like to mosaic the DOQ tiles at the highest resolution possible in order
to better blend the radiometric differences which may be present. Practically speaking,
however, the full resolution DOQs are too large to mosaic more than a few at a time. For
example, in one of the potential test sites (Landsat 4/5 path 28, row 30) the quarter quad
images are ~7700 lines by ~6900 samples and it takes ~200 of them to completely cover
the first 20 scans of the Landsat scene (see Fig. 3.1.526). As a compromise, the DOQs
will likely be mosaicked and reduced in resolution in stages.
For testing purposes, 16 quarterquads from four full quadrangles were lowpass filtered
and downsampled to 8m pixels. They were then mosaicked together and reprojected to
the SOM projection in a Landsat Pathorientation at a factor of the final desired pixel size
(7.5m). The mosaicked image was then lowpass filtered again and downsampled twice
more to the final 30m pixel size of the test image.
Version 3.2
135
7/9/98
IAS Geometric ATBD Version 3
Figure 3.1.5 26 DOQ quarterquad images covering about 20 TM scans
Since a precisionterrain corrected image was not available, a section of a Landsat 5
systematic image was correlated with the assembled DOQ imagery using hand picked
points. The TM image used for the correlation was a synthetic panchromatic image (at
30m pixel size) generated by a linear combination of bands 1, 2, 3 and 4. Using the
correlation offsets, a polynomial transformation was generated to remap the TM subscene into the frame and projection of the DOQ. Another set of features was hand picked
and the two images correlated. Without any sort of outlier rejection, the standard
deviation of the correlation was ~0.25 pixels. Using a regular grid of correlation points,
the standard deviation was ~0.35 pixels, again without any outlier rejection. Error
analysis indicates that a correlation accuracy of 0.25 pixels should be sufficient. From
these lessthanideal tests, it appears as though this goal will be achievable particularly if
the real panchromatic band ob the ETM+ is to be used.
A consideration to keep in mind is that the DOQ reference images are generally taken in
midspring. Therefore it would be best to choose test sites where seasonal variability will
have the least impact.
Polynomial Fitting
In order to test the ability of the algorithm to determine corrections to the scan mirror
coefficients, two systematic images were generated, one with the fifth order mirror
correction polynomial applied and one without. This should result in being able to back
out the original higher order coefficients. The two images were correlated using the first
96 scans. Each scan was broken down into 256 overlapping windows measuring eight
lines by 128 samples in output space. The inverse grid was used to make certain that the
correlation windows stayed within scans in input space.
The correlated (line, sample) pairs in the image and reference output space were run
through the inverse grid to find their corresponding (line, sample) in input space. In
practice, it is at this point that the terrain adjustment would be added to the input
samples. Since neither image had any terrain correction applied, it was not necessary for
this test. Knowing the input (line, sample) of the image and having initialized the forward
model, we can get a time in scan corresponding to each correlation point. We can also get
Version 3.2
136
7/9/98
IAS Geometric ATBD Version 3
the actual scan time, Ts, for each scan. These will allow us to determine the auxiliary
variable, x', as shown in Eq. 8. The difference between samples/lines in the two TM
images (in practice, a terrain corrected ETM+ image and a reference) in input space is
assumed to be an error which may be correctable by adjusting the coefficients of the scan
mirror profile.
It is this difference as a function of the auxiliary variable, x', which we attempt to fit with
a Legendre polynomial. In practice we would apply the constraints in Equation 19,
though this was not done for this test.
Error in Forward Scan for Higher Order Legendre Coefficients
3
2
Mirror Angle, ( µ rad)
1
0
1
1
0.5
0
0.5
1
2
3
4
5
6
7
Legendre Auxiliary Variable, x'
Image
Reference
ImageReference
Figure 3.1.5 27 Legendre Differences
The difference between the higher order Legendre terms (3rd, 4th, and 5th orders) of
the actual and estimated forward mirror profiles.
Version 3.2
137
7/9/98
IAS Geometric ATBD Version 3
Error in Forward Scan for Higher Order Mirror Coefficients
3
2
Mirror Angle, ( µ rad)
1
0
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
1
2
3
4
5
6
Time, t, (sec)
Image
Reference
ImageReference
Figure 3.1.5 28 Difference Between Actual Mirror Profiles
The difference between the higher order terms of the actual (reference) and
estimated (image) forward mirror profiles.
Figures 3.1.527 and 3.1.528 show the errors in the higher order terms of the correction
polynomial in Legendre terms and mirror terms respectively. In this test case, we expect
higher order terms of the polynomial we find to be the same as the actual along scan
mirror coefficients used in the base image. As we can see from the figures, the algorithm
was able to back out the polynomial to better than 1 µrad over most of the scan and
better than 2 µrad in all locations. The mirror profile repeatability specification for the
Landsat 7 ETM+ is 1.75 µrad 1σ, so these measurements using 42.5 µrad IFOV TM data
are of the same order as the random mirror variations.
3.1.5.5 Geometric Characterization Algorithms
3.1.5.5.1 Band to Band Registration
The Band to Band Registration Assessment segment of the Landsat 7 Image Assessment
System (IAS) is used to check the registration quality between bands of the same image.
Matching Resolutions
Version 3.2
138
7/9/98
IAS Geometric ATBD Version 3
Due to the difference in pixel sizes between bands, the higher resolution data will need to
be reduced so that it matches that of the lower resolution. The band to band
characterization will then be performed at the lower resolution. The higher resolution
imagery will be reduced using the Gaussian Pyramid approach [Reference 25]. This
approach works by low pass filtering the image data. The low pass filtering is
accomplished by convolving the image with an appropriate kernel or set of weights.
Applying this kernel will produce an image at a lower resolution than that of the original.
The lower resolution image can then also be reduced in size. The process takes on the
following form:
f r ( j, k ) =
2
2
∑ ∑ w( m, n) f (2 j + m,2 k + n)
j=1,…,M/2 k=1,…,N/2
m =−2 n =−2
Where f is the original image data, f r is the reduced subsampled image, and w(m,n) are
the low pass filter coefficients. Note that in the following case the imagery is reduced in
size by one half in both the line and sample direction. The low pass filter coefficients
can be represented as a set of separable one dimensional weights and are determined from
the following representation [Reference 25]:
w(2) = 1/4  a/2
w(1) = 1/4
w(0) = a
w(1) = 1/4
w(2) = 1/4  a/2
where a typically takes on a value between 0.3 and 0.6. More information on this
technique can be found in Reference 25.
In the case of Landsat 7 data this process will be used to reduce the 15 meter PAN band
to the 30 meter multispectral bands and to the 60 meter thermal band. It will also be used
to reduce the 30 meter multispectral bands to the 60 meter thermal band. In the case of
Landsat 4&5 test data, the 30 meter multispectral band will be reduced to the 120 meter
thermal band.
Correlation Routine
Normalized cross correlation will be used to produce a discrete correlation surface for an
array of points between two bands. This correlation routine helps in removing any gray
scale differences between two images. This will help to alleviate some of the problems
with the different radiometric responses between the bands. The correlation routine has
Version 3.2
139
7/9/98
IAS Geometric ATBD Version 3
the following formula when correlating between two image windowed subsets f and g,
each having a length and width of M and N respectively:
−
−
f
(
j
,
i
)
−
f
g
(
x
+
j
,
y
+
i
)
−
g
∑ ∑
j =− N / 2 i =− M / 2
N /2
R ( x, y ) =
M /2
[(
N / 2 M / 2
∑ ∑ f ( j, i) − f
j = − N / 2 i = − M / 2
)]
2
[(
N / 2 M / 2
∑ ∑ g ( x + j, y + i) − g
j =− N / 2 i =− M / 2
)]
2
1
2
where
−
N /2
−
1
g=
( M + 1)( N + 1)
−
f =
M /2
∑ ∑ g(x +
j, y + i)
j =− N / 2 i =− M / 2
1
( M + 1)( N + 1)
N /2
M /2
∑ ∑ f ( j, i)
j =− N / 2 i =− M / 2
Once the correlation surface is generated the peak of that surface must be determined.
This peak will represent the geometric shift between two windowed data sets and will be
calculated to a subpixel location.
Peak Finding Method
A subpixel location for the peak is found by using a polynomial fit method. This method
works by fitting a 2nd order polynomial around the peak of the discrete correlation
surface. The polynomial will represent the surface associated with a 3x3 area around the
peak. The discrete peak of the correlation surface is found first, then this peak and its
corresponding eight neighbors are used for fitting the polynomial. A least squares fit can
be used to calculate the polynomial coefficients once the 3x3 area and corresponding data
points are found. From this polynomial the location of the peak can be solved
analytically to a subpixel location.
Version 3.2
140
7/9/98
IAS Geometric ATBD Version 3
The polynomial equation has the following form:
P( x, y) = a 0 + a 1 x + a 2 y + a 3 xy + a 4 x 2 + a 5 y 2
The coefficients of P( x, y ) can be found by solving the matrix equation
[Y]=[X][a]
Where
[ Y ] is a 9x1 matrix corresponding to the eight correlation points found around,
and including, the peak.
[ X ] is a 9x6 matrix with each row consisting of:
1
xi
yi
xi y i
x2i
y 2i
where xi and y i are the x and y locations of the corresponding peak i in [ Y ]
[ a ] are the polynomial coefficients
The solution of [ a ] is then
[ a ] = ( [ X ]T [ X ] )1 [ X ] T [ Y ]
Where T refers to the transpose of a matrix and 1 refers to the inverse of a matrix.
Once the coefficients are solved for the subpixel peak can be found. The maximum
correlation peak can be found by taking the partial derivative of P( x,y ) with respect to x
and y:
∂
P ( x, y) = a 1 + a 3 y + 2a 4 x
∂x
∂
P ( x, y) = a 2 + a 3 x + 2a 5 y
∂y
By setting the partial derivatives to zero and solving for x and y the subpixel peak
location can be determined. This will lead to the following solution for the x and y
offsets:
x _ offset =
2a 1a 5 − a 2 a 3
2
a 3 − 4a 4 a 5
y _ offset =
2a 2 a 4 − a 1a 3
2
a 3 − 4a 4 a 5
Version 3.2
141
7/9/98
IAS Geometric ATBD Version 3
Outlier Detection
After computing a number of subpixel shifts between two bands, the data is checked for
outliers. These outliers are then removed from the data set. There are two steps involved
in detecting outliers. The first step involves looking at the measured offsets and
comparing these with a threshold established by an error budget associated with the
Landsat 7 satellite and the correlator. If any measurements are greater then this threshold
it is removed as an outlier. If a large percentage of the measurements are thrown out at
this point then the satellite is either behaving unexpectedly according to specifications or
none of the imagery corresponding to where the measurement was taken is of good
quality, a possibility for this would be a large amount of cloud cover. The second step
involves comparing the measured value with the bias estimate of all the measurements left
after applying the first threshold.
The first threshold is found from the following:
µrad = microradians
The known error associated with the satellite and correlator are:
Scan Mirror Repeatability 3.5 µrad
Field Angle 3.6 µrad
Jitter 1.45 µrad (spacecraft)
Jitter 0.65 µrad (ETM+)
Processing error 2.84 µrad
This leads to a total error for each band of
(35
. ) 2 + (3.6) 2 + (16
. ) 2 +(2.84) 2 = 6.0 µrad
The total error associated with two bands is then
(6.0) 2 + (6.0) 2 = 8.5µrad
Converting the microradian to pixels gives approximately 0.2 pixels. (The correlator
error budget is 0.25 pixel)
The total error budget is then
(0.2) 2 + (0.25) 2 = 0.32 pixels
The maximum band alignment error 0.3 pixel
The total measurement error is then
Version 3.2
142
7/9/98
IAS Geometric ATBD Version 3
(0.3) 2 + (0.32) 2 = 0.44 pixels
The threshold is then chosen as three times the total measurement error budget
3 * 0.44 = 1.32 pixels.
This value, 1.32 pixels, represents the first threshold that is to be applied to the data. All
measurements above this threshold are considered to be outliers and are removed.
The second threshold is found by calculating the average or bias and standard deviation
for the data set left after applying the first threshold outlier detection. The second
threshold is any point that is three standard deviations above the measured bias or
average. The data set left after applying the second threshold outlier detection represents
the valid measurements taken.
Statistical Measure
Once the outliers have been detected and removed the resulting data is used to calculate a
bias or average band offset between two bands. The standard deviation of the data set is
also calculated.
1 N
avg = ∑ m i
N i=0
(std ) 2 =
1 N
∑ ( m i − avg)
N − 1 i=0
2
where
avg = average
std = standard deviation
mi = measured offset
N = Number of valid measurements
Test Simulation
A test simulation was generated to determine the accuracy of the correlation and peak
finding routines by using an image of path 28 row 30 and resampling it with various
geometric offsets. A 16x16 truncated sinc [sin(x)/(x)] resampling function was used to
produce geometric shifts. Correlation windows of 512x512 were chosen. Sixteen windows
were chosen at random for measurement. The various images of geometric offsets were
then correlated with the original. This gave an indication of the accuracy of the correlator
itself. The results are listed below:
Version 3.2
143
7/9/98
IAS Geometric ATBD Version 3
Offset
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Average Offset
0.0
0.089817
0.152810
0.236563
0.344935
0.485531
0.632771
0.745469
0.831895
0.903604
Standard Deviation
0.0
0.003823
0.001479
0.002199
0.002791
0.006224
0.008459
0.008163
0.006718
0.004832
This indicates that the correlator works to within onetenth of a pixel. This represents an
ideal situation in that the image is being correlated to itself. Through work with different
band combinations an accuracy of 0.25 pixels can be expected across bands from the
correlator. A more thorough explanation of errors associated with correlation can be
found in Reference 21. The Image to Image characterization algorithm also performed a
similar test simulation; this can also be used as a reference on the accuracy of the
correlator.
Measured Landsat 4 & 5 Data
Correlation will be done in such a manner that all bands can be either directly or indirectly
referenced to each other. The test site that is chosen to perform band to band registration
should have features that are present throughout all bands. It was found that desert
scenes with little vegetation presented the best scenario. A scene corresponding to path
41 row 36 was chosen as a test site. Correlation windows were chosen at uniform
increments throughout the image. When correlating the multispectral bands an offset of
1000 lines was chosen from the top and bottom and 1000 samples from the right and left
sides of the image. Windows were then chosen at 500 lines and 500 sample increments
throughout the image. The window size chosen was 256x256. When correlating two
bands with one being the thermal band offsets of 250 lines, 250 samples, increments of
125 lines, 125 samples, and window sizes of 128x128 were chosen. The following values
were then calculated:
X Offsets
(Offsets are in terms of multispectral pixels)
Reference
1
1
1
1
1
2
Version 3.2
Search
2
3
4
5
7
3
average
0.003244
0.001516
0.025383
0.261743
0.252938
0.000870
144
standard deviation
0.010120
0.021785
0.113919
0.061539
0.053082
0.012164
7/9/98
IAS Geometric ATBD Version 3
2
2
2
3
3
3
4
4
5
4
5
7
4
5
7
5
7
7
0.014538
0.266884
0.249701
0.035892
0.267662
0.250002
0.227007
0.194449
0.012561
0.134436
0.055054
0.047483
0.079654
0.050735
0.048614
0.083736
0.145899
0.020213
average
0.786044
0.788873
0.808790
0.758903
0.757974
0.754215
standard deviation
0.177695
0.164740
0.159983
0.194022
0.108728
0.165351
(Offsets in terms of thermal pixels)
Reference
1
2
3
4
5
6
Search
6
6
6
6
6
7
Y OFFSETS
(Offsets in terms of multispectral pixels)
Reference
1
1
1
1
1
2
2
2
2
3
3
3
4
4
5
Search
2
3
4
5
7
3
4
5
7
4
5
7
5
7
7
average
0.003527
0.014177
0.024184
0.366080
0.382187
0.010344
0.014033
0.362300
0.374961
0.006703
0.350089
0.365271
0.319364
0.328139
0.007921
standard deviation
0.009219
0.019110
0.060681
0.062386
0.051700
0.011238
0.075062
0.047765
0.044400
0.067899
0.034363
0.036863
0.073798
0.097649
0.014011
average
0.737981
0.726864
0.718318
0.747355
0.633007
standard deviation
0.177695
0.197395
0.169007
0.191474
0.095764
(Offset in terms of thermal pixels)
Reference
1
2
3
4
5
Version 3.2
Search
6
6
6
6
6
145
7/9/98
IAS Geometric ATBD Version 3
6
7
0.611179
0.153688
It should be noted that the measured between band alignment does not meet the
specifications for the Landsat 7 system. This data will be used by the Focal Plane Band
to Band Calibration algorithm to relate these parameters back to the focal plane to adjust
the focal plane alignment, if deemed necessary.
3.1.5.5.2 Image to Image Registration
The Image to Image Assessment segment of the Landsat 7 Image Assessment System
(IAS) is used to check the registration quality between two images. The two images are
assumed to be from the same band and same (path, row) taken at different times,
processed to the same level either precision or terrain corrected products.
Accuracy Requirements
According to the Landsat 7 IAS Requirements Document, “The IAS shall provide the
capability to perform image to image registration to an accuracy of 0.4 multispectral
sensor GSD, 0.9p, in the alongtrack and crosstrack directions providing all inputs are
within specification.” The understanding is that this requirement refers to two
temporally spaced images processed to the same level using the same ground control.
Furthermore, since the accuracy requirements for systematically processed images are on
the order of ~250m (400m, 0.9p), this requirement is assumed applicable only to
precision and/or terrain corrected images.
As for requirements placed on the image to image correlation technique, “The IAS shall be
capable of digitally correlating common features in the…same bands of separate images to
an accuracy of 0.1pixel, 0.9p.” This is equivalent to requiring that the correlation
technique be capable of backing out a known offset to within 0.1 pixel, 0.9p (~0.06 pixels
at one sigma).
Procedure
The IAS Image to Image Assessment segment provides two options for selecting imagetoimage points. Correlation points can be generated either by creating a regular grid of
hundreds of points, or by manually selecting fifty to a hundred well distributed, “humanrecognizable” ground control features. Given that images only need to be processed to the
full precision and terrain corrected levels for a small number of test sites, the latter is the
more likely operational implementation. Using the grid of points, however, may be useful
in geometric analysis.
Version 3.2
146
7/9/98
IAS Geometric ATBD Version 3
The normalized cross correlation technique is defined as
−
−
f
(
j
,
i
)
−
f
g
(
x
+
j
,
y
+
i
)
−
g
j = − N / 2 i = − M / 2
N /2
M /2
∑ ∑
R ( x, y ) =
[(
N / 2 M / 2
∑ ∑ f ( j, i) − f
j = − N / 2 i = − M / 2
)]
2
[(
N / 2 M / 2
∑ ∑ g ( x + j, y + i) − g
j =− N / 2 i =− M / 2
)]
2
1
2
where
−
f =
−
g=
1
( M + 1)( N + 1)
1
( M + 1)( N + 1)
N /2
M /2
∑ ∑ f ( j, i)
j =− N / 2 i =− M / 2
N /2
M /2
∑ ∑ g(x +
j, y + i)
j =− N / 2 i = _ M / 2
The resulting correlation surface is searched for a peak. In order to determine the subpixel location of the true correlation maximum, a surface is fit to the peak and its
neighbors. The subpixel location of the peak value for the bestfit function is determined
analytically and stored. The correlation functions included options for fitting parabolic
and Gaussian surfaces using a 5 x 5 neighborhood of pixels around the peak. The
correlation function was modified to use a 3 x 3 neighborhood as testing showed these
results to be more accurate and reliable. This is discussed in greater detail below.
The verification routine uses the points and offsets determined by the correlation function
to calculate statistics. For each correlation point, it reports the line, sample and total
errors in the search image relative to the reference. It also calculates the mean, standard
deviation, median, range, and RMSE for those errors and determines the number (and
percentage) of points not meeting user specified accuracy requirements. All of these
values will be stored for trending.
Version 3.2
147
7/9/98
IAS Geometric ATBD Version 3
In order to assist geometric assessment, the statistics compiled by verification routine can
be plotted. Useful plots include, but are not limited to: surface plots of line or sample
error versus line and sample, vector plots showing line and sample error distributed on the
image plane, or color coded residuals overlaid on an image.
Some outlier detection and rejection is performed within the correlation function. The
routine flags and rejects correlation points if they are: too near the edge of the data
window, similar in magnitude to subsidiary peaks, or below a userspecified strength
threshold. There is also a flag which rejects points if the diagonal distance between the
reference and search subimages is greater than a userspecified value.
Whether further outlier rejection will be used on the final residuals is yet to be determined
since initially, analysis will primarily be concerned with determining whether outliers
indicate a problem in the data or a problem in the system. An exception flag should be
raised if the line and sample error mean and standard deviations are in excess of the
requirement value specified above.
Inputs
Image to image assessment requires a reference image which has been processed to the
precision or terrain corrected level and a search image processed to the same level with the
same ground control, (and the same DEM in the case of terrain corrected images). It also
requires a set of preselected image features or grid of whose locations are refined by the
correlation function.
Outputs
The correlation function generates a file containing the points and offsets which is used
by the verify routine to compile statistics. The verify routine generates a report file
containing information about all the points used for correlation including line, sample and
total errors. It also lists overall statistics including mean, standard deviation, median and
RMSE for line, sample and total errors and the number of points not meeting user
specified error thresholds. The information in the report is used for trending .
Algorithm Testing
To test the accuracy and precision of the correlation routine, a base image using band 1 of
a Landsat 5 test scene (path 28, row 30) of northwestern Iowa was generated. Artificial
shifts in the sample direction were then incorporated into the image in 0.1 pixel
increments using cubic convolution resampling accurate to 1/32 of a pixel. Image to image
correlation was then performed with artificial offsets ranging from 0.0 to 1.0 pixels using a
set of 51 manually picked points. The resulting plots of the mean and standard deviation
Version 3.2
148
7/9/98
IAS Geometric ATBD Version 3
values are shown in Figures 3.1.529 and 3.1.530. At an offset of 0.5 pixels, the
correlation routine using a 5 x 5 neighborhood about the correlation peak became
unreliable giving a bimodal error distribution surrounding the true offset. A possible
explanation for this peculiar behavior is that the objects used in the correlation are
typically linear features approximately one pixel wide (e.g., roads) surrounded by
relatively uniform backgrounds (e.g., fields). In these cases we would not expect a
parabola to be a good fit beyond one pixel to either side of the peak.
Correlation Error for 3 and 5 point
Parabolic and Gaussian Peak Fits
0.15
0.1
Prediction Error, (pixels)
0.05
3pt Parabolic Error
5pt Parabolic Error
0
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
3pt Gaussian Error
5pt Gaussian Error
0.05
0.1
0.15
True Displacement, (pixels)
Figure 3.1.529 Correlation Errors
Correlation errors using 3x3 and 5x5 pixel neighborhoods, and parabolic and
Gaussian peak surface fits for an artificially induced sample offset.
Version 3.2
149
7/9/98
IAS Geometric ATBD Version 3
Standard Deviation for 3 and 5 point
Parabolic and Gaussian Peak Fits
0.14
0.12
Standard Deviation, (pixels)
0.1
0.08
3pt Parabolic StDev
5pt Parabolic StDev
3pt Gaussian StDev
5pt Gaussian StDev
0.06
0.04
0.02
0
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
True Displacement, (pixels)
Figure 3.1.5 30 Correlation Standard Deviations
Standard deviation of the errors 3x3 and 5x5 pixel neighborhoods, and
parabolic and Gaussian peak surface fits for an artificially induced sample offset.
As an example, consider the case of a bright, onepixelwide road cutting through a
uniformly dark field. If this road were displaced to the right by 0.5 pixels, the resultant
image would have the road represented by two pixels of medium brightness surrounded
again by the dark field. This situation is shown schematically in Figure 3.1.531 as a one
dimensional cross section. Figure 3.5.132 shows the cross correlation (nonnormalized)
of these two data sets and the resulting parabolas which occur when trying to fit the data
with five points using both the left and right maxima as the center point. The resulting
peaks of the parabolas in Fig. 3.1.532 occur at 7.27 and 7.73 rather than at 7.5 where the
true peak should occur. So while the average value is correct, the standard deviation for
this case is unacceptable as are shown in Figs. 3.1.529 and 3.1.530.
By inspection of Fig. 3.1.531, it is evident that fitting a parabola using three points
rather than five, the peak will occur at 7.5 regardless of which maxima is taken as the
center point. Figure 3.1.529 shows the three point method to be more accurate than the
five point method. For these reasons, the correlation routine was modified to use a 3 x 3
neighborhood around the peak rather than a 5 x 5 neighborhood.
By looking at the shape of the cross correlation data in Fig. 3.1.532, one would expect
that a Gaussian curve might be a better fit than a parabola, thus giving reduced errors.
Version 3.2
150
7/9/98
IAS Geometric ATBD Version 3
The results of doing so are also shown in Fig. 3.1.529 and indeed they lead to a ~20%
reduction in error at the peak locations. The sinusoidal shape of the error curve shown in
Fig. 3.1.529 is consistent with the theoretical errors determined by Dvornychenko
(Reference 21).
10
9
8
Brightness
7
6
5
4
3
2
1
0
0
1
2
3
4
5
6
7
8
9
10
11
12 13
14
Sample
Original Road
Road Shifted by 0.5 pixel
Figure 3.1.5 31 Road Cross Section Example
Schematic crosssection of a bright road running through a dark field and
how it would appear when shifted 0.5 pixel to the right
Version 3.2
151
7/9/98
IAS Geometric ATBD Version 3
60
Cross Correlation
50
40
30
20
10
0
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
10
Sample
Cross Correlation
Left Peak
Right Peak
Figure 3.1.5 32 Correlation Errors
Cross correlation of original and shifted data and fitting parabolas to the
data using five point centered on either the left or right maxima.
Further Investigations
While it appears as though the accuracy of the modified correlation routine is sufficient to
meet the requirements, it must be noted that these tests were run under ideal conditions,
i.e., the same image data (although shifted and aliased by using a cubic convolution
resampler accurate to 1/32 pixel). Under real circumstances where there may be
significant changes in the scenes from image to image, the standard deviations can be
expected to increase.
Normalized crosscorrelation is designed to take some image differences into account and
has been used successfully in projects involving Landsat MSS and TM imagery.
However, it has not been used to measure the registration accuracy under such stringent
requirements as in the IAS.
Further investigations may include widening the reference chips to greater than 64x64.
This should help by having more data available and thus smoothing out the correlation
surface. Also, other techniques such as edge or phase correlation may be useful when
correlating images with differing spectral data (e.g., data from different seasons).
Version 3.2
152
7/9/98
IAS Geometric ATBD Version 3
3.1.5.5.3 Geodetic Accuracy
This section describes the algorithm, input data, output data, and proposed usage of the
output data applied in the Geodetic Characterization Module for the Landsat 7 Image
Assessment System.
The Algorithm
In the system data flow chart, this algorithm immediately follows the precision correction
solution module. Thus, this description should be read in the context related to the
description of precision correction solution module. The GCP data quality requirement
applies to the input data of precision correction solution module, and is based on the
assumption that GCP marking and correlation are done to 0.1 pixel accuracy.
Purpose:
The purpose of this algorithm is to verify the absolute accuracy of the systematic and
precision products by analyzing the Ground Control Point measurement residuals before
and after the precision correction solution.
Process:
This algorithm applies standard statistical error analysis on the residual file of the
precision correction solution. The residual file contains residuals for GCPs before and
after precision correction. The algorithm takes one input file (the residual file), and
generates two output files (the record file and the outlier data file). The following steps
are performed:
1. Calculate and write to the record file the following 3 values:
the total number of GCPs in the residual;
the mean latitude of all GCPs;
the mean longitude of all GCPs;
2. Calculate and write to the record file the following 5 values for the prefit residuals:
Mean value of the crosstrack residuals;
Mean value of the alongtrack residuals;
RMSE about the mean value for crosstrack residuals;
RMSE about the mean value for alongtrack residuals;
Correlation coefficient between along and crosstrack residuals;
3. Calculate and write to the record file the following 5 values for the postfit residuals:
Mean value of the crosstrack residuals;
Mean value of the alongtrack residuals;
Standard deviation for crosstrack residuals;
Version 3.2
153
7/9/98
IAS Geometric ATBD Version 3
Standard deviation for alongtrack residuals;
Correlation coefficient between along and crosstrack residuals;
4. Compare the postfit residuals with their Standard Deviations, mark those GCPs
with > 2 sigma residual as outliers;
5. Write to the outlier data file the information for each outlier GCP;
The Input Data
This module requires only one input data file, which is the output residual file by the
Precision Correction Module. This file contains the residual information from Iteration 0,
which is the prefit residual, through the Final Iteration, which is the postfit residual.
GCP Data Accuracy Requirement
Following properties of the GCP measurements are required to produce a meaningful
precision correction solution and geodetic characterization result:
• GCP coordinate should be in a system consistent with the that of satellite ephemeris.
WGS84 is preferred, NAD83 in North America is acceptable.
• The accuracy of the GCP coordinate should be on the level of 3 meters in horizontal
and 20 meters in elevation to guarantee 0.1 pixel accuracy;
• Number of GCPs: A minimum number of 15 valid GCPs is required for reliable
characterization (6/number of precision correction parameters + 9/reasonable number
degree of freedom). To allow certain portion of outlier editing, 25 to 50 input GCPs
are suggested.
Usage of the Output Data
•
•
•
•
The outlier information file will be used to remove the GCP from both the GCP data
file and the EPH data file for precision correction solution (see algorithm description
for precision correction solution module). A new round of precision correction
solution and geodetic characterization will be performed after the editing;
The record file will be used as a measure of the absolute accuracy of the systematic
product (prefit residual statistics) and precision product (postfit residual statistics).
Trending analysis can be done on the parameters provided in the record file.
Error characterization for check point data (GCPs not used for the precision
correction solution) will be done in a different round after the systematic model is
updated with the precision correction parameters.
Version 3.2
154
7/9/98
IAS Geometric ATBD Version 3
3.1.5.5.4 Geometric Accuracy
For the purposes of determining geometric accuracy, three tests are to be performed:
1. Visual inspection. The output image will be inspected visually for any geometric
distortions. These distortions could include scan misalignment, artifacts produced by the
resampling process such as ringing, or detector (line to line) discontinuity.
2. Plot of the residuals of the ground control points and their corresponding
locations in the terrain corrected image. This can be either a vector plot indicating the
lines/sample offsets or an image with the residuals color coded according to their
magnitude superimposed on a paper product of the image.
3. A polynomial fit involving the ground control points and their corresponding locations
in the terrain corrected output product.
Polynomial Fit
A polynomial fit will be used to estimate the amount of scaling in x and y, rotation, offset
in x and y, and nonorthogonality angle associated with a terrain corrected image. The
polynomial will take the form:
x o = a + b∗ x i + c∗ y i
y o = d + e∗ x i + f ∗ y i
Where
xo = output projection x coordinate (measured)
y o = output projection y coordinate (measured)
xi = input projection x coordinate (truth)
y i = input projection y coordinate (truth)
a,b,c,d,e,f transformation coefficients
The polynomial coefficients will be determined from a least squares fit of the x and y
coordinates of features located in the terrain corrected product and their corresponding x
and y coordinates for the same feature in the ground control source. The transformations
that “corrects” x o,y o to xi,y i is given by:
xi = S x ( x − x o ) cos(θ ) − S y ( y − y o ) sin(θ − α )
y i = S x ( x − x o ) sin(θ ) + S y ( y − y o ) cos(θ − α )
Version 3.2
155
7/9/98
IAS Geometric ATBD Version 3
The error or distortions can then be determined from the polynomial transformation
coefficients by the following:
e2 + f 2
Sx =
bf − ce
b2 + c2
Sy =
bf − ce
−e
theta = tan −1 ( )
f
− be − cf
alpha = tan −1 (
)
bf − ce
Sx = scale error in x direction
Sy = scale error in y direction
alpha = nonorthogonality angle
theta = rotation angle
a = offset in x direction
d = offset in y direction
a,b,c,d,e,f transformation coefficients
The scale error can be converted to pixels by the following:
(Sx1) * (image pixel size in x)
(Sy1) * (image pixel size in y)
Theta can be converted to pixels by the following:
theta * (size of image diagonal)
size of image diagonal = SQRT((image pixel size x)2 + (image pixel size y) 2 )
Alpha can be converted to pixels by the following:
alpha * (maximum pixel size)
maximum pixel size = MAX(image pixel size in x, image pixel size in y)
All of these values should be very close to zero except S x and Sy which should be close
to one. If any of these are not zero, or Sx and Sy are not close to one, then they represent
the amount of geometric distortions not corrected for by the system.
The residuals left after applying the polynomial coefficients to the terrain corrected x and
y coordinates will represent any higher order distortions. These values should also be
very small.
Version 3.2
156
7/9/98
IAS Geometric ATBD Version 3
Note: This algorithm requires a terrain corrected image produced from the IAS
system along with at least fifty well distributed ground control points. The
PAN band would give the highest precision when calculating errors, however
any band can be used in the process. All points, from both the image and the
control source, need to be in the same projection with the same datum and pixel
size. The control points and the points taken from the IAS terrain corrected
product should be to a subpixel accuracy. The absolute accuracy of the points
picked will depend on the method chosen for point matching.
3.1.6 Variance or Uncertainty Estimates
The following sections contain error budget analyses for the key Landsat 7 ETM+
geometric accuracy requirements. These analyses were used to predict the Level 1Gs
geodetic accuracy performance and to develop algorithm and support data accuracy
requirements for the geometric calibration and characterization activities.
3.1.6.1 Level 1 Processing Accuracy
The fundamental Level 1 Processing accuracy requirement is to be able to produce
systematically corrected (1Gs) products which are accurate to within 250 meters per
coordinate (one sigma) at nadir. This is stated in the IAS Element Specification as
requirement number 3.2.3.4. An error budget for 1Gs geodetic accuracy performance can
be built up using the expected accuracy of the data components which are part of the
Landsat 7 ETM+ geometric model. These components include spacecraft ephemeris and
attitude knowledge, spacecraft clock errors, knowledge and stability of the alignment of
the ETM+ instrument to the Landsat 7 spacecraft, and knowledge and stability of the
internal geometry of the ETM+ instrument. Accuracy bounds for most of these data
elements are specified in the Landsat 7 System Specification. Allocations for ground
processing and test point mensuration are included in the total geodetic accuracy error
budget.
Errors in each of the geometric model data components were assumed to be independent
zero mean random variables in the geodetic accuracy analysis. The relevant accuracy
values were extracted from the Landsat 7 System Specification, converted to 1 sigma
values if necessary, and then mapped to position errors in meters on the ground. For
along track and across track angular error values this mapping was done by multiplying
the angular error (in radians) by the nominal Landsat 7 height above ground of 705
kilometers. The spacecraft clock error was converted to ground meters by multiplying
the time error by a conservative estimate of the nominal Landsat 7 orbital velocity of 7.5
kilometers per second.
Version 3.2
157
7/9/98
IAS Geometric ATBD Version 3
Two cases are considered in the geodetic accuracy error budget: 1) the expected accuracy
atlaunch based on the prelaunch knowledge of the alignment between the Landsat 7
navigation base reference coordinate system and the ETM+ optical axis, and 2) the
expected accuracy in routine operations after the IAS has calibrated the instrument to
spacecraft alignment. The 250 meter accuracy specification applies to both the along
track and across track directions, but the presence of the spacecraft clock error makes the
along track direction the critical case. The atlaunch along track accuracy analysis is
presented in Table 3.1.6.11. The columns in the table contain the following information:
• Column 1: identifies the data element error source,
• Column 2: shows the error contribution mapped to ground meters at nadir (1
sigma),
• Column 3: shows the error contribution in the units in which it is specified (1
sigma),
• Column 4: identifies the specification units,
• Column 5: shows the original source of the specified value since some are rootsumsquare (RSS) combinations of lower level components,
• Column 6: shows where the error source was allocated at the system level (e.g., to
the spacecraft, instrument, or ground system).
Version 3.2
158
7/9/98
IAS Geometric ATBD Version 3
Source
Meters
1σ
154
Spec
1σ
45.0
Units
arcsec 45.0 arcsec (1 σ) (R,P,Y)
Satellite
Scan Mirror Repeatability
3
0.72
arcsec 1.75 microradians * 2 (1 σ)
ETM+
Field Angle
3
0.78
arcsec RSS of:
0.18 arcsec (1 σ) calibration
0.76 arcsec (1 σ) stability (0.2 pan)
IAS
ETM+
arcsec RSS of:
8.0 arcsec (1 σ) ETM+ to NBR
15.0 arcsec (1 σ) LOS to ETM+
Satellite
ETM+
arcsec RSS of:
720.0 arcsec (3 σ) prelaunch
60.0 arcsec (1 σ) prelaunch
Satellite
ETM+
Vehicle Attitude
Alignment Knowledge
(R,P,Y)  stability
58
Alignment Uncertainty
(R,P,Y)  calibration
846
Ephemeris (I,C,R)
133
Jitter (R,P,Y)
2
17.0
247.4
133.33
0.44
meter
Comments
400 meter (3 σ) from sum of:
375 m (3 σ) predict
20 m (3 σ) interpolation
Segment
FDF
Satellite
arcsec RSS of:
0.30 arcsec (1 σ) low frequency
0.30 arcsec (1 σ) high frequency
0.134 arcsec (1 σ) high frequency
Satellite
Satellite
ETM+
Satellite
Timing (along track)
38
5
msec
15 msec (3 σ) random error
Mensuration
15
15
meter
0.5 pixels (1 σ)
User
Processing Error
2
2
meter
Model budget
IAS
RSS Estimate
873
RSS Margin

Specification
250
Table 3.6.1.1 1 AtLaunch 1Gs Geodetic Error Budget
As Table 3.6.1.11 demonstrates the 250 meter geodetic accuracy requirement is not
expected to be met based on the atlaunch knowledge of the ETM+ instrument to Landsat
7 spacecraft alignment.
The Landsat 7 System Specification requires the IAS to estimate the onorbit alignment to
an accuracy of 24 arc seconds (as compared to the prelaunch accuracy of about 248 arc
seconds). This improvement in alignment knowledge accuracy allows the geodetic
accuracy budget to meet the 250 meter specification with significant margin. Table
Version 3.2
159
7/9/98
IAS Geometric ATBD Version 3
3.6.1.12 updates the along track accuracy analysis with the postcalibration alignment
knowledge accuracy. The table column definitions are the same as for Table 3.6.1.11.
Source
Meters
1σ
154
Spec
1σ
45.0
Units
arcsec 45.0 arcsec (1 σ) (R,P,Y)
Satellite
Scan Mirror Repeatability
3
0.72
arcsec 1.75 microradians * 2 (1 σ)
ETM+
Field Angle
3
0.78
arcsec RSS of:
0.18 arcsec (1 σ) calibration
0.76 arcsec (1 σ) stability (0.2 pan)
IAS
ETM+
arcsec RSS of:
8.0 arcsec (1 σ) ETM+ to NBR
15.0 arcsec (1 σ) LOS to ETM+
Satellite
ETM+
Vehicle Attitude
Alignment Knowledge
(R,P,Y)  stability
58
Alignment Uncertainty
(R,P,Y)  calibration
82
24.0
Ephemeris (I,C,R)
133
133.33
Jitter (R,P,Y)
2
17.0
0.44
Comments
Segment
arcsec 24.0 arcsec (1 σ) (R,P,Y)
meter
400 meter (3 σ) from sum of:
375 m (3 σ) predict
20 m (3 σ) interpolation
IAS
FDF
Satellite
arcsec RSS of:
0.30 arcsec (1 σ) low frequency
0.30 arcsec (1 σ) high frequency
0.134 arcsec (1 σ) high frequency
Satellite
Satellite
ETM+
Satellite
Timing (along track)
38
5
msec
15 msec (3 σ) random error
Mensuration
15
15
meter
0.5 pixels (1 σ)
User
Processing Error
2
2
meter
Model budget
IAS
RSS Estimate
231
RSS Margin
97
Specification
250
Table 3.6.1.1 2 1Gs Geodetic Error Budget After Alignment Calibration
3.1.6.2 Geometric Calibration Accuracy
The development of the three IAS geometric calibration algorithms included estimation
error analyses to ensure that the relevant geometric parameters could be estimated with
sufficient accuracy to justify updating the prelaunch knowledge in postlaunch releases
of the Calibration Parameter File. The results of these analyses are presented in the
following sections.
Version 3.2
160
7/9/98
IAS Geometric ATBD Version 3
3.1.6.2.1 Sensor Alignment Calibration Accuracy
As noted above, the IAS is required to calibrate the ETM+ instrument to Landsat 7
spacecraft alignment knowledge to an accuracy of 24 arc seconds in order to achieve a 1Gs
product geodetic accuracy of 250 meters. This is stated in the IAS Element Specification
as requirement number 3.2.3.12. The IAS sensor alignment algorithm uses the results of a
series of precision correction solutions to extract the systematic alignment bias from the
random scenetoscene attitude error, as described above. The goal of the sensor
alignment error analysis was to determine the number of calibration scenes required to
achieve the 24 arc second alignment knowledge accuracy as a function of the accuracy of
the supporting definitive ephemeris and ground control data.
The precision correction solution model was used to perform an error propagation
analysis to predict the accuracy of the pointing error (alignment plus attitude) estimates
produced by a single precision correction solution. This analysis was based on the
assumption that all inputs were accurate to their prelaunch specification values, and that
14 well distributed ground control points were successfully measured with an accuracy of
15 meters (ground position error plus image mensuration error) for each point. The input
accuracy assumptions and the resulting pointing error accuracy estimates are summarized
in Table 3.1.6.2.11.
Version 3.2
161
7/9/98
IAS Geometric ATBD Version 3
Input
Along Track Position
Across Track Position
Radial Position
Roll
Pitch
Yaw
Ground Control
Roll Estimate
Pitch Estimate
Yaw Estimate
1 σ error Rationale
63 meters RSS of 50 meter estimate of FDF definitive ephemeris accuracy
and 38 meter spacecraft clock error
50 meters Estimate of FDF definitive ephemeris accuracy
50 meters Estimate of FDF definitive ephemeris accuracy
252 arcsec RSS of 247 arcsec prelaunch alignment knowledge error, 17 arcsec
alignment stability error, and 45 arcsec attitude knowledge error
252 arcsec RSS of 247 arcsec prelaunch alignment knowledge error, 17 arcsec
alignment stability error, and 45 arcsec attitude knowledge error
252 arcsec RSS of 247 arcsec prelaunch alignment knowledge error, 17 arcsec
alignment stability error, and 45 arcsec attitude knowledge error
15 meters RSS of 12 meter position error (1:50,000scale map source) and
7.5 meter (0.5 panchromatic pixel) image mensuration error
13.1 arcsec A posteriori error estimate
16.6 arcsec A posteriori error estimate
16.7 arcsec A posteriori error estimate
Table 3.1.6.2.1 1 Pointing Error Estimate Accuracy Analysis  Inputs and Results
The roll estimate accuracy is somewhat better than the pitch and yaw estimates primarily
because the across track roll must be separated from the across track position error while
the pitch and yaw both must be separated from the along track position error which is
less accurate due to the spacecraft clock error.
The roll, pitch, and yaw pointing error estimates generated by the precision correction
process represent the combined effect of static alignment error, attitude error, and
dynamic alignment error. This can be expressed as:
pointing error observation = attitude + static alignment + dynamic alignment
It is the static alignment error which we seek to calibrate. Rearranging the previous
equation yields:
static alignment = pointing error observation  attitude  dynamic alignment
The attitude and dynamic alignment errors are assumed to be zero mean, so the best
estimate of the static alignment error is just the observed pointing error estimated by the
precision correction solution. The variance of the static alignment estimate is the sum of
the variance of the three terms on the right hand side:
σ2alignment = σ2observation + (45)2 + (17)2
where: 45 arc seconds = expected attitude standard deviation
Version 3.2
162
7/9/98
IAS Geometric ATBD Version 3
17 arc seconds = expected dynamic alignment standard deviation
Table 3.1.6.2.12 summarizes the resulting standard deviations for the three alignment
error components.
Component
Roll
Pitch
Yaw
Standard Deviation (Single
Estimate)
49.9 arc seconds
50.9 arc seconds
50.9 arc seconds
Table 3.1.6.2.1 2 Alignment Estimate Accuracy
As multiple test scene results become available the observed pointing errors will be
combined using a Kalman filter as described above. To a first approximation, the effect on
the resulting alignment estimate will be to improve the accuracy as the square root of the
number of test scenes processed. The resulting improvement in alignment knowledge
accuracy as a function of the number of calibration scenes processed is shown in Figure
3.1.6.2.11.
The figure suggests that the 24 arc second alignment accuracy requirement can be met
after five calibration scenes have been processed. The accuracy estimates based on this
relatively simple model of the alignment estimation process have been confirmed by
processing simulated data (with known alignment biases introduced) through the precision
correction and sensor alignment prototype software and analyzing the results. These tests
have shown that the alignment knowledge can be expected to reach the 24 arc second
accuracy threshold after five or six calibration test scenes have been processed.
3.1.6.2.2 Scan Mirror Calibration Accuracy
Although there is no formally specified accuracy requirement for the scan mirror
calibration process, an accuracy threshold of 1.75 microradians (in mirror space) was
selected as the goal for mirror profile estimation accuracy. This value corresponds to the
scan mirror repeatability specification called out in the Landsat 7 System Specification in
section 3.7.8.1.17.4. Determining the scan mirror profile to this accuracy will keep the
systematic mirror profile knowledge within the random variation of the scan mirror
motion.
Error propagation analysis (based on a set of relatively conservative assumptions
regarding the accuracy of the input data) was used to determine the number of ETM+
scans and the number of test points per scan which must be correlated with a DOQ
reference image to achieve a 1.75 microradian estimation accuracy. The input data and
imagetoimage correlation accuracy assumptions used are presented in Table 3.1.6.2.21.
Version 3.2
163
7/9/98
IAS Geometric ATBD Version 3
Source
DOQ Reference Image
ETM+ 1Gt Image
Correlation Error
Net Observation Error
Obs Error  Mirror Space
1 σ error
8.5 urad
21.25 urad
10.625 urad
25.23 urad
12.62 urad
A priori Coefficient Error
42.5 urad
Number of Points per Scan
60
Rationale
Corresponds to 6 meter DOQ accuracy spec.
Corresponds to 1 panchromatic pixel
Corresponds to 0.5 panchromatic pixel
RSS of DOQ, 1Gt, and Correlation Errors
One half the object space value above
Corresponds to 1 reflective band pixel
Corresponds to 1 per geometric grid cell
Table 3.1.6.2.2 1 Scan Mirror Calibration Accuracy Analysis Input Assumptions
The error propagation analysis assumed a set of test points which were well distributed
as a function of scan angle to form the normal equations which would be used to estimate
the mirror profile Legendre polynomial coefficients. The normal equations were then
inverted to form the predicted a posteriori error covariance matrix from which variance
estimates for the Legendre coefficients were extracted. These coefficient variances were
then summed to compute the net mirror profile uncertainty. The number of test points
used to construct the normal equation matrix was then varied until the net mirror profile
uncertainty was less than 1.75 microradians. Two separate cases were investigated:
1) estimation accuracy for an unconstrained fifth order polynomial with six unknown
coefficients; and 2) estimation accuracy for a fifth order polynomial constrained to meet
the beginning of scan, end of scan, and midscan constraints imposed by the ETM+ midscan correction process as described in the scan mirror calibration algorithm section.
Version 3.2
164
7/9/98
IAS Geometric ATBD Version 3
12
11
10
9
8
7
6
5
4
3
2
Roll
20
1
10
0
Scenes
Pitch Alignment Accuracy
60
50
40
30
12
11
10
9
8
7
6
5
4
3
2
20
Pitch
10
0
Scenes
Yaw Alignment Accuracy
60
50
40
30
12
11
10
9
8
7
6
5
4
3
20
Yaw
2
Yaw
50.9
36.0
29.4
25.5
22.8
20.8
19.2
18.0
17.0
16.1
15.3
14.7
30
1
Scenes
1
2
3
4
5
6
7
8
9
10
11
12
40
1
Pitch
50.9
36.0
29.4
25.5
22.8
20.8
19.2
18.0
17.0
16.1
15.3
14.7
50
Arc Seconds
Scenes
1
2
3
4
5
6
7
8
9
10
11
12
Roll Alignment Accuracy
Arc Seconds
Roll
49.9
35.3
28.8
25.0
22.3
20.4
18.9
17.6
16.6
15.8
15.0
14.4
Arc Seconds
Scenes
1
2
3
4
5
6
7
8
9
10
11
12
10
0
Scenes
Figure 3.1.6.2.1 1 Alignment Estimation Accuracy vs. Number of Scenes
In the unconstrained case it is necessary to process 40 scans (forward or reverse) to
achieve the 1.75 microradian accuracy target. This implies that the DOQ reference images
Version 3.2
165
7/9/98
IAS Geometric ATBD Version 3
must be a full scan wide (cross track) and at least 80 scans long (along track). Allowing
for some loss of data due to cloud cover and poor correlation targets leads to a
requirement for 90100 scans, or about onequarter of a WRS scene. Applying the scan
constraints simplifies the estimation problem by reducing the number of degrees of
freedom in the solution and, consequently, reduces the data requirements to 15 forward or
reverse scans. This can be satisfied by a DOQ reference image which is a full scan wide
and 30 scans long, or about oneeighth of a WRS scene.
Based on this analysis, we conclude that a single scan mirror calibration test site covering
at least onequarter of a WRS scene would be sufficient to provide usable results. The
goal is to provide at least two test sites each covering a full WRS scene, if possible. As
noted above, the input data accuracy assumptions are conservative, so using a full scene
reference image, a more accurate 1Gt product, and better imagetoimage correlation, the
performance of the algorithm with real Landsat 7 data is expected to exceed the 1.75
microradian design threshold.
3.1.6.2.3 Band Placement Calibration Accuracy
The IAS is required to determine band center field angles to an accuracy of 0.18 arc
seconds (1 σ per coordinate) per IAS Element Specification requirement number 3.2.3.10.
The bandtoband registration assessment algorithm uses test points which appear in all 8
ETM+ bands and performs crossband image correlation to measure any residual subpixel band to band offsets. These measurements are then used in the band placement
calibration algorithm to estimate band center offsets in the ETM+ focal planes, relative to
band 8 (panchromatic) which is treated as the reference band. These algorithms are
described above.
In the course of developing the bandtoband registration assessment algorithm it was
determined that, given a suitable test scene, it is possible to perform crossband image
correlation to an accuracy of approximately 0.1 pixels using preselected test points, and
to an accuracy of approximately 0.25 pixels using randomly generated test points. When
correlating bands of differing resolution (i.e., bands 6 or 8 to any of the others), the
measurements are performed at the coarser band resolution. Thus, the 0.1 or 0.25 pixel
measurement accuracy corresponds to a larger angular uncertainty for band offset
measurements for band pairs which include band 6. This makes band 6 the critical case
which drives the number of bandtoband test points required to meet the 0.18 arc second
accuracy specification.
An error propagation analysis was used to determine the number of test points required
to achieve a band center location accuracy of 0.18 arc seconds as a function of the
expected accuracy of the test point measurements. This analysis assumed that all test
points were measured in all band pairs and that the mapping between image pixel
Version 3.2
166
7/9/98
IAS Geometric ATBD Version 3
measurements and focal plane offsets could be approximately modeled as a constant scale
factor. The availability of measurements from seven band pair combinations makes the
band center location solution over determined so a least squares procedure was set up to
form the normal equation matrix which would be used to estimate band center offsets
from test point measurements. These normal equations were then inverted to compute
the predicted a posteriori covariance matrix for the band center estimates. The square
roots of the diagonal terms of this matrix were taken to be estimates of the band center
location standard deviations. The number and accuracy of the input test point
measurements were then varied to examine the sensitivity of the band center estimate
standard deviations to the input data volume and accuracy. As expected, the predicted
standard deviation for band 6 location was highest. Figure 3.1.6.2.31 shows a plot of the
number of test points needed to achieve a predicted standard deviation of 0.18 arc
seconds versus the accuracy of the input test point measurements.
As the figure shows, the accuracy requirement can be met with 2025 high quality hand
picked test points or with 125150 automatically generated test points. The accuracy of
hand picked points should be more reliable since the analyst selecting the test points can
ensure their quality and consistency. Based on the relatively small number of test points
needed (if they are sufficiently accurate), it was decided to use approximately 25 analyst
selected test points to perform the bandtoband registration assessment and subsequent
band placement calibration. The error propagation analysis predicts a band center location
accuracy of 0.09 arc seconds for bands 15,7 and 0.15 arc seconds for band 6. Band 8 is
the reference and is assumed to be known exactly. This meets the requirement for all
bands if the test point measurements are accurate to 0.12 pixels and for all but the thermal
band if the measurements are accurate to 0.19 pixels.
Version 3.2
167
7/9/98
IAS Geometric ATBD Version 3
Test Point
Accuracy
(pixels)
Test
Points
Needed
0.50
407
400
0.45
332
350
0.40
263
0.35
202
0.30
149
0.25
103
0.20
66
0.15
37
0.10
17
Band Placement Calibration
# Points
450
300
250
200
150
100
50
0
0.50
Figure 3.1.6.2.3 1
0.45
0.40
0.35
0.30
0.25
0.20
0.15
Accuracy
BandtoBand Test Point Number vs. Accuracy
3.1.6.3 Geometric Characterization Accuracy
Only two of the four geometric characterization algorithms implemented in the Landsat 7
IAS have quantifiable accuracy requirements. The geometric accuracy assessment is
primarily a qualitative evaluation of the internal distortions present in ETM+ images.
Geodetic accuracy assessment provides a verification of the Level 1Gs product geodetic
accuracy requirement analyzed above. Test point mensuration error was included in the
error budget for Level 1 processing geodetic accuracy to account for the contribution of
the assessment process. The bandtoband and imagetoimage registration accuracy
assessments have specific performance requirements (called out in the IAS Element
Specification) which are analyzed in the following sections.
3.1.6.3.1 BandToBand Registration Accuracy
The IAS is required to be capable of performing bandtoband registration to an accuracy
of 0.28 pixels (at the ground sample distance of the coarser of the two bands) per
coordinate at the 90% confidence level. This requirement is stated in the IAS Element
Specification as requirement number 3.2.3.6. The error budget for bandtoband
registration includes the high frequency jitter displacements which can occur in the time
interval over which the eight ETM+ bands sample the same ground location. Since the
Version 3.2
168
7/9/98
0.10
IAS Geometric ATBD Version 3
maximum time between band samples is approximately 2 milliseconds, which
corresponds to the jitter sensor sampling period, the bandtoband jitter error is
essentially the error in the difference between two jitter measurements. Other
components of the bandtoband registration error include band field angle errors,
processing model error, and test point mensuration error.
Table 3.1.6.3.11 shows the error components combined into an error budget for the
bandtoband registration assessment process. The table columns contain the following
information:
•
•
•
•
•
•
•
Column 1: identifies the data element error source;
Column 2: contains the total effect of the error source in units of 30 meter pixels
at the 90% confidence level (this is equal to the square root of two times the error
for one band for the errors which are independent for all bands);
Column 3: contains the contribution of the error source to a single band in units of
30 meter pixels at the 90% confidence interval;
Column 4: contains the original specified or estimated error source one sigma
accuracy;
Column 5: contains the units of the original error specification;
Column 6: provides a brief text explanation of the error estimate source;
Column 7: shows where the error source was allocated.
The field angle stability and test point mensuration contributions are relative errors which
are only counted once. High frequency jitter, field angle knowledge, and processing errors
apply independently to all bands and are therefore counted twice in the bandtoband
error budget.
Based on the analysis shown in Table 3.1.6.3.11 the 0.28 pixel (90%) requirement can be
met but there is little error budget margin. Manually selected test points must be used to
ensure the highest correlation accuracy since relaxing the correlation accuracy to 0.11
pixels exceeds the registration budget.
Version 3.2
169
7/9/98
IAS Geometric ATBD Version 3
Source
Both
Bands
0.079
30m pixels
(90%)
0.056
Spec
1σ
0.3
Units
Comments
Segment
arcsec
0.30 arcsec (1 σ)
Satellite
ETM+ Jitter
(high frequency)
0.035
0.025
0.134
arcsec
0.134 arcsec (1 σ) high freq
ETM+
Field Angle
Knowledge
0.048
0.034
0.18
arcsec
0.18 arcsec (1 σ) calibration
IAS
Field Angle
Stability
0.143
0.143
0.76
arcsec
0.76 arcsec (1 σ) stability
(0.2 pan)
Processing Error
0.140
0.099
1.8
meter
Model budget
IAS
Mensuration
0.164
0.164
0.1
pixel
Correlation accuracy
IAS
RSS Estimate
0.277
(90%)
RSS Margin
0.041
(90%)
Specification
0.280
(90%)
High Frequency
Jitter
ETM+
Table 3.1.6.3.1 1 BandtoBand Registration Error Budget
3.1.6.3.2 ImageToImage Registration Accuracy
The IAS has a requirement to perform imagetoimage registration to an accuracy of 0.4
pixels (at the multispectral band pixel ground sample distance of 30 meters) per
coordinate at the 90% confidence level. This requirement is stated in the IAS Element
Specification as requirement number 3.2.3.8. The error budget for imagetoimage
registration includes the high frequency components of the overall Landsat 7 ETM+
geometric error budget which model the internal geometric variability unique to each
image, the expected accuracy of the precision correction process which uses imagetoimage tie points to perform the image registration, and the expected mensuration accuracy
for the registration assessment test points.
The high frequency internal errors are a subset of the geodetic error budget components
presented in Table 3.6.1.12 and include the jitter, scan mirror variability, and
modeling/processing errors. The imagetoimage registration correction model accuracy is
required to be 3.6 meters (1 σ) according to the IAS Element Specification requirement
number 3.2.3.9. The same estimation error covariance propagation model used to analyze
the sensor alignment estimation accuracy was used to verify this requirement. In this
case, the expected accuracy of the PCD ephemeris and attitude data was combined with
20 ground control points assumed to be accurate to 10 meters (corresponding to control
from a 1:24,000scale map and 0.5 pan pixel identification/mensuration error) to yield
Version 3.2
170
7/9/98
IAS Geometric ATBD Version 3
predicted postfit ephemeris and attitude accuracy estimates. These estimates were then
propagated to ground coordinate accuracy for a range of scan angles and then rootsumsquared to establish an accuracy estimate for the precision correction in ground space.
This procedure resulted in an accuracy estimate of 3.5 meters (1 σ) for the precision
correction process. The specification value of 3.6 meters was used for the imagetoimage error budget.
The final component of the error budget is the imagetoimage test point correlation
accuracy. Testing with Landsat 5 data has shown that correlation accuracy of 0.1 pixel is
achievable for selected targets. Preselected test points will be used to evaluate the
imagetoimage registration accuracy.
Table 3.1.6.3.21 shows these components combined into an error budget for the imagetoimage registration assessment process. The table columns contain the following
information:
•
•
•
•
•
•
•
Column 1: identifies the data element error source;
Column 2: contains the total effect of the error source in units of 30 meter pixels
at the 90% confidence level (this is equal to the square root of two times the error
for one image for those errors which are contained in both the reference and object
images);
Column 3: contains the contribution of the error source to a single image in units
of 30 meter pixels at the 90% confidence interval;
Column 4: contains the original specified or estimated error source one sigma
accuracy;
Column 5: contains the units of the original error specification;
Column 6: provides a brief text explanation of the error estimate source;
Column 7: shows where the error source was allocated.
Note that the high frequency errors are independent for the two images and therefore
contribute twice whereas the precision correction and mensuration errors apply to the
process of relating the two images to each other and are therefore only counted once.
Version 3.2
171
7/9/98
IAS Geometric ATBD Version 3
Source
Both
Images
0.079
30m pixels
(90%)
0.056
Spec
1σ
0.3
Units
Comments
Segment
arcsec
0.30 arcsec (1 σ)
Satellite
Low Frequency
Jitter
0.079
0.056
0.3
arcsec
0.30 arcsec (1 σ)
Satellite
ETM+ Jitter
(high frequency)
0.035
0.025
0.134
arcsec
0.134 arcsec (1 σ) high freq
ETM+
Scan Mirror
Repeatability
0.191
0.135
0.72
arcsec
1.75 microradians * 2 (1 σ)
ETM+
Processing Error
0.140
0.099
1.8
meter
Model budget
IAS
Correction Model
Error
0.198
0.198
3.6
meter
Estimation error
IAS
Mensuration
0.164
0.164
0.1
pixel
Correlation accuracy
IAS
RSS Estimate
0.318
(90%)
RSS Margin
0.243
(90%)
Specification
0.400
(90%)
High Frequency
Jitter
Table 3.1.6.3.2 1 ImagetoImage Registration Error Budget
Based on the analysis shown in Table 3.1.6.3.21 the 0.4 pixel (90%) requirement can be
met but there is little error budget margin. Preselected test points must be used to ensure
the highest correlation accuracy since relaxing the correlation accuracy to 0.14 pixels
exceeds the registration budget. Similarly, it is necessary to use the precision correction
process to fit the object image directly to the reference image rather than fitting both to a
separate control source, leading to independent correction model errors which, when
combined, would exceed the 3.6 meter imagetoimage correction model budget.
Version 3.2
172
7/9/98
IAS Geometric ATBD Version 3
4. Constraints, Limitations, Assumptions
Atmospheric Refraction Correction
The atmosphere causes a ray of light to bend towards the satellite as it propagates
through the atmosphere. The amount of bending depends on the angle of propagation and
the density of the atmosphere at a given instant. For Landsat 7, the magnitude of this
affect is less than one meter at end of scan. Currently, this affect is not accounted for.
Resampling
Resampling Weight Tables (i.e. for cubic convolution, MTF, or other resampling kernels)
are calculated at 1/32 pixel increments. These weights assume evenly spaced pixels and
are used in the sample direction.
Spline weights (used for resampling over the scangap where the pixels are not necessarily
evenly spaced) are also calculated at 1/32 pixel increments. These weights are applied in
the line direction in the area between two scans. The weights are set up for a shift of 0.0
to the maximum gap possible from the system. A value of 3.0 covers the maximum gap
possible in the LANDSAT 5 satellite. A value of 6.0 will be needed for the LANDSAT
7 satellite, to cover the panchromatic band. When the panchromatic band is used for
LANDSAT 7 this could cause some anomalies in cases of large varying grey values
between pixels that are located in an area with a large scan gap.
Speed of Light Correction
Due to the noninfinite velocity of light, the velocity of a point on the surface of the Earth
and the velocity of the satellite, cause pixel location errors if not accounted for. The
speed of light correction due to the Earth's rotation is submeter and is currently not
accounted for.
Version 3.2
173
7/9/98