High-resolution digital elevation models derived ... - Pages perso de

1 downloads 0 Views 3MB Size Report
Ius-Tithonium Chasma Cook et al., 1999] has been produced by the Extraterrestrial Orbital DEMs for. Understanding Surfaces (EXODUS) project at Uni-.
High-resolution digital elevation models derived from Viking Orbiter images: Method and comparison with Mars Orbiter Laser Altimeter Data

1

D. Baratoux, C. Delacourt and, P. Allemand

Laboratoire des Sciences de la Terre, Universite Lyon I, Centre National de la Recherche Scienti que et Ecole Normale Superieure de Lyon, Villeurbanne, France

Abstract

The purpose of this paper is to present a method to derive local and high-resolution digital elevation models (DEM) of Mars from Viking Orbiter images. We focus on two aspects that appear to be new in photogrammetry applied to Mars: the correlation method and the relative orientation method based on the analysis of the coplaneity. We demonstrate that a DEM from two Viking Orbiter images selected according to the criteria of parallax height ratio, resolution, and illumination conditions can be obtained with a better spatial lateral resolution than the interpolated Mars Orbiter Laser Altimeter (MOLA) data. Two images on the Deuteronilus Mensae area satisfy these criteria. The algorithm of reconstruction of the topography is divided into two steps. First, a correlation method at subpixel level based on Fourier transform has been adapted to extract image location of homologous points. As initial camera positions and angles are known with poor accuracy, an adjustment is then carried out by the method of relative orientation which requires no absolute control points and takes advantage of the subpixel accuracy of our matching algorithm. Finally, the DEM is projected into a Martian reference frame and registered to MOLA data for validation. The spatial resolution depends on the roughness of the area and is locally close to the distance between individual MOLA footprints along one pro le (300 m), as demonstrated from the comparison with MOLA pro les. These rst results of high-resolution DEM from Viking Orbiter images should be used with MOLA data for many geological studies that require high-resolution topographic data.

1. Introduction

Mars topographic mapping is essential to understand the physical processes in geology. Long-wavelength topography re ects processes that have shaped Mars on a global scale. Topography and gravity measured by the Mars Global Surveyor have enabled the determination of the global crust and upper mantle structure of Mars and give insights in the thermal evolution of the planet [Zuber et al., 2000]. The northern hemisphere depression is also demonstrated to be the result of an internal process [Smith et al., 1999]. Topographic data on both poles of Mars have been used for volatile inventory [Schenk and Moore, 2000; Clifford et al., 2000; Zuber et al., 1998] and investigations concerning the past extent of the north polar cap [Fisbaugh and Head, 2000]. Short-wavelength topography is controlled by crustal, impact, and climatic processes. Applications of topographic data at those scales are numerous. Topography across extensional and compressional features can be used to estimate magnitudes and distributions of strain. The shape of volcanism constructs and lava ow constrains volumes, styles, and mechanics of Martian volcanism. The detailed shapes of simple and complex impact craters provide information on the mechanics of impact, nature of the impacted substrate, and subsequent modi cation processes. The dynamic of erosional processes on Mars involving wind, water, and ice requires kowledge of the topography and topographic gradients. Thus, there is great interest in realizing both local and global high-resolution digital elevation models (DEM). The objectives of this work are de ned with respect to previous published works concerning the Martian topography obtained by photogrammetry. At the global scale of the planet, a global topographic map has been derived from stereophotogrammetric analysis of Viking Orbiter images [Wu et al., 1986]. The resolution of this topographic map is 1 km at the equator, and the vertical precison ranges from 500 m at the equator to 2.5 km at the poles [Zuber et al., 1992]. The U.S. Geological Survey control points network has been recently updated from the Mars Path nder tracking data [Zeitler and Oberst, 1999a]. The accuracy of these control points is quoted as 800 m at 1. Frow these new data, they have reanalyzed the shape of Mars and have made a global DEM of the whole planet with 50-km grid spacing [Zeitler and Oberst, 1999b]. The main data set of Mars topography is provided by the Mars Orbiter Laser Altimeter (MOLA). In fall 1997, MOLA acquired the

2 rst 18 pro les across the northern hemisphere of the planet. The acquisitions are now more than 800,000 measurements a day with a vertical absolute accuracy of 30 m and a measurement spacing along track of 300 m [Zuber et al., 1998]. MOLA is still in the acquisition phase, and the data acquired to February 2001 have been released. While MOLA will provide extremely well controlled high-precision topographic pro les over the surface of Mars, spatial coverage is limited by the number of ground tracks. A 0:005o by 0:005o global topographic grid has been released. However, signi cant gaps of up to several kilometers are present at low latitude between pole-to-pole ground tracks. Few stereo-DEMs from Viking Orbiter images have been produced in the last decades. Two attemps of stereophotogrammetric mapping of the south pole of Mars have been found in the literature [Dzurisin and Blasius, 1975; Herkenho and Murray, 1990]. In these studies, sparse and low-accuracy data were obtained. Stereophotogrammetry of Viking Orbiter images was used to map the Ma'adim Vallis area [Thornhill et al., 1993]. This rst automated method requires an initial selection of tie points, and the vertical resolution of the DEM was about 700 m. A DEM of the Ius-Tithonium Chasma [Cook et al., 1999] has been produced by the Extraterrestrial Orbital DEMs for Understanding Surfaces (EXODUS) project at University College London. Owing to the low resolution of the images and the large patch size used for matching, the e ective resolution of Cook's model [Cook et al., 1999] is comparable to the interorbit spacing of MOLA tracks. In order to characterize both the general and detailed topography of the landing site zone for the Mars Polar Lander Mission, a stereo-DEM of the south pole of Mars was generated recently from Viking Orbiter images. However, a smoothing function is applied on these DEMs to remove noise of several hundred meters [Schenk and Moore, 2000]. The purpose of this paper is to present a method for computing a local high-resolution DEM derived from Viking stereoimages. The method is divided into two steps (Figure 1): a speci c correlation method adapted to Viking stereopairs and the restitution of the DEM. After choosing suitable stereopairs, homologous points in the pair of images have to be extracted by a phase correlation algorithm based on bidimensional Fourier transform. Then adjustment of the camera parameters and computation of the DEM have to be achieved from the pixel positions of homologous points on each image. Then we propose a

