IMPROVEMENT OF DEM GENERATION FROM ASTER IMAGES USING SATELLITE JITTER ESTIMATION AND OPEN SOURCE IMPLEMENTATION

The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) system embarked on the Terra (EOS AM-1) satellite has been a source of stereoscopic images covering the whole globe at a 15m resolution at a consistent quality for over 15 years. The potential of this data in terms of geomorphological analysis and change detection in three dimensions is unrivaled and needs to be exploited. However, the quality of the DEMs and ortho-images currently delivered by NASA (ASTER DMO products) is often of insufficient quality for a number of applications such as mountain glacier mass balance. For this study, the use of Ground Control Points (GCPs) or of other ground truth was rejected due to the global “big data” type of processing that we hope to perform on the ASTER archive. We have therefore developed a tool to compute Rational Polynomial Coefficient (RPC) models from the ASTER metadata and a method improving the quality of the matching by identifying and correcting jitter induced cross-track parallax errors. Our method outputs more accurate DEMs with less unmatched areas and reduced overall noise. The algorithms were implemented in the open source photogrammetric library and software suite MicMac.


INTRODUCTION 1.1 Motivation
The Terra (EOS AM-1) satellite was launched in December 1999 on a Sun-synchronous orbit.Aboard this satellite is the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) system.For more than 15 years, pairs of stereo images were collected by ASTER globally at a 15m resolution in the near infrared band, making its data the largest consistent multi-temporal dataset of stereo images available worldwide.
The ICEMASS project aims at analyzing glacier thickness changes using satellite data at a global level.Our aim is to use the incredible goldmine of data produced by ASTER for 2D/3D movement estimation and change assessment in glaciological research, providing several consistent snapshots through a time series.

Challenges
Even if the existence of such raw data is a great starting point, DEMs generated by NASA with SilcAst [SILC, 2015] (ASTER DMO products) do not provide a sufficient geometric quality for glacier volume change estimation over short periods, the expected change being significantly smaller than the accuracy of the product (a few meters against ±30m, see Figure 2 of the DMO over a sea ice (flat) scene).This is due to high frequency satellite jitter inducing attitude perturbation.Proposed causes for this disturbance include the mechanical cooling system, the rotation of mirrors and movements of the high gain antenna, but a source is yet to be formally identified.
The jitter yield at least three superimposed sinusoidal signals in both cross-track and along track directions in DEM comparison between the ASTER DMO and ground truth [Nuth and Kääb, 2011] or between ASTER and QuickBird images [Ayoub et al., 2008].The low acquisition frequency of the GNSS/IMU platform and star camera prevents the on-board estimation of jitter using the meta-data alone.The jitter induced waves were measured by [Kääb et al., 2013] as: • For the cross track : -Frequency around 4.6 km (average amplitude 1 m).
Figure 2: ASTER DMO product over the "Sea Ice" scene, which should be uniformly flat in the sea area, presents missing areas and high noise, the data peaks at several hundred meters.
A solution was proposed by [Iwasaki, 2011] to correct the jitter in the SWIR system (pseudo-nadir, short wave infra-red band) using the very short time delay between the bands S4 to S9.However, it cannot be applied to correct the VNIR system (pseudonadir and back-looking, green, red and near infra-red band) since "The VNIR subsystem is free from [band to band] parallax error because a dichroic filter divides incident light into each VNIR band."[Iwasaki, 2011].Also, only one band (near infra-red) is captured by the back-looking telescope.
Due to our global objectives and the focus on changing glaciers, it is not possible to use ground control points (GCPs) or ground truth obtained by other methods on stable terrain because such ground truth not being available (the SRTM mission only covering 60 • N to 56 • S for instance).

General description
We aim to identify and correct the impact of the composition of the jitters of the two bands on one of the bands while considering the other unafected.In our study, we chose to apply corrections to the band 3B since conserving the geometric integrity of nadir bands allows for color images to be generated.Our method evaluates the inaccuracies in co-registration due to imprecise orbital and orientation data of the Level 1A ASTER products without the use of ground control points or ground truth.After having extracted the images and meta-data, we start by applying radiometric corrections to the images.We then estimate RPC models (Rational Polynomial Coefficient) for both images of the stereo pair (bands 3N and 3B) using seed points taken at different altitudes from the lines of sight between the satellite and the lattice points, a method similar to the one described in [Lutes and Grodecki, 2008].From these RPCs and the images, we can estimate jitter corrections to apply to the images to correct for the jitters.We then apply these corrections to the backlooking image which leads to improvement in the correlation when finally computing an improved DEM.
2.2 From the raw data to images and RPCs 2.2.1 Destriping the images Each ASTER image is delivered with a calibration of the striping.It consists of three coefficients per column of the image, describing a linear function to be applied (see equation 1).The parameters (α, β , γ) of the equations are calibrated regularly in flight using a halogen lamp which light is uniformly radiated to the sensors [Abrams et al., 2002].This process is however not perfect since the sensors' response is not perfectly linear and their aging of the sensor is degrading the response even further.Some residual striping is indeed still present after correction.A more potent problem caused by such corrections is the further reduction of the image's dynamic range that was already limited because of the 8 bits bit depth, resulting in important posterization in the highlights and shadows.

