DSM BASED ORIENTATION OF LARGE STEREO SATELLITE IMAGE BLOCKS

High resolution stereo satellite imagery is well suited for the creation of digital surface models (DSM). A system for highly automated and operational DSM and orthoimage generation based on CARTOSAT-1 imagery is presented, with emphasis on fully automated georeferencing. The proposed system processes level-1 stereo scenes using the rational polynomial coefficients (RPC) universal sensor model. The RPC are derived from orbit and attitude information and have a much lower accuracy than the ground resolution of approximately 2.5 m. In order to use the images for orthorectification or DSM generation, an affine RPC correction is required. In this paper, GCP are automatically derived from lower resolution reference datasets (Landsat ETM+ Geocover and SRTM DSM). The traditional method of collecting the lateral position from a reference image and interpolating the corresponding height from the DEM ignores the higher lateral accuracy of the SRTM dataset. Our method avoids this drawback by using a RPC correction based on DSM alignment, resulting in improved geolocation of both DSM and ortho images. Scene based method and a bundle block adjustment based correction are developed and evaluated for a test site covering the nothern part of Italy, for which 405 Cartosat-1 Stereopairs are available. Both methods are tested against independent ground truth. Checks against this ground truth indicate a lateral error of 10 meters.


INTRODUCTION
In May 2005 India launched its IRS-P5 satellite with the CARTO-SAT-1 instrument which is a dual-optics 2-line along-track stereoscopic pushbroom scanner with a stereo angle of 31 • and the resolution of 2.5 m.The CARTOSAT-1 high resolution stereo satellite imagery is well suited for the creation of digital surface models (DSM).In this paper, a system for highly automated DSM generation based on CARTOSAT-1 stereo scenes is presented.More details about CARTOSAT-1 are given in (Srivastava et al., 2007).
CARTOSAT-1 stereo scenes are provided with rational polynomial functions (RPC) sensor model (Grodecki et al., 2004), derived from orbit and attitude information.The RPC have a much lower accuracy than the ground resolution of approximately 2.5 m.Subpixel accurate ground control points (GCP) have been used in previous studies to estimate bias or affine RPC correction parameters required for high quality geolocation of HRSI images (Lehner et al., 2007).Such highly accurate GCP are usually derived from a DGPS ground survey or high resolution orthoimages and digital elevation models.For large scale and continent wide processing, establishing a highly accurate GCP database with the density required for processing the relatively small CARTOSAT-1 scenes (900 km 2 ) requires significant resources.
The georeferencing algorithms developed in this paper are part of a larger processor for the generation of a high resolution, European DSM from Cartosat-1 stereo scenes (Uttenthaler et al., 2011).Several tens of thousands stereo scenes need to be processed for a complete coverage.The DSM generation process includes affine RPC correction of all scenes, dense stereo matching using Semiglobal Matching (d'Angelo, 2010), reprojection and merging of the DSMs derived from each stereo pair, closing remaining holes in the DSM with delta surface fill (Grohman et al., 2006) and interpolation of remaining voids and orthorectifaction of the near nadir images.
We propose the use of widely available lower resolution satellite data, such as the Landsat ETM+ and SRTM DSM datasets as a reference for RPC correction.The accuracy of these datasets is low compared to the high resolution CARTOSAT-1 images.The absolute lateral error of ETM+ Geocover is 50 m (LE90), the absolute lateral error of SRTM is between 7.2 m and 12.6 m (LE90, depending on the continent), with an absolute height error of 4.7 m to 9.8 m (Rodriguez et al., 2005).The traditional method of collecting lateral position from a reference image and interpolating the corresponding height from the DEM ignores the higher lateral accuracy of the SRTM dataset.Our method avoids this drawback by using a RPC correction based on DSM alignment, resulting in improved geolocation of both generated DSM and ortho images.
The first part of the paper describes a scene based and a bundle block adjustment based georeferencing.The second part evaluates the proposed orientation algorithms on a large test area with 414 Cartosat-1 stereo scenes.

GEOREFERENCING USING DSM AS MAIN GROUND CONTROL
The use of height models, such as DSM as only, or major ground control for the orientation of images has been evaluated by various previous authors.Bundle block adjustment of aerial frame imagery with DSM as ground control has been performed by (Strunz, 1993).Here, the differences between the height of the object coordinates and the corresponding bilinarly interpolated height from the DSM are used as additional observation in the bundle adjustment.A similar approach was used to orient Mars Express HRSC imagery (M.Spiegel, 2006).In these works, only height differences between the DSM surface and the object points are used.A refined method by (Jaw, 2000) considers the euclidean distance between a surface spanned by the image tie points and the surface points.A similar approach has been used by (Akca, 2007) for point cloud alignment with Helmert transfor-  (Grodecki and Dial, 2003), it describes the relationship between object coordinates φ,λ,h and image coordinates r,c: The RP C col and RP Crow functions are determined by the satellite operator using the measured attitude and ephemeris data and the physical sensor model.The location error of the Cartosat-1 RPCs is approximately 100-250 m.These errors can be corrected by using image space RPC correction parameters (Grodecki and Dial, 2003), which capture the effect of roll, pitch error and ephemeris errors.Previous studies have shown that Cartosat-1 Orthokit imagery requires the use of a fully affine RPC correction (Lehner et al., 2007): The a and b coefficients can be determined using a set of accurate, well distributed GCPs.A bundle block adjustment can be used to estimate RPC cofficients for more than one image at a time.

Preprocessing
The orientation procedure starts with matching tie points between the stereo pair and a transfer matching to the reference image (for approximate GCPs).Hierarchical intensity based matching is used for matching the stereo pairs and the reference image.It consists of two major steps, hierarchical matching to derive highly accurate tie points and dense, epipolar based stereo matching.
The initial matching step uses a resolution pyramid to cope with large stereo image distortions originating in carrier movement and terrain.Large local parallaxes can be handled without knowledge of exterior orientation.The selection of pattern windows is based on the Förstner (Förstner and Gülch, 1987) interest operator which is applied to one of the stereo partners.For selection of search areas in the other stereo image, local affine transformations are estimated based on already available tie points in the neighborhood (normally from a coarser level of the image pyramid).Tie points with an accuracy of one pixel are located via the maximum of the normalized correlation coefficients computed by sliding the pattern area all over the search area.These approximate tie point coordinates are refined to subpixel accuracy by local least squares matching (LSM).The number of points found and their final (subpixel) accuracy achieved depend mainly on image similarity and decrease with increasing stereo angles or time gaps between imaging.Strict thresholds on correlation coefficient and bidirectional matching differences are used to select reliable and highly accurate stereo tie points.
After high quality tie points within the stereo pairs have been found, they are transferred to the lower resolution reference image using local least squares matching.Due to large acquisition time differences between the ETM+ and the CARTOSAT-1 imagery, only a small fraction of the stereo tie points can be matched reliably between all three images.Three-dimensional GCP are found by bilinear interpolation of the SRTM DSM at the geographic positions found in the ETM+ reference scene.These GCPs are called Image GCP in the following sections.

Stereo model alignment
The first method for RPC correction is a two step process.First, the Image GCPs are used to estimate bias RPC correction parameters.This results in good forward intersection behavior, but the lateral and vertical accuracy is still limited by the ETM+ reference.The higher lateral accuracy of the SRTM dataset is used in a second RPC correction step.A 3D point cloud is calculated by forward intersection of approximately 1 million stereo tie points.This point cloud is aligned to the SRTM DSM using 3D least squares matching.A subset of the aligned point cloud is used as GCPs for an estimation of the final affine RPC correction parameters.
It is assumed that the height zi of a point Pi located at (xi, yi, zi) equals the reference DSM height hD(xi, yi) at the corresponding position (xi, yi): A 3D affine transformation is used to align the initial stereo point cloud to the SRTM DSM: where pi = (xi, yi, zi, 1) is the original point, A is a 3×4 matrix, and pti = (xti, yti, zti) is the transformed point.
The affine transformation matrix A is estimated using an iterative least squares algorithm.Using Eq. 3 and 4, the following observation equation is obtained: Since the model is non-linear, the solution is obtained iteratively.An identity transform is used as initial approximation, since the stereo points are not far from the reference.It is likely that the stereo point cloud, and to a smaller extend the DSM contains outliers or changes due to the 10 year acquisition difference of the data, which cannot be handled by a standard least mean squares algorithm.After the initial estimation, points with a residual larger than 3 times the standard deviation are removed and a new transformation is estimated.This procedure is repeated until less than 0.3 % outliers are detected and the squared sum of the outlier residuals accounts for less than 5 % of the squared sum of all residuals.
As CARTOSAT-1 requires an affine RPC correction, well distributed GCPs are essential for high quality results.The 3D alignment process works well in areas where significant terrain is available throughout the whole scene.In practice, this methods works very well for images with enough terrain (Uttenthaler et al., 2011).Scenes that contain no usable relief for DSM alignment can only be corrected by using the Image GCP.Problems also occur for hilly scenes that contain large flat areas, as the 3D matching does not provide lateral constraints in flat areas.This is similar to using an uneven distribution of ground control points.A method that can bride some flat areas is required to avoid the need for additional GCP acquisitions.

Bundle block adjustment
An improved method should be able to handle multiple stereo pairs at once in an single estimation step.Furthermore, flat areas should be bridged either by constraints from neighboring images, or with additional ground International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XXXIX-B1, 2012 XXII ISPRS Congress, 25 August -01 September 2012, Melbourne, Australia control points for flat areas.The two step method also lacks a direct error modeling, and does not necessarily yield the optimal parameter estimates.
A bundle block adjustment of larger blocks is thus required for the seamless processing of larger blocks.We follow the RPC bundle block adjustment of (Grodecki and Dial, 2003), which is based the following observation equations: where φi, λi, hi are the object coordinates, a,b are the RPC correction parameters, ri,ci are the measured image coordinates in image j.If the ground coordinates (φi, λi, hi) for some points are known, for example by GPS measurement or from matching against a reference image, an additional observations is added: Prior information about the rpc correction terms can be added, if available: An extension to (Grodecki et al., 2004) is the DEM height difference observation term, similar to the approch by (Strunz, 1993): where D(φ, λ) is the bilinarly interpolated height in the reference DEM at position (φ, λ).
Using the observations equations ( 6)-( 9), and corresponding weights W , based on prior knowledge, an iterative least squares estimation is used to estimate tie point object coordinates and the R correction parameters, see (Grodecki et al., 2004) for details.In classical bundle block adjustments, only a small number of well distributed tie points are required for each image.If only DEM observations are used, many well distributed tie points are required.These tie points should capture the main relief in the scene to add horizontal constraints to the adjustment.This requires the use of computationally efficient estimation algorithms, as millions of unknowns need to be estimated for larger blocks.
Our implementation first eliminates the object coordinates from the normal equations, using the Schur complement (Triggs et al., 2000) and uses sparse cholesky factorisation techniques to solve for the unknowns.This is efficient, as no inversion of the normal matrix is required, but currently does not yield a theoretical covariance of the estimated object point coordinates.

EVALUATION
The above process is evaluated using a block of 414 Cartosat-1 Stereopairs covering the northern part of Italy.This is a challenging test area for the DSM based georeferencing, as the flat Po plain covers a major part of the test area.Each stereo pair covers an area of approximately 28 by 30 km.An overview of the test area and the location of the scenes is shown in Fig. 1.Multiple stereo pairs are available for many areas.The scenes were acquired in the years 2008-2010, and are mostly cloud free, but some scenes contain a cloud cover of up to 20%.Out of the 828  images, 18 were very cloudy, contained only water areas or were otherwise unsuitable for matching.1927389 well distributed tie points have been matched between the remaining 810 images using SIFT (Lowe, 2004) and local least squares matching.The large number of tie points is required for the DSM alignment term, Eq. ( 9).
The SRTM-C Band DSM,version 2, unfilled and edited by NGA is used as primary reference.

Stereo pair alignment
In a first step, the performance of single stereo pair alignment, as described in Section 2.2 is evaluated against the Euro-Maps 2D mosaic.Approximately 1 million tie points from the dense

Dependency of bundle block adjustment accuracy on the number of tie points
While the number of tie points is not a concern for the stereo pair alignment, using 1 million points per stereo pair would result in approximately 1.2 billion unknowns for the adjustment of the complete block.This cannot be handled by the block adjustment program, and thus a reduced number of tie points is required.For every test area, two main kinds of tie points are provided for the block adjustment: tie points derived from SIFT image matching and from LLSM image matching.As the DSM alignment depends on well distributed points at various slopes and aspects, a large number of tie points is required.By varying the number of used tie points, the relationship between the number of tie points and the accuracy and stability of the adjustment is evaluated.Different numbers of SIFT and LSM tie points are derived from image matching by thinning with a grid.Only a single tie point is selected from each grid cell.We select Cartosat-1 stereo pairs overlapping with the RealItaly test areas shown in Fig. 2, and compute the planimetric RMSE after adjustment computing a separate block adjustment for each overlapping scene.As an example, the results for SIFT tie points are shown in Fig. 4 and LSM tie points in Fig. 5.It can be seen that large fluctuations appear when a low number of tie points is used.The results stabilize somewhat when more than 1000 tie points are used.We did perfom these tests for all 15 RealItaly test areas, and the LSM points seem to result in less fluctuations than using the SIFT tie points.It should be kept in mind that the resolution of the SRTM DSM is only 90m, and most results are thus stable to 1/10 of the reference data resolution.This agrees with the accuracy achievable by classical local least square matching, which is based on a similar estimation procedure.We thus strive to use more than 1000 tie points per image pair.

Bundle block adjustment with DSM as ground control
After analysis of single stereo pairs, the bundle block adjustment is performed for all 810 images.1.927.389tie points and SRTM as a reference were used in the adjustment, resulting in a least squares problem with 5.787.027unknowns and 10.158.271observations.Note that no GCPs are used in this experiment, and the process is fully automatic.The mean errors for checkpoints extracted from the Euro-Maps 2D mosaic are shown in Fig. 6.Except for the Adria coast on the left, the errors are reasonably small, considering that the accuracy of the Euro-Maps 2D mosaic is specified with CE90 15 m, and SRTM with CE90 10 m.A trend is seen when evaluating the checkpoints extracted from the 15 RealVista test sites, cf.Fig. 7.The flat Adria coast does not provide lateral constraints for the right side of the block, thus resulting in larger errors, as the DSM term does not provide lateral constraints.It is thus useful to stabilize these areas with additional ground control points.

Bundle block adjustment with DSM and GCPs
In order to stabilize the block at the eastern edge, some GCPs have to be used in the bundle block adjustment.As the Landsat Mosaic has errors larger than 20 m in this area, it cannot be used as additional ground control.Checkpoints derived from Euro-Maps 2D are used as GCPs for two scenes in the upper right and lower right corner of the block.The results are shown in Fig. 8.The large errors on the right edge of the block, as visible in Fig. 6 have been somewhat reduced, but are still significantly larger than in the remaining block.The bundle block results show a wavy pattern, which seems to be correlated with the Cartosat-1 orbits, and which is not visible in the scene based results, cf.Fig. 3.The reason for this wavy pattern is still unclear, it might be related to the lower number of tie points used in the bundle block adjustment (approx.5000 per stereo pair), when compared to the single stereo pair alignment (1 million tie points per pair).Also, the a priori weights for tiepoint and DEM observations also play an important role.The waves can be further reduced by stronger weighting of the DEM measurements, at the expense of slightly larger image residuals.The block is nearly unconnected at latitude 10 • , as no cloud free pass was available for this area.The neighboring scenes contain some clouds and barely touch each other, leading to problems for the tie point matching, and a weak connectivity, and thus higher errors.A summary of the errors for each adjustment scenario is given in Table 2.It is important to note that the Euro-Maps 2D mosaic and the SRTM DSM have similar absolute accuracy, cf.Table 1, thus Euro-Maps 2D is not the best ground truth dataset for this evaluation.While the stereo pair alignment results in a consistent error, it fails for the 101 scenes in flat terrain.The block adjustment succeded in orienting the whole block, but shows slightly higher errors than the scene based method.
TODO: Add errors from aerial images.

CONCLUSIONS
Two methods for the orientation of high resolution stereo satellite imagery using a DSM as major ground control have been developed in this paper.This approach avoids or significantly reduces the need for manually measured GCP, allowing a highly automated orientation of continent wide Cartosat-1 coverage.The first method is a two step approach for the correction of single stereo scenes, which has been shown to work very well for hilly and mountainous areas, but fails in very flat regions.To bridge these problematic areas, a bundle block adjustment, containing a DEM pseudo observation is used to avoid these problems.A flexible and efficient adjustment software for handling large blocks with millions of tie points has been developed, and is used for the evaluation of a block of 405 Cartosat-1 stereo scenes.The test area is the Po plain with the surrounding mountainious regions in northern Italy, a challenging test area due to large flat areas.Using only DEM information results in sufficient accuracy in most of the block, except for large errors on the flat Adria coast, which can be reduced by adding GCPs in the upper right and lower right corners of the block.
The current results indicate, that some GCPs are needed for a stable block adjustment.Further work will include stability analysis of the block and the use of more GCPs inside the block.Currently, all CARTOSAT-1 stereo pairs are treated individually, even if they originate from the same image strip.The bundle adjustment will be extended with additional pseudo observations for larger strips.

Figure 1 :
Figure 1: North Italy test area.The is covered by 414 Cartosat-1 stereo pairs, indicated by the green squares.

Figure 3 :
Figure 3: Mean errors of the Euro-Maps 2D checkpoints after single stereo pair alignment.Each arrow represents the mean X and Y error of a stereo pair.101 stereo pairs were located in very flat areas and could not be oriented using DSM alignment.

Figure 4 :
Figure 4: Evaluation of lateral accuracy vs number of SIFT tie points.Each line corresponds to a scene overlapping with Test area 7.

Figure 5 :
Figure 5: Evaluation of lateral accuracy vs number of LSM tie points.Each line corresponds to a scene overlapping with test area 7.

Figure 6 :
Figure 6: Evaluation of the bundle block adjustment against checkpoints extracted from the Euro-Maps 2D mosaic.Each arrow represents the mean X and Y error of a stereo pair.RPC corrections for all scenes were computed.

Figure 7 :
Figure 7: TODO: create nicer figure.Evaluation of the bundle block adjustment against checkpoints extracted from the Re-alVista aerial images.Each arrow represents the mean X and Y error of a stereo pair.

Figure 8 :
Figure 8: Evaluation of the bundle block adjustment with DEM and GCPs in the upper and lower right corners (indicated by the black boxes).
Summary of errors against the Euro-Maps 2D.Note that the stereo pair alignment failed for the more demanding scenes with flat terrain.

Table 1 :
Reference and check datasets and their horizontal and vertical accuracy.