3 method to register the Viking DEM and MOLA data. As a rst application, a DEM is presented on the area of Deuteronilus Mensae. Mesas, scarps, and craters present on this area have strong topographic signatures. These results are registered to MOLA data for validation of the method and estimation of its accuracy.

2. Digital Elevation Model Processing Mathematical relations of photogrammetry [Kraus and Waldhausl, 1998] can be applied on Viking Or-

biter images, after correction of the distortion. DEMs generation algorithms from aerial or satellite images have been developed for terrestrial applications for decades. The correlation method and the adjustment method we present are not new in photogrammetry. However, these methods have never been applied to Viking Orbiter images.

2.1. Viking Orbiter Images Selection and Preprocessing

The quality of a DEM derived from a stereopair depends on six limiting factors: orbital data, pointing, geometric calibration, illumination conditions, parallax/height ratio, and resolution. In the case of Viking missions, pointing and orbital positions are poorly known and must be updated. The visual imaging subsystem (VIS) of Viking Orbiters consists of two vidicon framing televisions. Each telescope has a 475-cm focal length and a vidicon of 37 mm-diameter [Snyder, 1979]. The major e ect that a ects the geometric calibration for narrow angle optics and vidicon camera is beam bending or electronic distortion. In a vidicon camera, illumination ground point information is transferred into pixel value through a charge process. In the case of a bright area, the electron beam is de ected away from the bright pixel: it is the beam bending e ect. Errors in location measurements on the images can be up to 4 pixels (T. Duxbury, personal communication, 1999). A method to address this problem is presented at the step of restitution of the DEM. The illumination angles should be similar on each image to have the same shadow pattern and contrast. The resolution of both images has to be close enough to exhibit the same level of details. The parallaxheight ratio, which is de ned as the di erence in parallax by unit height between both images, has to range from 0.01 to 1 in order to have both stereoe ect and a good correlation of the ground features.

Viking Orbiters 1 and 2 collectively acquired more than 48,000 images of Mars' surface between 1976 and 1980 [Snyder, 1979]. The acquisition conditions of Viking Orbiter images cover a broad range of illumination conditions, incidence angles, and resolutions. The sizes of areas imaged range from several kilometers to several thousands of kilometers. This heterogeneous data set is explained by the numerous scienti c objectives, technical reasons, and consequently changes in acquisition trajectories [Snyder, 1977]. All these elements imply that Viking Orbiter images are not optimized for stereophotogrammetric mapping. Thus, from all the images with stereocoverage, in the case of perfect knowledge of pointing, position, and beam bending, the expected vertical precision (EP) ranges from less than 10 m to more than 1 km [Kirk, 1999]. About 10% of the surface of Mars is covered by stereopairs with EP of less than 100 m. This is also true at low latitude where the topographic grid derived from MOLA data has the lowest resolution. Two images have been selected (Table 1) on Deuteronilus Mensae where the knowledge of the topography can be used for the geological interpretation of landforms probably due to ice in the subsurface [Mangold and Allemand, 2001]. A rst processing, using the Integrated Software for Imagers and Spectrometers [USGS, 1999] is applied to the images in order to remove the reseaus, to calculate a model of the camera distortion for each image, and to remove most of the electronic noise and to enhance the contrast.

2.2. Initial Homologous Points Processing

In order to avoid problems with the curvature of the planet, the DEM is processed in a Cartesian planetocentric frame. Hereinafter, we call homologous points the corresponding pixels on the two images and absolute control points the points known with absolute coordinates. Assuming that each point is on the reference ellipsoid (a two-axial ellipsoid is used with an equatorial radius equal to 3396.0 km and a polar radius equal to 3378 km [Zuber et al., 1998]), the geometric problem is to nd homologous points on the two images. This rst set of pixel positions corresponds to approximate homologous points computed without elevation to be used as an initial guess to nd exact pixel positions of homologous points on images. A pixel position on the rst image corresponds to a line of sight and can be expressed by a unit vector in the camera frame. The camera frame is de ned

