of 16
JournalofGeophysicalResearch:Planets
RESEARCHARTICLE
10.1002/2014JE004678
Key Points:
The location of Mars prime meridian
has been updated
New IAU expressions are needed
for the Mars spin state and
prime meridian
Unique processing techniques and
results are presented for THEMIS,
HRSC, and MOLA
Correspondence to:
T. C. Duxbury,
tduxbury@gmu.edu
Citation:
Duxbury, T. C., P. Christensen,
D. E. Smith, G. A. Neumann, R. L.
Kirk, M. A. Caplinger, A. A. Albee,
N. V. Seregina, G. Neukum, and
B. A. Archinal (2014), The location
of Airy-0, the Mars prime meridian
reference, from stereo photogrammet-
ric processing of THEMIS IR imaging
and digital elevation data,
J. Geo-
phys. Res. Planets
,
119
, 2471–2486,
doi:10.1002/2014JE004678.
Received 18 JUN 2014
Accepted 3 NOV 2014
Accepted article online 11 NOV 2014
Published online 5 DEC 2014
The location of Airy-0, the Mars prime meridian reference,
from stereo photogrammetric processing of THEMIS IR
imaging and digital elevation data
T. C. Duxbury
1
, P. Christensen
2
, D. E. Smith
3
, G. A. Neumann
4
, R. L. Kirk
5
, M. A. Caplinger
6
,
A. A. Albee
7
, N. V. Seregina
1,8
, G. Neukum
9
, and B. A. Archinal
5
1
School of Physics, Astronomy and Computational Sciences, George Mason University, Fairfax, Virginia, USA,
2
School of Earth and Space Exploration, Arizona State University, Tempe, Arizona, USA,
3
Department of Earth,
Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA,
4
NASA Goddard Space Flight Center, Greenbelt, Maryland, USA,
5
Astrogeology Science Center, U.S. Geological Survey,
Flagstaff, Arizona, USA,
6
Malin Space Science Systems, San Diego, California, USA,
7
Department of Geological and
Planetary Sciences, California Institute of Technology, Pasadena, California, USA,
8
Geology Faculty, Department of
Geocryology, Moscow State University, Moscow, Russia,
9
Department of Planetary Sciences and Remote Sensing,
Freie University, Berlin, Germany
Abstract
The small crater Airy-0 was selected from Mariner 9 images to be the reference for the Mars
prime meridian. Initial analyses in the year 2000 tied Viking Orbiter and Mars Orbiter Camera images of
Airy-0 to the evolving Mars Orbiter Laser Altimeter global digital terrain model to update the location of
Airy-0. Based upon this tie and radiometric tracking of landers/rovers from Earth, new expressions for the
Mars spin axis direction, spin rate, and prime meridian epoch value were produced to define the orientation
of the Martian surface in inertial space over time. Since the Mars Global Surveyor mission and Mars Orbiter
Laser Altimeter global digital terrain model were completed some time ago, a more exhaustive study has
been performed to determine the accuracy of the Airy-0 location and orientation of Mars at the standard
epoch. Thermal Emission Imaging System (THEMIS) IR image cubes of the Airy and Gale crater regions were
tied to the global terrain grid using precision stereo photogrammetric image processing techniques. The
Airy-0 location was determined to be about
0
.
001
east of its predicted location using the currently defined
International Astronomical Union (IAU) prime meridian location. Information on this new location and how
it was derived will be provided to the NASA Mars Exploration Program Geodesy and Cartography Working
Group for their assessment. This NASA group will make a recommendation to the IAU Working Group on
Cartographic Coordinates and Rotational Elements to update the expression for the Mars spin axis direction,
spin rate, and prime meridian location.
1. Introduction
Knowledge of the Mars orientation in inertial space is required when targeting landers to specific locations,
observing its surface from Earth and spacecraft, for tracking landers on the surface of Mars from orbiters
and directly from Earth, for determining the gravity field from orbiting spacecraft, etc. The orientation is
defined by the Mars spin axis direction, the location of the prime meridian, and the spin rate to propagate
the orientation over time. The Mars orientation originally had been determined from Earth-based telescopic
observation where the large crater Sinus Meridiani was adopted as the prime meridian reference. The loca-
tion of Mars surface features in inertial space when propagating ahead in time had an accuracy at the tens
of kilometers level due to errors in spin axis direction and spin rate. After the Mariner 6 and 7 flybys and
Mariner 9 orbiter explorations, the knowledge of the orientation of Mars in inertial space was improved to
the 10 km level.
de Vaucouleurs et al.
[1973] refined the prime meridian reference by using Mariner 9 images
to fix it to the small crater Airy-0 within the large crater Airy.
The spin axis direction and spin rate continued to be improved by the Earth-based tracking of the Viking
and Mars Pathfinder landers. Improvements to the location of the crater Airy-0 came after the orbiting Mars
Global Surveyor (MGS) [
Albee and Arvidson
, 1991] started producing a global digital terrain model (DTM)
from its Mars Orbiter Laser Altimeter (MOLA) [
Smith et al.
, 2001] in the late 1997. By the end of 2001, MOLA
had amassed over 600 million altimeter points, distributed globally with an average spacing of 300 m in
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2471
Journal of Geophysical Research: Planets
10.1002/2014JE004678
Figure 1.
In 2001, the Airy-0 location was determined by registering nine
Viking Orbiter images to an early MOLA DIM with features A–H being
MOLA control points.
latitude and 2 km spacing in longi-
tude at Mars equator. With each MGS
orbit crossing over thousands of other
orbits and since there were thousands
of orbits, millions of MOLA crossover
points were available for detailed pro-
cessing. The differences between all
individual observed and predicted
altimetry points where orbits crossed
are now at the meter level, obtained
by constraining the MGS spacecraft
orbits and pointing to fit mill ions
of orbit crossover altimetry points
together with a precision gravity field
[
Neumann et al.
, 2001]. However, tying
imaging data to the individual MOLA
points is limited at the 100 m level
due to the larger spacing of the MOLA
points in latitude and longitude. The
MOLA global DTM was the most sig-
nificant development for precision
Mars cartography since the beginning
of Mars exploration. A detailed history
of the Mars prime meridian is given in
Archinal and Caplinger
[2002].
In 2000, the International Astronomical Union (IAU) Working Group on Cartographic Coordinates and
Rotational Elements (WGCCRE) was preparing its publication that would include the spin axis, spin rate,
and prime meridian expressions for Mars. Even though the MOLA data coverage was far from complete
and the precision cartographic processing had only begun, the MOLA data set was sufficient to register a
Viking Orbiter-controlled photomosaic to a MOLA-derived digital image model (DIM) by one of the authors
(Duxbury) in 2001 (Figure 1). Also, a single Mars Orbiter Camera [
Malin et al.
, 1992] image was registered to a
single MOLA ground track by two of the authors (Caplinger and Neumann) in 2001. These two efforts gave
solutions of the areocentric coordinates of Airy-0 that agreed within a few hundred meters. Those were aver-
aged to produce the location of Airy-0 relative to the MOLA grid that was adopted within the NASA Mars
Exploration Program Office by its Geodesy and Cartography Working Group (MGCWG) [
Duxbury et al.
, 2002].
Airy-0 was found to be 0.25
or 15 km from the IAU reference position. The fundamental epoch angle (W
0
)of
the IAU prime meridian expression was adjusted as were the longitudes of all MOLA points, so when images
of Airy-0 were registered to MOLA, Airy-0 would have a computed longitude of 0
. The spin axis direction,
spin rate, and prime meridian expressions based upon this Airy-0 location and the spin axis/rate work by
Yuan et al.
[2001] are referred to as the “IAU 2000 pole and prime meridian of Mars” [
Seidelmann et al.
, 2002].
Since this WGCCRE publication, no work had been performed to improve the location of Airy-0 until the
MGCWG recently initiated studies of the Mars orientation to support a new IAU WGCCRE publication to
update the Mars orientation expressions. Based upon this revived interest, new results of the Airy-0 loca-
tion are presented in this article. Stereo photogrammetric data processing was used to register multispectral
infrared image cubes to MOLA and High Resolution Stereo Camera (HRSC) [
Neukum et al.
, 2009] DTMs. The
2001 Mars Odyssey (ODY) [
Saunders et al.
, 2004] Thermal Emission Imaging System (THEMIS) IR [
Christensen
et al.
, 2004] image cubes of the Airy crater and Gale crater comprised the image data set. The precision
global MOLA DTM, the local HRSC DTM, and THEMIS IR image cubes did not exist in 2001; therefore, there is
significant independence of the results presented here to the earlier 2001 results.
The spatial resolution of the THEMIS image cubes is 100 m/pixel. MOLA DTMs surrounding Airy and Gale
craters were constructed from all of the individual final/corrected MOLA ground tracks [
Neumann et al.
,
2001] within these areas and illuminated under similar conditions as each of the THEMIS IR image cubes
to produce DIMs at 230 m/pixel. Additionally, a DTM and DIM of Gale crater were produced by the German
Space Center (DLR) [
Gwinner et al.
, 2010] that were derived by registering HRSC stereo images to MOLA with
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2472
Journal of Geophysical Research: Planets
10.1002/2014JE004678
Figure 2.
THEMIS IR microbolometer area array layout with 10 spectral
bands looking down on the array.
a spatial resolution of 50 m/pixel
(http://europlanet.dlr.de/node/index.
php?id=380GaleDTM). The Airy crater
data reduction gave the direct loca-
tion of Airy-0, while the Gale crater
data reduction gave an independent
verification of the THEMIS IR geo-
metric parameter values versus those
derived with the Airy-0 location and
the stability of the THEMIS IR param-
eter values over time and over the
Martian surface.
By-products of this analysis
included precision 10-band THEMIS
IR-controlled photomosaics of
the Airy and Gale crater regions
with thousands of precision con-
trol and tie points. The controlled
photomosaic of Gale crater was
used by the MGCWG to provide an independent validation of the map products used by the Mars Sur-
face Laboratory (MSL) [
Grotzinger et al.
, 2012] flight operations team to plan and target its successful
landing. Precision THEMIS IR geometric and pointing calibrations, archived as NAIF SPICE [
Acton
, 1996;
Semenov et al.
, 2005] Instrument (IK) and Frames (FK) kernels, were produced as part of the stereo pho-
togrammetric reduction. These kernels were used in producing the THEMIS IR 100 m/pixel, 10-band
global Mars digital image mosaic (MDIM) [
Edwards et al.
, 2011]. Additionally, the precision and noise
characteristics of the reconstructed ODY spacecraft attitude were determined over an 8 year time period.
The following sections describe the cartographic mapping properties of THEMIS IR, how local DTM control
point and THEMIS tie point networks were established in the Airy and Gale crater regions, how the THEMIS
IR image cubes were measured to produce the tens of thousands of linear equations used in the stereo
photogrammetric reduction process, how the THEMIS IR images were registered to the MOLA/HRSC refer-
ence surfaces for absolute control, how the observations were validated and weighted, the results of the
data reduction in determining the location of Airy-0, the THEMIS IR geometric and alignment parameter
values, the precision/noise characteristics of the observations, the control and tie point networks, and the
spacecraft attitude data over the two crater sites and over 8 years of time.
Figure 3.
THEMIS IR active pixel area.
2. Cartographic Mapping Properties
of THEMIS IR
THEMIS produces both multispectral infrared
(IR) push broom and multispectral visual (VIS)
area array images. To date, these images have
been mostly taken along the spacecraft ground
track, but recently, some images have been taken
off nadir. The THEMIS IR uses a 320
×
240 pixel
microbolometer area array with different narrow-
band spectral filters covering 10 sections of the
detector (Figure 2). Each filter covers 16 lines that
are read out at a rate designed to match the ground
track speed from the near-circular and near-polar
ODY orbit. Line summing of the 16 lines within
each spectral band, using time delay integration,
produces one effective image line for each band
to increase the signal-to-noise ratio. Therefore, in
normal IR imaging mode, 10 equivalent lines,
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2473
Journal of Geophysical Research: Planets
10.1002/2014JE004678
Figure 4.
The NASA Mars Odyssey spacecraft-fixed
̄
X
̄
Y
̄
Z
coordinate system.
covering narrow IR spectral bands, form the push
broom, multispectral image cubes.
A THEMIS IR pixel within the microbolometer area
arrayis50
×
50
μ
m
2
, but only a 45.4
×
24.6
μ
m
2
area is active (Figure 3) to image Mars, thus limit-
ing image metric measurement residuals to the few
tens of a pixel level. This occurs because each of the
16 pixels summed to produce one effective pixel per
spectral band and the effective pixels in the other
nine different bands have each seen a slightly dif-
ferent view of the surface due to spacecraft orbital,
attitude, and Mars rotational motion during image
construction. An automated digital image match-
ing/correlation technique cannot achieve greater
accuracy because the pixels will not match exactly,
even if they have the same spectral band pass.
Each of the 10 IR bands share optical distortion and
spacecraft mounting alignment/pointing parame-
ters but have slightly different values for time tag,
pointing, spacecraft position, and focal length.
In order to relate individual THEMIS IR image pixels to Mars-fixed surface coordinates, multiple reference
systems are used that include the following.
2.1. Mars-Fixed System
The Mars Mean Equator and prime meridian of-date reference system has axes
̄
X
m
̄
Y
m
̄
Z
m
where
̄
X
m
and
̄
Y
m
define the Mars equator and
̄
Z
m
defines the Mars pole (spin axis): (1)
̄
X
m
is in the of-date Mars equator and
along the prime meridian tied to the crater Airy-0; (2)
̄
Z
m
is along the of-date Mars north pole; and (3)
̄
Y
m
is
in the of-date Mars equator, completing the orthogonal, right-handed system.
The orientation of these axes in inertial space at any Barycentric Dynamical Time, TDB, referred to as
“of-date,” is recommended by the IAU WGCCRE [
Archinal et al.
, 2011] based on
Duxbury et al.
[2002] and are
included in the NAIF SPICE Planetary Constants kernel f
tp://naif.jpl.nasa.gov/pub/naif/M01/kernels/pck/
pck00010.tpc.
A Mars-centered and fixed surface position vector
̄
u
having areocentric coordinates of latitude,
positive
north, longitude,
positive east, and radius,
u
in kilometers, is given by
Figure 5.
THEMIS IR microbolometer area array image plane coordi-
nates looking down on the array.
̄
u
=
u
cos
cos
cos
sin
sin
=
u
x
u
y
u
z
(km) (1)
where
u
x
,
u
y
, and
u
z
are the position
components of
̄
u
along the Mars-fixed
axes
̄
X
m
̄
Y
m
̄
Z
m
.
The Mars-centered and fixed posi-
tion vector
̄
u
is translated to the
spacecraft-centered and Mars-fixed
position vector
̄
v
using
̄
v
=
̄
u
̄
w
(km)
(2)
where
̄
w
is the Mars-centered and
fixed spacecraft position vector at TDB
obtained from a NAIF SPICE spacecraft
ephemeris kernel ftp://naif.jpl.nasa.gov/
pub/naif/M01/kernels/spk/*.bsp.
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2474
Journal of Geophysical Research: Planets
10.1002/2014JE004678
Figure 6.
Simulated THEMIS IR image in Marwth by
illuminating a MOLA-derived DTM to produce a DIM
under the same lighting conditions as the THEMIS
IR cube. (left) MOLA DIM, (middle) IR band 5 image,
and (right) simulated IR image from MOLA DIM.
The .bsp files having “mgs095j” in the file name have abso-
lute position accuracies at the meter level [
Konopliv et al.
,
2011] and were used in this analysis. At this level of accu-
racy, the ODY positions at the times of the image cubes
were not estimated but were held as fixed parameters.
2.2. The 2001 Mars Odyssey Spacecraft System
The ODY spacecraft-fixed coordinate system has axes
̄
X
̄
Y
̄
Z
(Figure 4) where during normal orbital operations,
̄
Z
is
nadir pointed toward Mars;
̄
X
is canted 17
from the orbit
plane, orthogonal to
̄
Z
; and
̄
Y
completes the orthogonal,
right-handed system.
The 3
×
3 unitless transformation matrix
C
defines the
orientation of the Odyssey spacecraft in Mars-fixed coordi-
nates and is reconstructed from telemetered quaternions
based upon star camera and inertial measurement unit
observations. The
C
matrix allows one to transform from
the of-date Mars body-fixed system to the Odyssey
spacecraft-fixed coordinate system using:
̄
X
̄
Y
̄
Z
=
C
̄
X
m
̄
Y
m
̄
Z
m
(3)
where
C
is obtained from the 2001 Mars Odyssey NAIF
SPICE archive ftp://naif.jpl.nasa.gov/pub/naif/M01/
kernels/ck/*.bc.
Given a
C
kernel [
Acton
, 1996], the spacecraft-centered
and Mars-fixed position vector
̄
v
is transformed to a
spacecraft centered and fixed position vector
̄
r
by using
̄
r
=
C
̄
v
(km)
(4)
For 2001 Mars Odyssey, the noise on the
C
kernels was found to be about 0.01
(3
).
2.3. THEMIS IR-Fixed System
The THEMIS IR base/optics (
TO
) are hard mounted to the spacecraft. The alignment of this base system rel-
ative to the ODY spacecraft axes was determined prelaunch and will be held fixed for this analyses. The
transformation
T
ODY
TO
from the ODY-fixed
̄
X
̄
Y
̄
Z
to the THEMIS base/optics system
TO
is defined by
T
ODY
TO
=
[
0
.
02378
o
]
1
[
0
.
06119
]
2
[
16
.
91498
]
3
(5)
where the
≈−
17
rotation about the third axis removes the spacecraft cant angle to align the THEMIS IR
lines normal to the orbit plane. The form
[
]
i
for
i
= 1, 2, or 3 represents a unitless 3
×
3 rotation matrix
given by
[
]
1
=
10 0
0cos
sin
0
sin
cos
,
[
]
2
=
cos
0
sin
010
sin
0cos
,
and
[
]
3
=
cos
sin
0
sin
cos
0
001
.
(6)
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2475
Journal of Geophysical Research: Planets
10.1002/2014JE004678
Figure 7.
(left) MOLA- and (right) HRSC-derived DIMs.
From this base/optics system, one can then transform to the THEMIS IR-fixed microbolometer area array
coordinate system that has axes
̄
x
̄
y
̄
z
(Figures 2 and 5) where
̄
x
is in the direction of increasing IR sample
numbers and normal to the orbit plane;
̄
z
is along the nadir direction during nadir operations; and
̄
y
com-
pletes the orthogonal, right-handed system, is along the decreasing line direction, and is in the general
direction of the spacecraft velocity vector.
The following 3
×
3 unitless transformation matrix
T
TO
IR
from THEMIS base/optics to THEMIS IR fixed
̄
x
̄
y
̄
z
coordinates is expressed as
T
TO
IR
=
[
w
]
3
[
y
]
2
[
z
]
3
(7)
From ground calibrations,
z
= 0.0
,
y
=90
, and
w
= 0.672
.
As previously mentioned, the angles in equation (5) will be held fixed but the angles in equation (7) will be
solved for in the stereo photogrammetric reduction process where there will be three estimated angles for
each 10-band IR image cube. The averages of these three angles over all IR cubes will be referred to as the
fixed mounting alignment angles, and the differences between the fixed biases and the individual values for
each cube will be used to asses the spacecraft attitude pointing accuracy and stability.
Figure 8.
An IR image was registered to the HRSC DIM of the same feature to measure control point locations. Image
misregistration produces red/blue-green ghost images as seen on the left.
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2476
Journal of Geophysical Research: Planets
10.1002/2014JE004678
Figure 9.
Measured and predicted image locations were displayed as
plus and cross overlays, respectively, to validate measurements in all
bands and all cubes. Usually, the measured and predicted overlays were
on top of each other so only the plus was seen as is the case here where
the stereo reduction process was completed and all bad image location
measurements were corrected or removed. These residuals, the differ-
ence between the measured and predicted image locations, verified that
the measurement accuracy was at a subpixel level. Note that subpixel
measurement accuracy was achieved even with the significantly different
appearance of the features due to the thermal properties in the different
spectral bands and different viewing/lighting conditions.
Equations (5) and (7) can be concate-
nated to give the transformation
T
ODY
IR
from ODY-fixed to THEMIS IR-fixed
coordinates by
T
ODY
IR
=
T
TO
IR
T
ODY
TO
(8)
The spacecraft centered and fixed posi-
tion vector
̄
r
in equation (4) can then
be transformed to a THEMIS centered
and fixed position vector
̄
p
by
̄
p
=
T
ODY
IR
̄
r
=
p
x
p
y
p
z
(km) (9)
where
p
x
,
p
y
, and
p
z
are the compo-
nents of position vector
̄
p
in THEMIS
IR-centered and fixed
̄
x
̄
y
̄
z
coordinates.
2.4. THEMIS IR Focal Plane/
Image Coordinates
There are multiple methods to map
from the THEMIS-fixed position vec-
tor
̄
p
through the optics to the focal
plane and then to THEMIS IR image coordinates. In mapping through the optics to the image plane, a square
shape in object space becomes a “keystone” shape in the image. THEMIS uses a beam splitter to focus some
light onto the IR detector with the rest of the light onto the VIS detector. One explanation for the keystone
effect is that the beam splitter is not exactly aligned 45
going to the IR detector. The net effect, from what-
ever the cause, is that each of the 10 bands has a different scale along the 320 samples within a line. The first
model used in this analysis to map from
̄
p
to image coordinates will give each of the 10 bands its own focal
length along the 320 samples.
Selecting a scale across the 10 bands also has its problems. Since each band has its own focal length along
each line, no single focal length can be used to represent the scale across all 10 lines and obtain subpixel
model accuracy. However, a single scale in the line direction was implemented in the first model, and then a
bias was added to each band locations to account for the differences between the single line scale and the
actual band locations. This model provides parameters that are easily visualized and gives quick insight into
the THEMIS IR geometric properties.
The nominal image plane
x
,
y
mm coordinates (Figures 2 and 5) on the THEMIS IR area array detector, asso-
ciated with
̄
p
, are computed using the colinear equation of photogrammetry [
Thompson
, 1966] where the
camera optics are modeled as a pinhole
[
x
y
]
=
1
p
z
[
f
i
p
x
fp
y
]
(nominal, mm)
(10)
where
f
i
is the THEMIS IR effective focal length along the sample direction
̄
x
for each of the IR bands
(i=1,...,10) and
f
is an average focal length along the line direction
̄
y
, across all of the IR bands.
If the optics produce radial distortion, the actual image plane coordinates
x
and
y
are shifted relative to the
nominal
x
,
y
mm coordinates and are computed from
[
x
y
]
=
[
x
y
]
+
r
2
[
x
y
]
1
+
r
4
[
x
y
]
2
(distorted, mm)
(11)
where
r
2
=
x
2
+
y
2
mm
2
(12)
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2477
Journal of Geophysical Research: Planets
10.1002/2014JE004678
and where
1
and
2
are the coefficients of radial optical distortion. An optical distortion term propor-
tional to
r
would not be distinguishable from a change in focal length and is not included in the model.
Equation (11) assumes that the optical principal point is at the center of the THEMIS IR area array detector,
an approximation that this analysis found was sufficient to maintain subpixel model accuracy.
The predicted THEMIS IR sample
s
and line
l
image coordinates, associated with the position vector
̄
p
,are
computed from
[
s
l
]
=
[
s
0
l
0
]
+
K
[
x
y
]
+
[
0
b
i
]
(raw image coordinates)
(13)
For THEMIS IR,
s
0
= 160.5,
l
0
= 120.5, the IR area array center coordinates, and
K
= 20.000 pixels/mm since a
THEMIS IR pixel is 0.050 mm across. These three CCD parameter values are fixed and are not estimated. The
minus sign is needed in equation (13) for the
y
coordinate since
̄
y
is in the direction of decreasing line num-
ber. The
b
i
terms absorb errors in approximating the scale in the line direction
̄
y
with a single focal length
f
and are a fraction of an IR pixel.
Because of systematic residual errors when using the first model, a second model was additionally evaluated
to map from position vector
̄
p
to THEMIS IR image coordinates
s
and
l
. Given the THEMIS IR centered and
fixed position vector
̄
p
to a surface feature, one can compute the IR focal plane coordinates
x
and
y
[
x
y
]
=
f
k
p
z
[
p
x
p
y
]
(nominal, mm)
(14)
where
f
k
is the average focal length across the entire IR focal plane. THEMIS IR image coordinates
s
and
l
are
then computed directly from
[
s
l
]
=
[
s
0
l
0
]
+
K
[
K
11
K
12
K
13
K
14
K
15
K
21
K
22
K
23
K
24
K
25
]
x
y
xy
x
2
y
2
(15)
For the case of no optical or keystone distortions,
K
11
=
1
,
K
22
=−
1
, the rest of the
K
terms have values
of zero, and
K
,
s
0
, and
l
0
have the same values as in equation (13). The values of
f
k
and
K
ij
are
estimated parameters.
As mentioned previously, the IR area array detector is divided into 10 areas under different narrowband IR
filters where 16 lines within these bands are summed to produce better signals, yielding one equivalent line
observation per band (Figure 2). The center line number of each band is adopted here as the observed line
coordinate of any feature (e.g., a control point) within that band. Therefore, the observed image coordinate
of a feature can have a sample value
s
ranging between 1 and 320, but the observed line coordinate
l
in
equations (13) and (15) is fixed for each of the 10 IR bands that have the following values:
l
1
=
8
.
5
for IR band #1
l
2
=
24
.
5
for IR band #2
l
3
=
50
.
5
for IR band #3
l
4
=
76
.
5
for IR band #4
l
5
=
102
.
5
for IR band #5
l
6
=
128
.
5
for IR band #6
l
7
=
154
.
5
for IR band #7
l
8
=
180
.
5
for IR band #8
l
9
=
205
.
5
for IR band #9
l
10
=
231
.
5
for IR band #10
(16)
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2478
Journal of Geophysical Research: Planets
10.1002/2014JE004678
Figure 10.
A seamless controlled THEMIS IR photomosaic of the Martian
Airy crater area showing all control and tie points. The more dense popu-
lation of control and tie points in the north-south direction are in areas of
IR cube overlap. Tie point MED near the center and identified by the circled
plus is Airy-0.
3. Time Tag
When one is measuring an IR
image cube file (.QUB), one mea-
sures the record number
L
in the
file of a feature location. The value
of
L
is not used as the line obser-
vation associated with equation
(13) or (15). The line observation
is always one of the 10 values of
l
i
listed in equation (16). The value
of
L
is used to determine the time
tag of the
s
i
,
l
i
observation.
The readout rate within the IR
area array is constant and is not
affected by
f
i
,
f
,
b
i
,훼
1
,훼
2
,
f
k
,
K
ij
and
is
30 lines/s. The start time of
an IR image, defined in the QUB
image file label record, gives the
time of the first line in the entire
area array detector. The time tag
of any sample
s
in one of the 10
equivalent lines must include
the time needed to read out the
detector down to
l
i
, i=1,...,10.
Given an observed
s
i
,
l
i
image loca-
tion of a control point in band
i
,
i
=
1
, ...,
10
at QUB image record
L
, the timing offset
Δ
t
from the start time given in the QUB label record
for the observed image coordinates is computed from
Δ
t
=
l
i
+
L
1
RATE
s
(17)
For
L
= 1 of band
i
= 1, the time offset for this band is
0.55 s from the image start time defined in the QUB
label record.
The start time of an image cube could have an uncertainty of up to four IR image lines. An error in the cube
start time and an error in angle
y
are similar, only effecting the line direction, and cannot be separated in
the stereo reduction process. Therefore, an error in
Δ
t
for each IR image cube is not estimated, so its effect
will be included in the estimated value for
y
.
An example of using equations (1)–(5) is shown in Figure 6. A MOLA DTM in the Marwth region of Mars was
illuminated under the same lighting conditions as a THEMIS IR cube (middle) to produce the map gridded
DIM (left). Going backward through equations (1)–(15) and using the associated NAIF SPICE kernels, band
5 of the IR cube was simulated (right). One could then vary the modeled surface properties in an attempt
to bring the simulated image into agreement with the
actual image brightness to infer surface scatter-
ing/thermal properties (not included in the scope of this analyses). In producing this simulated THEMIS IR
image that was registered to the MOLA DIM, a readout rate of 30.0147 lines/s was determined, consistent
with the prelaunch calibrated value of 30.0477 lines/s. Over a typical 15 s THEMIS IR image, the difference
in readout rates maps to less than 0.5 lines. This in-flight determined value was held constant in the stereo
photogrammetric reduction processing across all THEMIS IR cubes.
4. Stereo Photogrammetric Linear Equations
Equations (1)–(15) allow one to map a location on the surface of Mars to a sample/line image location in any
band of a THEMIS IR cube given the START_TIME from the image label record and the associated NAIF SPICE
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2479
Journal of Geophysical Research: Planets
10.1002/2014JE004678
Figure 11.
The difference is shown between the two THEMIS IR
geometric models. The systematic differences between these two
models are due to the lack of freedom in the model 1 parameters that
model 2 provides.
kernels. One can also use these
equations to go backward from image
location
s
,
l
to Mars-fixed surface coordi-
nates, given a local Mars radius or digital
terrain model of the area imaged.
For the stereo photogrammetric
analyses of THEMIS IR cubes performed
here, the first-order partial
derivatives of
s
,
l
with respect to
f
i
,
f
,
b
i
,훼
1
,훼
2
,
f
k
,
K
ij
,휃
z
,휃
y
,휃
w
,휙,휆,
u
are used to construct linear equations
where the measured image locations of
surface features in IR cubes are used as
observables in a square root, weighted
least squares estimation process
[
Bierman
, 1977] to determine the
THEMIS IR geometric, alignment, and
pointing parameters and the control/tie
point locations on the surface of Mars,
including Airy-0.
MOLA DTMs and DIMs and HRSC-derived DTMs and DIMs (precisely tied to MOLA) provide “absolute con-
trol” in registering image cubes to the Mars reference surface. The term absolute control is used to reflect
that the millions of MOLA points have been tied together globally at the meter level but is limited by their
300 m spacing in latitude between points and a few kilometer spacing in orbital longitude at Mars equator.
Therefore, the DTMs and DIMs provide similar absolute information as a star catalog does for astromet-
ric reductions. Finally, when all camera parameters and surface feature locations are determined, one can
produce precision (subpixel accuracy) multiband, controlled photomosaics using equations (1)–(15) by
running through a range of
and
at a fixed spatial resolution using a standard map projection and pop-
ulating the mosaic with the appropriate IR image pixel brightness values. The MOLA DTM was used as the
reference surface to produce the Airy crater-controlled photomosaic, and the HRSC DTM was used as the
reference surface to produce the Gale crater-controlled photomosaic. The map-projected THEMIS IR image
cubes and the bands within each cube were registered seamlessly to better than 0.3 pixels (1
), the tie point
measurement accuracy.
Figure 12.
The solved for and derived THEMIS IR equivalent focal lengths
are displayed as a function of band and CCD location. The focal lengths
indicated by plus were derived from the
K
term values using equations (20)
and (21).
5. THEMIS IR Control and Tie
Point Measurements
Photomosaics were produced in sinu-
soidal map projections having spatial
resolutions of 100 m/pixel for both
Airy and Gale craters using THEMIS
IR band 5 images, a priori THEMIS
IR parameter values, the associated
SPICE kernels, and equations (1)–(15).
Misalignment in the areas of image
cube overlap was obvious because of
errors in the image projections due
to errors in the THEMIS IR parameter
values, random noise in the THEMIS
pointing (
C
kernel), and errors in
start times of the image cubes. A
low-resolution version of a mosaic
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2480
Journal of Geophysical Research: Planets
10.1002/2014JE004678
was displayed on a large cinema monitor that allowed the user to poke on any position within the mosaic
and display all THEMIS IR cubes/bands that covered the point. Because of the large number of overlapping
cubes being measured, up to six THEMIS IR cubes with all 10 bands within each cube could cover a point.
To easily and quickly transform between photomosaic
휙, 휆
space and THEMIS IR cube
s
,
l
space, equations
(1)–(15) were evaluated 300 times, evenly spaced between samples 1 and 320 and between all of the IR
image lines that were within the map projection limits for each band of all cubes. These sets of 300
휙, 휆,
s
,
l
groups for each band of all cubes were used to determine the coefficients in the following approximate
transformations, assuming a constant Mars radius within the photomosaic:
[
s
l
]
=
[
s
0
l
0
]
+
[
s
1
s
2
s
3
s
4
l
1
l
2
l
3
l
4
]
2
2
(18)
[
]
=
[
0
0
]
+
[
1
2
3
4
1
2
3
4
]
s
l
s
2
l
2
(19)
A priori values for equations (18) and (19) had initial accuracies at the 5–10 pixel level, sufficient to start the
iteration process. The coefficients in equations (18)
and (19) were continually updated as the evolving data
processing produced more accurate parameter values
, reducing the errors in the approximate transforma-
tion to a few pixel level. Once a location
휙, 휆
was selected in a photomosaic, equation (18) would be used
to compute which bands within all cubes contained this point to begin the image location measurement
process. Hundreds to thousands of locations throughout the photomosaic were selected, with concentra-
tions in the cube overlap areas to yield accurate THEMIS IR geometric and pointing parameters needed for
precision map projection.
Figure 7 shows two Gale crater DIMs derived from the MOLA and HRSC DTMs where the higher spatial res-
olution HRSC DIM is on the right. Features measured in both the IR image cubes and the DIMs are called
control points because their absolute locations (
휙, 휆,
u
) were determined precisely from the DIMs/DTMs.
Their Mars-fixed positions were held as constants throughout the stereo photogrammetric reduction pro-
cess, similarly as star positions are held as fixed/known parameters in an astrometric reduction process. To
provide complete and redundant coverage over the 320 samples and thousands of lines in each band of all
image cubes, and to precisely tie all bands together within a cube and between cubes, thousands of addi-
tional features only seen in the image cubes were measured and are referred to “tie” points. The tie point
Mars-fixed positions were not known and were estimated in the stereo photogrammetric process. Airy-0
was one of these tie points. Selecting and measuring tie points are described later.
Control points were measured by producing a color image with the red channel being filled with the MOLA
or HRSC DIM and the green and blue channels being filled with band 5 of an IR image (Figure 8). When the
two data sets are not registered, the MOLA DIM appears as a red ghost image of the surface. By shifting the
IR image in the blue and green channels with respect to the DIM in the red channel, the image becomes
more grey as the two data sets become registered. The residual colors on the right after registration repre-
sent spatial resolution differences and actual surface thermal variations from the simple DTM-derived DIMs.
Once registered, the IR
s
,
l
image coordinates become observables in the stereo photogrammetric process.
The DIM image coordinates are used to compute the value of
휙, 휆
for the control points from the sinusoidal
map projection equations relating mosaic image location to Mars coordinates. These same DIM coordinates
are used to extract the value of
u
from the associated DTM. For a tie point that was not identified in a DIM,
usually because it was small relative to the DIM spatial resolution, its approximate DIM location was used to
determine the approximate (a priori) values of
휙, 휆,
u
of the point. Then equation (18) was used to determine
which other bands in all image cubes contained the control or tie point.
Once a control or tie point was identified, a 21
×
21 pixel
2
area surrounding the point was extracted from
band 5 of a selected cube and used to make a 21
×
21 weighted digital Wiener filter [
Wiener
, 1949] to con-
volve with all bands in all cubes that contain the feature. The maximum output of the filter convolution
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2481
Journal of Geophysical Research: Planets
10.1002/2014JE004678
process was used to automatically measure the fea
ture location, using interpolation to obtain locations to
sub-THEMIS IR pixel resolution, and also used to compute a signal-to-noise ratio for the measured image
location. Figure 9 shows a THEMIS tie point that was attempted to be measured in all bands of six differ-
ent cubes. Only 46
s
,
l
observation pairs were extracted as the control point went out of the field of view
of some of the bands or the noise on the measured image location was too high to include in the stereo
photogrammetric analyses.
Ideally, the convolution filter should be changed to reflect the differences in the spectral bands within an
image cube and for the different lighting and viewing conditions between cubes. This was not done, simpli-
fied by the fact that the image cubes were all taken in nadir orientation and within a few hours of the same
time of day. Even with the images changing between the spectral bands because of the surface thermal
properties, the difference between the observed and predicted image locations was at the subpixel level,
and therefore, no further work was needed to account for surface thermal properties.
The elimination of bad measurements had to be implemented in stages in order to not eliminate accurate
measurements. In the early processing with a much smaller data set, measurements could appear to be bad
when in reality the model parameters were not accurate. Therefore, as more and more measurements were
added to the process and the parameter estimates became mor
e accurate, tighter constraints were placed
on defining bad measurements. The first step was to identify all potential measurement errors using the size
of the postfit residuals, the differences between the measured and the computed image locations. Points
with residuals larger than two pixels were remeasured if the measurement error was obvious; otherwise,
the measurements were removed. Usually, the measurements that were removed had high image noise
or the measured sample coordinates were outside of the 1–320 sample range. This first step removed all
measurement blunders at the few pixel level.
To identify remaining bad measurement that were harder to identify by just looking at residuals, the point-
ing for each image cube was reestimated by removing one control or tie point at a time. The postfit residual
measurement statistics of an image cube would improve by tens of percent or more when a single bad mea-
surement was removed. These bad points were remeasured or deleted as before. This process was repeated
by removing image location measurements from combinations of two points at a time in a cube since mul-
tiple measurement errors could have been made in a single cube. This process could identify when only one
or two bands of a single cube were measured in error. After this process, the postfit residuals over the entire
control and tie point data set were better than the one pixel level (1
).
The final technique was to display the measured and predicted image locations as overlays on the images
(Figure 9). The observed location is indicated by a “plus” and the predicted location as a “cross.” Displaying
all measurements in all bands and cubes together provided a powerful visual tool to identify measurement
location discrepancies relative to the feature within the individual images. This display was also very use-
ful in identifying measurement errors in noisy image bands or when the point was near the ends of a line.
Once these measurements were corrected or removed, the differences between the measured and pre-
dicted image locations, referred to as residuals and measurement accuracies, were at the 1.5 pixel level for
the MOLA control points and better than 0.3 pixels for the THEMIS tie points.
There were 89 MOLA control points and 501 THEMIS tie points in the Airy region, and there were 142 HRSC
control points and 380 THEMIS tie points for Gale crater.
One of the tie points was the crater Airy-0. There
were 3942
s
,
l
observations of the MOLA control points and 23,022
s
,
l
observations of the THEMIS tie points
to estimate 1562 tie point location and THEMIS IR parameters in the Airy-0 region. There were 8780
s
,
l
observations of the HRSC control points and 21,222
s
,
l
observations of the THEMIS tie points to estimate
1227 tie point location and THEMIS IR parameters in the Gale cr
ater region. Parameter estimates were pro-
duced for the Airy-0 data set only, the Gale crater data set only, and the combined Airy-0 and Gale crater
data sets. The solutions were iterated 4 times in the least squares parameter estimation process for all three
cases. No change to the postfit residual statistics were seen after the third iterations, and a forth iteration
was performed to verify this.
6. Results
6.1. Airy-0 Location
Four solutions were obtained for the location of Airy-0 from the Airy only and Airy-Gale combined data
sets using both THEMIS IR geometric models 1 and 2. The four solutions varied by only 3 m on the Martian
DUXBURY ET AL.
©2014. American Geophysical Union. All Rights Reserved.
2482