Metadata to RPC
The metadata attached to the images contains for each band a number of records for the satellite position (from 12 to 16 depending on the scene and the band considered, usually more for the backlooking image than for the nadir image), to each of which is attached a line of 11 lattice points in image coordinates (column, row) and geographical coordinates (Longitude, Latitude) (see Figure 4).
For each line of sight defined by the association of a satellite position and one of its lattice point, we can define a collection of points regularly spaced on the line.This allows for the creation of layers of points forming a 3D grid of points in geographical coordinates (Longitude, Latitude, Ellipsoid Height) that are linked to points in the image (see Figure 4).
Figure 4: Satellite positions and associated lattice points in image and geographical coordinates and 3D GRID created from the lattice points.
These points are then used to estimate direct (from image to geographical coordinates, see Equations 2 and 3) and inverse RPC (from geographical to image coordinates, see Equations 4 and 5) models for the image through mean square matching ; they are rational function polynomial equations of the normalized image and geographical coordinates (scaled to a unit cube) defined as : With : The number of original lattice points (minimum 146) insures that the system is solvable.The 3D density of the generated grid makes the solution robust, with residuals in the order of 10 −6 in the unitary cube unit (10 −6 * (Coord max −Coord min ), equivalent to 10 −6 degrees or 10 −3 pixels) for both the direct and inverse RPC.

RPC to DEM
The RPC models are then used to determine the images coordinates for different candidates altitude for the points of a grid over the area of interest (or in other words describing the epipolar lines).For each candidate altitude (see Figure 5), the normalized cross-correlation score with a 5x5 window size is computed and the best candidate is kept.It yields DEMs of relatively good quality, but presenting significant noise (see the example in Figure 8 (top)).
Figure 5: Search for the altitude yielding the best correlation for a given position.

Error detection and quantification
The jitter of the satellite can be divided into two components : the cross-track and the along-track.Since the two images of each set are taken in the same track, the epipolar lines are almost parallel to the along track direction (y axis of the images).This results in the separation of both components of jitter in image space.
The effect of the cross-track component of the jitter is directly observable by performing a bi-dimensional correlation around the epipolar lines (see Figure 6), the potency of the vibration in pixels being equal to the distance between the point of maximum correlation (the homologous point) and the epipolar line.It is therefore possible to correct this element of the jitter.In Figure 7, we can see that the perspective ray from the nadir image can only cross the one from the backlooking image if the attitude of the satellites is known (the ray does not intersect with the rest of the red triangle).
The effect of the along-track jitter component however does not create any degradation in the correlation and is therefor not measurable in the raw line-of-site data.In Figure 7, we can see that the perspective rays from both images are crossing for all the possible attitudes represented by the blue triangle.Figure 7: Effect of the jitters on the direction of the perspective ray of the backlooking image (Nadir image considered stable).

Image correction
The cross-track parallax errors are exported as a grid (see Figure 6) with a correction factor for each pixel of the backlooking image, it can therefore be directly applied to the image using bilinear interpolation : Image Corr (x, y) = Image init (x − Grid ParallaxErr (x, y), y) (7) Once the image is corrected, the DEM generation can be performed using the method presented in section 2.3.The process returns better correlation scores over the entire area, and the DEM produced is significantly less noisy compared to the DEM obtained from the uncorrected images (see Figures 8 (bottom) and 9).
Figure 8: DEM product over the "Sea Ice" scene.Top : using only the metadata.Bottom : after correcting the cross-track parallax.

Implementation
Our implementation uses the MatLab HDF tools to extract the data from the files provided by NASA as two destriped tiff images and two xml files containing the positions of the satellites and the lattice points in J2000.00ECEF Cartesian coordinates as well as the lattice points in image coordinates.

CONCLUSION AND PERSPECTIVES
For the need of this study and because the Silcast software is a closed source, "black box"-type software, we have developed an open source processing chain that takes ASTER L1A scenes in HDF4 format, extract the images from the VNIR bands 3N and 3B as well as the necessary meta-data and computes a DEM and ortho-image of the scene.From this baseline, we can conduct fine analysis of the raw data and apply corrections to it.
The cross-track jitter correction method presented here offers a significant improvement in the quality of the correlation leading to significantly reduced noise, more accurate local geometry and facilitation of the terrain analysis compared to both the products obtained without corrections and the products obtained through previous ASTER DEM generation software.The ground sample distance (GSD) is also halved compared to the ASTER DMO product (from 30m to 15m).Both factors also improve the quality of the ortho-images that can be computed with the new DEMs.
However, as seen in Figure 8 (bottom) where the scene should be flat (apart from the small pieces of land on the sides), the along-track jitter still creates significant distortion in the DEM in the form of a low-frequency wave of ca 10m amplitude and a high frequency wave of ca 2.5m amplitude.If this wave is visible on a flat surface, or in comparison with ground truth, it is more challenging to identify and correct it in more standard data containing variations in terrain.The somehow regular frequencies of the signals combined with their noticeable amplitude are however encouraging future work that might lead to their correction.

Figure 1 :
Figure 1: False color composite of the nadir bands of the "Sea Ice" scene used for most figures of this paper.

Figure 3 :
Figure 3: Extract of the VNIR band 3N of an ASTER scene.Top : Original image.Bottom : After destriping.

Figure 6 :
Figure 6: Measure of the combined across-track vibration visible in the backlooking image through bi-dimensional correlation (scene "Sea Ice") in image space coordinates.

For
the jitter estimations and corrections, an implementation in the free open source MicMac photogrammetric library (developed at the French National Institute of Geographic and Forest Information -IGN -[Pierrot-Deseilligny et al., 2015] and [Pierrot-Deseilligny, 2015]) was conducted, hence offering a ready to use software to produce DEMs and ortho-images from Level 1A ASTER products.

Figure 9 :
Figure 9: Extract of a shading of a scene in central Norway.Top : DEM computed from original images.Bottom : DEM computed after corrections.We can easily see the reduction of the noise, revealing the geomorphological features of the terrain.