4 by the line of sight at the center of the image and the sides of the image. Camera angles and spacecraft positions are given in the equatorial reference frame (i.e., the J2000 frame), de ned by the projection on the celestial sphere of the equatorial plane and Earth to Sun vector on spring equinox 2000 [Davies et al., 1992]. These values are given in the Spacecraft Planet Instrument C-matrix Event (SPICE) les [USGS, 1999] processed by the Jet Propulsion Laboratory. A rotation matrix (the C-matrix, Appendix A), is computed from camera angles and allows the transformation of pixel position into vectors expressed in the J2000 frame. This line of sight is then projected into the planetocentric frame using the rotation matrix derived from planet angles de ned in the J2000 frame. The intersection between the line of sight and the reference ellipsoid is nally computed in the planetocentric frame. Thus, a coordinate on the ellipsoid is virtually attributed to each pixel of the rst image. The inverse procedure is applied to nd the position of virtual homologous points on the second image. At the end, a rst set of homologous points is obtained as an initial guess for the correlation algorithm.

2.3. Phase Correlation Algorithm

We present here an adaptation of the phase correlation method and discuss the speci c problems with the utilization of Viking Orbiter images. This method gives the best results and requires less computation time than most conventional correlation methods. The phase correlation method appeared in the literature in 1975 [Kuglin and Hines, 1975]. Its superiority to conventional correlation was noted and demonstrated [Schaum and McHugh, 1991]. Its superiority is double when it is compared with other methods: the location of the correlation peak is more precise and the correlation peak is narrower. A translation of the coordinate frame between the two images is re ected in the Fourier domain purely as a linear shift in the phase. If xn describes the intensity of the nth sample on the rst image and yn = xn;s represents the intensity of the same sample on the second image frame, which is shifted by s samples, then Yk = Xk exp( ;i2sk (1) N ); where Xk and Yk are the discrete Fourier transform of xn and yn , and Xk =

X

N ;1 n=0

xn exp( ;i2nk N );

(2)

where N is the number of samples in the image. Equation (1) holds generally as N ! 1 because the discrete Fourier transform becomes the continuous Fourier transform. For nite N, equations (1) and (2) are strictly true only if images are periodic. Although a real image seldom has this property, the phase correlation method is still an accurate method, even if ksk is a large fraction of N [Schaum and McHugh, 1991]. The shift is obtained by computing the phase correlation function pn: N X Xk Yk exp( i2nk ); pn = N1 kXk kkYk k N k=0

(3)

where X is the complex conjugate of X. The normalization of the conjugate product is the key element of this technique. It produces a narrow correlation peak and enhances the peak to noise ratio. With this operation, all frequencies in the Fourier domain have the same weight for the computation of the phase shift. Thus the e ect of the phase shift is not limited to high-amplitude frequencies. The method is easily extended to the two-dimensional case in order to measure shift in both directions. In the case of stereopairs with the same geometry and resolution, small windows are de ned around the rst approximation of homologous points on each image. The correlation method estimates the lateral shift between the windows due to mean topography of the area. The smaller the window is, the better will be the spatial resolution of the DEM. However, the limiting factors will be the signal to noise ratio which decreases inversely with the number of points in the windows.

2.4. Application to Viking Stereopairs

The images of the selected area have di erent geometry and resolution (Figure 2). In order to apply this particular matching algorithm, the two images must be locally projected in the same geometry (resolution and azimuth). Expected shifts, due to the topography, are about a few pixels: windows with width of 16 pixels are used to satisfy overlapping area condition for these two images. A rst shift between the two windows is estimated by the phase correlation method. From the shift, location values of the homologous point on the second image is then obtained. However, the signal to noise ratio (SNR) of the correlation function is low for many windows,

5 but only a part of the Fourier spectrum is assumed to be a ected by some alterations and thus can be ltered. The pixel noise is expressed in the highest frequencies, and so high frequencies are removed. Even if radiometric correction have been processed [USGS, 1999], ground seasonal variations and varying response across the eld of view of the camera a ect the image at the scale of the window. Thus, largest wavelengths at the mean size of the window must also be removed. Finally, topographic signal is expressed at medium wavelengths which correspond to features that can be correlated on the window and producing the highest signal to noise ratio (Figure 3). An adaptative band-pass lter is applied to extract only signi cant wavelengths between 3 and 7 pixels, about one third of the spectrum. Correlation peak is clearly enhanced by this process (Figure 4). In order to attain a subpixel accuracy, the shift value is estimated by a weighted average of the 3x3 neighbors centered on the correlation peak [Michel, 1997]. If the shift of the two windows is larger than 1 pixel, nonoverlapping areas produce decorrelation of the signal. A second window is selected on the second image, shifted by the integer part of the rst estimated shift. The resulting overlapping area of the two windows becomes as large as possible, and the signal to noise ratio increases. The patch size at this nal step is 6-pixel radius on the st image, corresponding to 960 m. This operation is also a test of the correlation. The second shift has to be less than 1 pixel, and the signal to noise ratio has to increase. If not, the point is removed from the correlation grid.

3. Digital Elevation Model Restitution At the end of the correlation process, a data set of pairs of homologous points is obtained. The restitution of the DEM consists in the derivation of the three-dimensional coordinates of ground points from these two-dimensional locations on the images. A given pair of homologous points de nes two homologous rays. A ray is the direction de ned by the focal point of the camera and the considered ground point. As two homologous points represent the same ground point, two homologous rays have to intersect on the planetary surface. Both camera angles, positions and ground points are unknown. Three parameters for spatial location and three parameters for the angles have to be estimated for each camera. These parameters can be approximated from SPICE [USGS, 1999] les, but a DEM computed from the data set

of homologous points and approximated navigation data contains large errors. The accuracy of camera pointing and position depends on the accuracy of the gravity eld of Mars and the accuracy of the rotational state of the planet. These parameters de ne the transformation between the inertial frame and the cartographic coordinates. In worst cases, camera angles are known with a standard deviation of 0:25o. In order to increase the accuracy of the DEM, a global adjustment of camera angles and positions must be computed. An absolute orientation method would require absolute control points and homologous points [Kraus and Waldhausl, 1998]. The accuracy of control points in the case of Viking Orbiter images is about several hundred meters [Zeitler and Oberst, 1999b]. Absolute control points are generally topographic features which have high-frequency brightness variations. The images of these points on the vidicon are then a ected by a strong beam bending, and their positions are then not precise. There is no sense to compute intersection of homologous rays from subpixel accuracy measurements if the orientation has been achieved from pixel locations with lowest precision. From these statements, a relative orientation method is developed for the Viking Orbiter images. The method takes advantage of the subpixel correlation results and is able, as we demonstrate below, to select points that are devoid of beam-bending residuals (this is theoretically possible close to the reseaus and when brightness variation is low enough). It is emphasized here that this orientation process does not allow us to obtain ground points in absolute coordinate, longitude, latitude and altitude. An absolute orientation must be processed later, anyway, to evaluate the ground points positions (digital terrain model (DTM)) on the planetary reference surface.

3.1. Relative Orientation of the Camera

In the relative orientation process, ground points and navigation data are the unknowns of the problem. The equations of intersection between homologous rays are solved by the least squares minimization method [Kraus and Waldhausl, 1998]. The minimization method requires linearization of equations. The derivatives of pixel position have to be estimated. The pixel positions on each image are expressed in the planetocentric frame as functions of all the unknowns: camera angles, spacecraft positions, and coordinates of ground points:

6 1 = f (X1 ; Y1; Z1; 1; 1 ; !1; X; Y; Z) (4) 1 = f (X1 ; Y1; Z1; 1; 1; !1; X; Y; Z) (5) 2 = f (X2 ; Y2; Z2; 2; 2 ; !2; X; Y; Z) (6) 2 = f (X2 ; Y2; Z2; 2; 2; !2; X; Y; Z); (7) where j , j are pixel positions on images taken by camera j of ground point (X; Y; Z) and Xj , Yj , Zj , j , j , !j are camera positions and angles. The general linearized function is expressed: j ; j0 =f (Xj0 ; Yj0; Zj0 ; 0j ; j0 ; !j; X 0 ; Y 0 ; Z 0) + @f + @f + @f + @f + @f + @f + @Xj @Yj @Zj @ j @j @!j @f + @f + @f @X @Y @Z j ; j0 =f (Xj0 ; Yj0; Zj0 ; 0j ; j0 ; !j; X 0 ; Y 0 ; Z 0) + @f @f @f @f @f @f @Xj + @Yj + @Zj + @ j + @j + @!j + @f + @f + @f ; @X @Y @Z

where exponent 0 holds for initial approximated values. With no control points, only ve parameters can be found by the relative orientation [Kraus and Waldhausl, 1998]. The six parameters of the rst camera and one parameter of the second camera (one coordinate of location) are held xed. The ve remaining parameters of the second camera (two coordinates of location and three angles) and ground points are estimated by the minimization method. The inverse method minimizes the deviation between pixel positions estimated from the correlation process and computed from newly estimated ground points at xed orientation parameters. Mathematically, the square matrix, in which the number of elements is proportional to the square of the number of ground points, must be inverted. Owing to computer limitations, the inversion can only be computed for about a hundred points. A criterion must be applied for the selection of the most accurate homologous points. This is the key point for a successful relative orientation. In the case of Viking Orbiter images, acquired with vidicon camera, this mainly implies nding homologous points correlated with subpixel accuracy while beam bending is low enough to not a ect the measurements. The reseaus are used to correct the beam-bending e ect, but if one has highfrequency brightness variations between reseaus, the

errors can easily vary from 0 to 4 pixels (T. Duxbury, personal communication, 2000). The method implemented for point selection follows. Homologous rays intersect on the surface model. The unit vectors of the two homologous rays and the vector de ned by the camera positions are in the same plane and follows

; ! ;;;! ;;;! O;; 1O2:(O1P1 ^ O2P2) = 0;

(8)

! where;;; O!1 and O2 are the camera positions and ; O;; 1 P1 and O2P2 are the unit vectors corresponding to the two homologous rays. When camera parameters contain small errors, equation (8) does not apply. This relation in this case can be expressed as

! ;! ;;;! ;;;! coplaneity = (; O;; 1 O2 + 12)(R1 O1P1 ^ R2 O2P2); (9)

where ; ! 12 is the translation vector due to errors in camera positions and R1 and R2 are rotations matrix due to small errors in positions and angles. This relation is a function of the intersection angle of homologous rays which depends on pixel position. The general trend of the coplaneity along the line of the image can be observed in Figure 5. Distribution of these points is explained both by camera pointing and errors in pixel locations. From (8), the general trend along the line de ned at the rst order by a straight line is produced by the errors of orientation and pointing. Deviations around this straight line are interpreted as errors in pixel locations. The origins of these last errors are explained by the failure of the correlation algorithm and badly detected reseau and are in most cases due to beam bending. For those points the distribution above and below the straight line is the result of brightness variations: a bright spot will deviate the beam in opposite directions for two points symmetrically located from its position and will result in errors of opposite sign. Thus, the electronic distortion produces errors randomly distributed around values according to equation (9). From this statement, a least squares t with a second-degree polynomial model of the coplaneity is processed along each line of the image. Points which are best tted by the model are expected to be the most accurate ones. This statement has been checked by moving pixel locations: when the values of coplaneity are very low, any small modi cations result in a fast increase of many orders of the coplaneity. Failure

7 of the coplaneity algorithm produces points with the largest errors. Thus, the signal to noise ratio obtained from the correlation of each point is used as weight in the tting process. The relative orientation is nally accomplished with about a hundred points uniformly distributed on the image. We can observe residuals on the camera (j ; j0 ; j ; xi0j ) in pixel units in Figure 6, before and after the orientation process. It has been noted that with only two images, errors in pixel locations due to beam bending can produce errors in ground point positions without an increase of the error computation (T. Duxbury, personal communication, 2000). Simulated beam bending of one pixel has been introduced randomly in the whole data set, and it is checked for how residuals increase for points initially at the subpixel accuracy level. We can say nothing for the points with initial high residuals (more than 1 pixel). This means that residuals increase with beam bending, only if the orientation has been achieved from points that only contain subpixel (0-0.3) beam bending. Intersection distances between lines of sight are computed for the whole data set from initial orientation parameters and adjusted ones (Figure 7). This demonstrates the e ect of the low accuracy on the pointing and positions of cameras. This demonstrates also the consistency between the data set of homologous points and the adjusted orientation parameters. The location where distance between lines of sight is minimum is an approximation for a ground point position. These values are inputs for the inverse method used to compute the nal ground point positions.

3.2. Absolute Orientation

DEM obtained by this process is not referenced and has to be projected into the MOLA reference frame (aerocentric, longitude positive west (International Astronomic Union, 1991)). Mathematically, 7 parameters, held xed during the relative orientation, from the 12 have to be determined. Again, a minimization method is applied and the geometric transformation is linearized:

;;! ;! ;! OM = T + m  R(; OM); (10) where M 0 is the transform of M, ;! T is the translation 0

vector, R is the rotation matrix, and m is the scale factor. This transformation can be expressed as:

;;! ;!; Xu; Yu; Zu ; !; ; ; m); OM 0 = f(; OM

(11)

where (Xu ; Yu ; Zu) are the coordinates of the translation vector; (!; ; ) are the angles of rotation relative to each axis; and m is the scale factor. Absolute control points are required to estimate the transformation. The orientation is achieved in two steps. With all the MOLA points available in the area, a DEM is computed. A rst set of absolute control points is derived from the comparison of topographic features on the Viking Orbiter image and on the MOLA DEM. A rst transformation is computed that minimizes deviation between MOLA points and transformed points on the Viking DEM. This rst adjustment gives good results concerning the translation and the scale factor and allows a rough comparison between the two DEMs. The residuals in elevation can be up to 400 m because of the interpolation of MOLA data. To improve the absolute orientation, individual MOLA points must be used. Thus, in a second step, a new set of absolute control points is extracted from the comparison of MOLA topographic pro les and the projected Viking DEM. To do this, topographic pro les are computed from the Viking DEM at the exact location of MOLA ground tracks and MOLA individual points. Absolute control points are then selected directly on each pair of MOLA and Viking topographic pro les, and a new transformation is processed. In the second case, the residuals at the end of the computation of this second, more accurate transformation are all below 200 m. The nal residuals are explained by the errors of the Viking DEM and the errors coming from the rst transformation. However, the nal errors in positioning are estimated to be less than 1 km. The DEM is nally processed by interpolation and gridding with Delaunay triangulation. High-frequency noise is removed by smoothing the DEM slightly.

4. Discussion The DEM derived from images 529a05 and 675b62 is presented in this section. Seven topographic pro les are extracted for comparison with MOLA data. Other topographic pro les are extracted in roughly the E-W direction, where MOLA data have to be interpolated to produce a DEM.

4.1. Visual Inspection of the DEM and Estimation of Errors

The correlation algorithm has been applied to the two selected images. The process of correlation is achieved for each pixel of the rst image (529a05).

8 The SNR of the correlation is signi cant for 96% of the pixels of the overlapping area. Around 45% of this set of points are used to process the DEM. However, the density of points is close to 100% on the rough mensae, while many points are missing in the

at area, because of the lack of texture for the stereomatcher eciency. Relative orientation is processed according to the method described above. 120 control points selected from the coplaneity have been used in this process. The DEM is then projected on the MOLA reference frame, and the Viking image is warped on it (Figure 8). An estimation of the accuracy is obtained from the minimization step (see Appendix B). Furthermore, maps of errors and resolution can be processed. The resolution map is de ned by the mean distance between points used for interpolation (Figure 9). We can observe a high concentration of points on the mensae which are the roughest areas. Two error maps are then computed: altitude errors and planimetric errors. The errors displayed are estimated by averaging errors for each point inside sliding windows. Estimation of the standard deviation of errors (generally about 2%) in the sliding windows con rms that variation of the errors in location measurements over the image are relatively smooth. These maps show mainly that the DEM's accuracy ranges from 10 m to 200 m and mean accuracy is 60 m, which is slightly greater than the apriori value estimated from the parallax-height ratio (58 m).

4.2. Comparison With MOLA Data

In order to validate our method, topographic pro les are extracted from the DEM (Figure 10) along ground tracks of MOLA pro les. Deviations are consistent with the error maps. In order to quantify the agreement between the Viking data and the MOLA data, we compute the correlation coecient for each pro le presented. The average correlation coecient is 0.96. The better agreement is obtained along the scarps and generally along slopes, while at areas appear more noisy. These errors are seen as random oscillations around the MOLA pro le. The beam bending and other sources of errors in image location measurements explain errors both in altimetry and planimetry. However, actually, it is not possible to identify the part of the errors due to the uncertainty in the absolute orientation. In conclusion, good agreement is obtained between the two data sets, considering the expected vertical precision.

4.3. Examples of High-Resolution Features of Our DEM

Four pro les are extracted along small features in the W-E direction (Figure 11), perpendicular to MOLA pro les. The reference (P1, P2, P3, P4) corresponds to the starting point of the cross section. The arrows underline the structures of a few hundred meters scale seen on the image and on the topographic pro les. On P2, a small crater of 200-m depth on the mesa is detected. A knob, which may have been carried with the deposits, appears to have the same elevation as the mesa (P3). At the limit of accuracy of the DEM, a little partially lled trough(P4) is also identi ed in our DEM.

5. Conclusion Both correlation process and adjustment of orientation parameters have been adapted to the speci cs of Viking Orbiter images, in order to derive a DEM of the Martian surface. The large proportion of points correlated allows interpolation of topographic grids with a spatial resolution close to the image resolution. However, topographic details smaller than the patch size of correlation should be theoretically blurred. The relative orientation method requires no absolute control points and takes advantage of the subpixel accuracy of the correlation algorithm where beam bending is low. In other methods, absolute ground coordinates were used in the orientation process. Uncertainties on these ground control points of several hundred meters and possible uncertainties on locations of these control points on images produce less reliable orientation than relative orientation from homologous points known at the subpixel level. Thus, the accuracy obtained by the method de ned in this paper is impossible to obtain from a one-step absolute orientation method using ground control points, even if correlation were possible at the subpixel level. Consequently, the accuracy in relative elevations has never been attained before from Viking images. The registration between MOLA data and Viking images is now possible as topographic data derived from Viking images are used for this. From these promising rst results, we plan to compute new DEMs according to the possibilities given by the Viking stereocoverage and applications to geologically relevant questions. This technique applied in the special case of Viking imagery can be relevant to many other planetary missions and Earth topics.

Appendix A: Mathematical Expression of Camera Frame Transformation Matrix (C-Matrix)

The camera angles de ned orientation of the line of sight in the J2000 reference frame. The C-matrix allows the transformation from the camera frame (pixel positions) and the J2000 frame:

0r r r 1 @r1121 r1222 r1323A ;

where L is the deviation between the initial model and observations, xb = (dx1; dx2; dx3, ... ,dxn) is the correction to apply to the initial model and A is the matrix of error equations. If v is the system of error equations, v = Axb;

(B3)

the condition for minimization is as follows min(vT v) = (Axb ; L)T (Axb ; L) = xbT AT Axb ; 2LT Ax + LT L:

r31 r32 r33

The solution is computed by adding a small correction to the initial point:

where r11 = ; sin  cos ! ; cos  sin sin!; r12 = cos ! ; sin  sin sin !; r13 = cos sin!; r21 = sin  sin ! ; cos  sin cos !; r22 = ; cos  sin ! ; sin  sin cos !; r23 = cos cos !; r31 = cos  cos ; r32 = sin  cos ; r33 = sin ; is instrument right ascension;  is instrument declination and, ! is instrument twist:

xb = (AT A);1 AT L:

(B4)

Taking into account the precision of measure is possible by computing the following weight matrix: 0 12 1 ::: 1 BB 12 CC ::: 2 B C 1 P =B B@ ::: ::: :::32 :::::: ::: CCA ; ::: 1n2 where i is the standard deviation. The weighted solution is expressed by: xb = (AT PA);1 AT PL:

Appendix B: Minimization Method

Given the following set of nonlinear equations

li = f(x1 ; x2; x3; :::; xn); (B1) this system is rst linearized by a rst-order Taylor development of an approximate solution: li =f(x01 ; x02; x03; :::; x0n) + @f )0dx + ( @f )0 dx + ( @x 1 @x2 2 1 @f )0dx + ::: + ( @f )0 dx : ( @x 3 @xn n 3 Then these linear equations can be written in a matrix form: L = Axb;

9

(B2)

With n equations and u unknowns the standard deviation for all the unknowns can be computed from the remainders:

p

xk = 0  Qkk; where

s

T 0 = Vn ;PVu

Q = (AT PA);1:

Acknowledgments. We are grateful to T. Duxbury for the review of an earlier manuscript about this method. We acknowledge three anonymous referees for their thoughtful remarks. This work was supported by the Programme National de Planetologie of the Institut National des Sciences de l'Univers from France.

References

Cli ord, S.M., et al., The state and future of Mars polar science and exploration, Icarus, 144, 210-242, 2000. Cook, A. C., T. Day, and J.-P. Muller, The UCL IusTithonium Chasma Digital Elevation Model in Proceedings of the Fifth International Conference on Mars, Lunar and Planetary Institute, Houston, Pasadena, Calif., 1999 Davies, M. E., R. M. Batson and S. S. C Wu, Geodesy and cartography in Mars, in Mars, pp. 321-342, Univ. of Ariz. Press, Tucson, 1992. Dzurisin, D., and K. Blasius, Topography of the polar layered deposits of Mars, J. Geophys. Res., 80, 32863306, 1975. Fisbaugh, K. E. and J. W. Head, North polar region of Mars: Topography of circular deposits from Mars Orbiter Laser Altimeter (MOLA) data and evidence for asymmetric retreat of the polar cap, J. Geophys. Res., 105, 22,455-22,486, 2000. Herkenho , K., and B. Murray, High-resolution topography and albedo of the south polar layered deposits, J. Geophys. Res., 95, 14,511-14,529, 1990. Kirk, R., A database of Viking orbiter image coverage of Mars for cartographic and scienti c use, in Proceedings of the Fifth International Conference on Mars, Lunar and Planetary Institute, Houston, Pasadena, Calif., 6217, 1999. Kraus, K., and P. Waldhausl, Manuel de Photogrammetrie: Principes et Procedes Fondamentaux, 407 pp., Hermes, Paris 1998. Kuglin, C., and D. Hines, The phase correlation image alignment method, in Proceedings of the Conference on Cybernetics Society, IEEE Press, Piscataway, N.J., 163-165, 1975. Mangold, N., and P. Allemand, Topographic analysis of features related to ice on Mars, Geophys. Res. Lett., 28, 407-410, 2001. Michel, R., Les mesures de mouvements par imagerie SAR et leur application en glaciologie et en sismotectonique, Ph.D. thesis, Univ. de Paris-Sud, Paris, 1997 Schaum, A., and M. McHugh, Analytic method of image registration: Displacement estimation and resampling, NRL Tech. Rep. 9298, Nav. Res. Lab., Washington, D. C., 1991. Schenk, P.M., and J. M. Moore, Stereo topography of the south polar region of Mars: Volatile inventory and Mars Polar Lander landing site, J. Geophys. Res., 105, 24,529-24,546, 2000. Smith, D. E., et al., The global topography of Mars and implications for surface evolution, Science, 284, 14951503, 1999. Snyder, C. W., The missions of the Viking Orbiters, J. Geophys. Res., 82, 3971{3983, 1977.

10 Snyder, C. W., The extended mission of Viking, J. Geophys. Res., 84, 7917{7933, 1979. Thornhill, G. D., D. A. Rothery, J. B. Murray, A. C. Cook, T. Day, J.-P. Mullerand J. C. Ili e, Topography of Apollinaris Patera and Ma'adim Vallis: Automated extraction of digital elevation models, J. Geophys. Res., 98, 23,581-23,588, 1993. U.S. Geological Survey, Isis documentation, Washington D.C., 1999. Wu, S. C. C., R. Jordan and F. J. Schafer, Mars global topographic map, scale 1:15,000,000, NASA Tech. Memo., TM-88383, 614-617, 1986. Zeitler, W., and Jurgen Oberst, The Mars Path nder landing site and the Viking control points network, J. Geophys. Res., 104, 8935-8941, 1999a. Zeitler, W., and J. Oberst, The shape of Mars before global surveyor: Results from reanalysis of the Viking control points network, J. Geophys. Res., 104, 14,05114,063, 1999b. Zuber, M. T., D. E. Smith, S. C. Solomon, D. O. Mulheman, J. W. Head, J. B. Garvin, J. B. Abshire and J. L. Bufton, The Mars Observer Laser Altimeter investigation, J. Geophys. Res., 97, 7781-7797, 1992. Zuber, M. T., D. E. Smith, R. J. Philipps, S. C. Solomon, W. B. Banerdt, G. A. Neumann, and O. Aharonson, Shape of the northern hemisphere of Mars from the Mars Orbiter Laser Altimeter (MOLA), Geophys. Res. Lett., 25, 4393-4396, 1998. Zuber, M. T., et al., Internal structure and early thermal evolution of Mars from Mars Global Surveyor topography and gravity, Science, 287, 1788-1793, 2000.

D. Baratoux, C. Delacourt, and P. Allemand, Laboratoire des Sciences de la Terre, UMR 5570, UCBL Lyon I, CNRS et ENS-Lyon, 69 622 Villleurbanne, Cedex, France. ([email protected];[email protected],[email protected])

Received January 15, 2001; revised July 6, 2001; accepted August 29, 2001.

This preprint was prepared with AGU's LATEX macros v5.01. File jgrrep formatted February 27, 2002.

11

Table 1. Image Parameters with Expected Vertical Precision Estimated (With the Assumption of 0.2-Pixel Accuracy Attained From the Correlation) of 56 m Image

Resolutiona,

Latitude Rangeb, m

Longitude Range, deg

Incidence, Anglea, deg

Emission, anglea, deg

North azimuthc, deg

529a05 675b62

162 200

38:78 ; 42:76 38:55 ; 43:75

339:85 ; 344:37 337:91 ; 343:90

67:74 72:77

18:31 20:99

186:63 168:28

a Values computed at the center of the image. b Latitude is aerographic. c North azimuth is the angle between the north

and the line de ned by the horizontal line of the image that joins the center of the image and the middle of the right side.

12

Global adjutsment process

Correlation process

Pair selection of Viking images

Initial set of approximated homologous point Correlation algorithm based on Fourier transform Data set of homologous point sub-pixel accuracy

Coplaneity fitting

Selection of the most accurate homologous points Adjustment by relative orientation

New adjusted camera positions and angles Ground points computation

Absolute orientation with MOLA Interpolation and gridding

Figure 1. General overview of DEM processing method applied to Viking images .

13 N

50km

50km

N

A.

B.

Figure 2. (a) Viking image 529a05; resolution of 160 m per pixel (b) Viking image 675b62; resolution of 200 m

Signal To Noise Ratio

per pixel. The overlapping area is depicted by a bold line.

Wave length (pixel)

Figure 3. General trend of signal to noise ratio as a function of the center of the band pass- lter. It is generally

observed that a maximum signal to noise ratio of around 3 to 7 pixels corresponds to a meaningful correlation; inverse Fourier transform of these frequencies shows features on images that are expected to be correlated.

14 Phase correlation function with band pass filter

Filtered second image

Filtered first image

Phase correlation function

Filtered second image

Filtered first image

A.

B.

Figure 4. Enhancement of the correlation peak by band-pass ltering. (a) Correlation peak over rugged area. Image 2 f675b62

Image 1 f529a05

(b) Correlation peak over featureless area. In this case, the ltering has no e ect, and the correlation failure will produce a gap in the DEM or a lower spatial resolution. 4.0 2.0

0.0

0.0

Residuals (pixel unit)

Coplaneity

2.0

2.0 2.0 4.0 Pixel along the line of the image

Figure 5. General trend of the coplaneity along one line of the image.

15 Residuals before adjustement

Residuals in pixel units on image 2 (f675b62)

1.0 0.8 0.6 0.4 0.2 0.0 0.0

0.6 0.4 0.8 0.2 Residuals in pixel units on image 1 (f529a05)

Residuals after adjustement

Residuals in pixel units on image 2 (f675b62)

0.5 0.4 0.3 0.2 0.1

0.0

0.1 0.2 0.3 0.4 Residuals in pixel units on image 1 (f529a05)

Figure 6. Residuals before and after the relative adjustment in pixel units, for both images. 35000

Before adjustment

30000

Number of points

After relative orientation 25000

20000

15000

10000

5000

0

0

1

2

3

Standard deviation of the positions of ground points (km)

Figure 7. Intersection distance between lines of sight before and after the minimization is applied.

16

50 km Figure 8. Overlay of the image on the digital elevation model.

17 < 250m 250m-300m 300m-350m 350m-400m 400m-450m 450m-500m 500m-750m 750m-1m > 1km

Spatial resolution of the DEM 0 km

50 km

Elevations above the reference areoid

A. -3900 m

B. 0m

< 50m 50m-60m 60m-70m 70m-80m 80m-90m 90m-100m 100m-125m 125m-150m > 150m

C.

D.

Figure 9. (a) DEM in gray levels obtained from Viking images 529a05 and 675b62. (b) Map of the resolution of the digital elevation model. (c) Map of the errors in altimetry. (d) Map of the errors in planimetry.

18 -1.5

13073

10740

11262

12639

Elevations (km)

14777

42.5

Orbit n˚ : 13903

Orbit n˚ : 11262 -1.0

C.C : 0.948

C.C : 0.964

-1.5

-2.0

-2.0

-2.5

-2.5 -3.0

-3.0

-3.5

Mola profile Viking DEM

-4.0 39

-3.5

Mola profile Viking DEM

-4.0 41 39

40

40

41

13384

Aerocentric latitude Aerocentric latitude -1.0 Elevations (km)

13903

Orbit n˚ : 14777

-1.5

C.C : 0.975

-2.0 -2.5 -3.0 -3.5

38.5 344.0

A.

340.0

Elevations (km)

Elevations (km)

Orbit n˚ : 12639 C.C : 0.956

-2.0 -2.5 -3.0 Mola profile Viking DEM

39

40

41

42

43

C.C : 0.976

-2.5 -3.0

39

40

41

42

43

-1.0 Orbit n˚ : 13384

Orbit n˚ : 13073

C.C : 0.968

Elevations (km)

Elevations (km)

Mola profile Viking DEM

Aerocentric latitude

-2.0 -2.5 -3.0

-4.0 38

42

-2.0

-4.0 38

-1.0

-3.5

41

Orbit n˚ : 10740

-1.5

Aerocentric latitude -1.5

40

Aerocentric latitude

-3.5

-3.5 -4.0 38

Mola profile Viking DEM

-1.0

-1.0 -1.5

-4.0 39

40

C.C : 0.943

-2.0 -2.5 -3.0 -3.5

Mola profile Viking DEM

39

-1.5

41

42

Aerocentric latitude

43

-4.0

Mola profile Viking DEM

38

39

40

41

42

43

Aerocentric latitude

Figure 10. Comparison between MOLA pro le and pro le extracted from the Viking DEM. (top left) MOLA

ground tracks on the MOLA-derived DEM. The cylindrical projection is used, latitude is aerocentric and longitude is west positive. C.C., correlation coecient computed between each MOLA and Viking topographic pro le.

19 2

Elevations (km)

P1

P1

P3

P2 P4

-4

-6 0

20

2

P4

P3

-2

-4

20

40

Distance (km)

Elevations (km)

Elevations (km)

0

0

60

40

2 P2

Elevations (km)

-2

Distance (km)

2

-6

0

0

-2

-4

-6 0

20

40

0

-2

-4

-6 0

Distance (km)

Figure 11. High-resolution topographic pro les extracted from the Viking DEM.

20

40

Distance (km)