ASSESSMENT OF THE MAPPING POTENTIAL OF PLÉIADES STEREO AND TRIPLET DATA

The Pléiades satellites provide very high resolution optical data at a swath width of 20 km and a ground sampling distance of about 0.7 m at nadir direction. The sensors are remarkable agile as their pointing angle can be changed in a range of ± 47 degrees. Thus, they are able to collect three images in one over flight representing tri-stereo data. In the presented work the mapping potential of Pléiades stereo and tri-stereo data is assessed in detail. The assessment is performed on two test sites and contains discussions on 2D initial geo-location accuracy, sensor model optimization, 3D geo-location accuracy, and a novel workflow for dense reconstruction of digital surface models (DSMs). The main outcomes are that the sensor accuracy is within the range as defined by Astrium, however a sensor model optimization is obligatory when it comes to highly accurate 3D mapping. The derived DSMs show a high level of detail thus enabling varying applications on a large scale, like change detection or forest assessment.


INTRODUCTION
The Pléiades satellite system is a dual system comprising the two identical satellites Pléiades-1A and Pléiades-1B, which have been launched in December 2011 and in December 2012, respectively, both providing very high resolution (VHR) image data.The satellites operate in the same orbit with an offset of 180 degrees to offer a daily revisit capacity.They are supplied with a remarkable agility, as the pointing angles can be triggered in a range of ± 47 degrees (standard mode ± 30 degrees).The sensors are capable of acquiring a panchromatic band (470-830 nm) with 0.7 m ground sampling distance (GSD) at nadir and four multi-spectral bands (blue: 430-550 nm; green: 500-620 nm; red: 590-710 nm; near infrared: 740-940 nm) with 2.8 m GSD.The swath width of Pléiades image data is 20 km on ground (Astrium, 2012;Poli et al., 2013;Gleyzes et al., 2013).
The satellite agility is of importance with respect to 3D data extraction from Pléiades image data.Similar to the present very high resolution missions of Ikonos, Quickbird, Worldview or GeoEye, stereo data can be acquired during one over flight (single pass) through appropriate forward and backward arrangement of the sensor.A significant innovation and advantage of Pléiades, however, is provided through the capability to acquire even three images for an area, taken from the same orbit at along track forward-, nadir-and backwardview of the sensor and through the possibility of an across-track swipe.Such image triplets are also denoted as tri-stereo data sets.First assessment of the benefit of such arrangements was made by the authors with respect to image triplets composed from multi-sensor and multi-temporal acquisitions, respectively (Raggam, 2006).
Obviously the swath length that can be acquired for a predefined base-to-height (B/H) ratio is limited and shorter for a triplet than for a stereo acquisition (e.g. 60 km vs. 175 km at a B/H-ratio of 0.4).However, large B/H-ratios tentatively in the order of 0.7 to 1.0 should be aspired with regard to 3D mapping applications.Then swaths with a length of up to about 200 km can be well acquired in tri-stereo mode.
In the context of (3D) mapping applications the geo-location accuracy, which is inherent to or can be achieved for the image data being used, is of importance.For Pléiades products an apriori geo-location accuracy in the order of 8.5 m CE90 only was proposed.Improvement on the one hand was expected from a recalibration of the Pléiades-1A satellite, and can be on the other hand achieved through sensor model optimization using ground control points (GCPs).This paper deals with a detailed accuracy assessment of panchromatic Pléiades image data for two test sites, one being covered by a tri-stereo and the other by a pure stereo acquisition.The study contains an assessment of the 2D geolocation accuracy, employing the initial sensor models as well as optimized sensor models, and an assessment of the 3D mapping accuracy involving both stereo as well as tri-stereo data sets.Another focus of this paper deals with the extraction of dense digital surface models (DSMs) from Pléiades image data in theory and practice.DSMs extracted from stereo as well as tri-stereo data sets are presented, and first validation results achieved from a comparative analysis with reference LiDAR DSMs are included.The whole assessment is performed with a software package designed and implemented by the authors at JOANNEUM RESEARCH (JR).

Test Sites and Pléiades Image Data Sets
Two test sites were considered for this study, located in Trento (Italy) and Innsbruck (Austria), both of them showing mountainous topographic characteristics.The Trento test site covers rural as well as mountainous terrain, the ellipsoidal heights ranging from 175 to 1550 meters a.s.l. and spans over an area of 220 km 2 .A Pléiades-1A triplet was acquired, which however is far from optimal for 3D reconstruction, since the first two images have a very small intersection angle of only 5.5°, while the second stereo pair has a huge intersection angle of 27.3°.
For the Innsbruck test site only a stereo acquisition from the Pléiades-1A satellite was available.This set was taken after the launch of Pléiades-1B and thus after the recalibration of the Pléiades-1A satellite.In theory the 2D geo-location accuracy for this data set thus should have improved.The test site covers urban, rural and mountainous terrain, with ellipsoidal heights ranging from 560 to 2750 meters a.s.l. and spans over an area of 435 km 2 .
All images where acquired with platform PHR 1A, in acquisition mode PX and spectral processing PA.Acquisition parameters of the image data, like along, across and overall incidence angle are summarized in Table 1.In this study only the VHR panchromatic bands are used, which originally are acquired with a GSD of 0.7 to 0.77 m, but always are delivered with an upsampled GSD of 0.5 m (2.0 m for multi-spectral bands).Figure 1 shows ortho-rectified products generated from the multi-spectral Pléiades IMG1 data, visualized as colorinfrared (CIR) illustration for Trento and as true color RGB for Innsbruck.Figure 2

Reference Data
For both test sites, airborne LiDAR DSMs were available as reference data sets, collected in 2006 and 2007 at 1.3 points/m 2 for the Trento and at 2.0 points/m 2 for the Innsbruck test site, respectively.The LiDAR DSMs are available as raster data in UTM 32 North map projection and WGS 84 datum at 1 m GSD with ellipsoidal heights and could be further used to measure GCPs.3D discontinuities in the LiDAR data can be used to define and manually measure ground coordinates with a sufficiently high accuracy, i.e. with an estimated measurement error of less than 1 m.The purpose of GCP measurement was • to validate the quality of the initial Pléiades sensor models as provided along with the image data, • to optimize the Pléiades sensor models by means of least squares parameter adjustment procedures and • to assess the 3D mapping potential of Pléiades stereo and triplet acquisitions.
In order to get a representative feedback on these issues, a fairly high number GCPs, homogeneously distributed in planimetry and height, was measured for both test sites.Targets like road intersections, corners of houses, water bodies or field boundaries, served as control point candidates.Thus, 21 and 30 points, respectively, were acquired for Trento and Innsbruck.

SENSOR MODEL VALIDATION AND OPTIMIZATION
Validation of the geo-location accuracy can be made via residuals of GCPs, which can be determined by means of the available sensor model.Therefore, in general a backward transformation, i.e. transformation of map to image coordinates, of the GCPs is first made.Point residuals in along and across track direction then can be determined as the differences between transformation results and measured pixel coordinates.

Initial Geo-location Accuracy
The geo-location accuracy inherent to Pléiades panchromatic imagery is reported to be 8.5 m CE90 at nadir direction and 10.5 m CE90 within 30 degrees off-nadir, when applying the provided rational polynomial coefficient (RPC) model (Astrium, 2012).
For validation, these RPCs were used to determine "initial" residuals for all image data sets, representing the initial geolocation accuracy.Mean values as well as standard deviation of across (x) and along (y) track pixel residuals are summarized in Table 2. Considering the 0.5 m GSD of the Pléiades image data, the geo-location accuracy of the Trento data set is up to about 18 pixels or 9 meters, which fairly well corresponds to the stated 8.5 m CE90.For the Innsbruck data set, a geo-location accuracy of up to 10 pixels, corresponding to about 5 meters, is achieved.Sensor recalibration of Pléiades-1A could be the reason for this improvement, which on the other hand is not sufficient to fulfil precision mapping requirements.However, it should be noted that the standard deviation residual values in along-track direction (y) have remarkably improved for the Innsbruck data set.This is of importance with respect to 3D mapping, as height accuracy in general is more severely affected by along-track pixel errors.
It further should be noted that the initial mean error values are diverse for the two test sites with respect to sign and order of magnitude.Thus the a-priori geo-location accuracy of Pléiades cannot be calibrated globally, e.g. by an appropriate pixel shift.

Parameter Adjustment
A need to optimize the Pléiades sensor models is evident in order to achieve improved geo-location accuracy, e.g. in the sub-pixel/sub-meter range.Therefore, in general least squares parameter adjustment procedures based on GCPs are used.With respect to the rational polynomial models being used for these test data sets, the intention is to apply optimization to as few parameters as reasonable in order to reduce need for GCPs and to avoid drifting or even oscillating of the RPCs.With respect to follow-on extraction of 3D information therefore optimization of only the linear terms of the RPC nominators was performed, involving all measured GCPs.While only 4 GCPs would be required as a minimum, over-determination is reasonable anyway in order to mitigate the impact of measurement errors.
The RMS values of across/along track residuals (x/y) as well as residuals length (|xy|) as resulting after optimization of the sensor models, i.e.RPCs, are also given in Table 3.Now these values to the majority are in the sub-pixel range (0.6 to 0.9 pixels) and show a widely homogeneous behaviour in along and across track direction.According to the nominal GSD of 0.5 meters these pixel values correspond to geo-location errors in the range of 0.3 to 0.45 meters on ground.(Dial and Grodecki, 2002).An even more simplified approach is suggested by (Fraser and Hanley, 2003)

ASSESSMENT OF 3D MAPPING ACCURACY
Using point measurements made in the stereo or triplet image data, the corresponding 3D ground coordinates can be calculated through a least squares spatial point intersection.
Corresponding point residuals can then be determined in comparison to the reference ground coordinate measurements.This analysis was made for the different stereo pairs and the triplet included in the test data sets based upon utilization of initial as well as optimized sensor models.The RMS values of the 3D point residuals which were achieved for East, North and height are summarized in First, it is obvious, that the initial accuracy of the rational polynomial models yields 3D displacements in the order of several meters, which is clearly beyond aspired precision.Next, the values given in the table show, that the stereo intersection angle ∆θ -as an equivalent for the base-to-height ratio -has a predominant impact onto the 3D geo-location accuracy.Higher accuracies are achieved for image pairs covering larger stereo intersection angles and vice versa.Thus, the Trento stereo pair with the small intersection angle of 5.5 degrees yields significant worse 3D geo-location than all other stereo (and tristereo) constellations.
For With respect to the 3D mapping approach being applied to the Pléiades data sets it has to be emphasized, that optimization of the rational polynomial models is indispensable.In the respective workflow (see chapter 5) accurate sensor models are a demand to generate strict epipolar image pairs (cf.Gutjahr et al., 2014) with mutually corresponding image lines.This is a pre-requisite to apply 1D matching algorithms, like the semiglobal matching approach (Hirschmüller et al., 2008).

SURFACE MAPPING PROCEDURES
In the following the 3D mapping procedures of JR, which have been applied to the Pléiades test data, are described.

Stereo-processing Workflow
Starting from a pair of images to be used for the generation of a digital surface model (DSM), the following procedures are applied: • Epipolar rectification of both stereo images based on the optimized sensor models, such that a pre-defined point in the reference image can be found along a horizontal line in the search image.While the concept of epipolar geometry was first realized for perspective images, appropriate implementations were further made for Pléiades-like pushbroom geometries (Wang et al., 2011) and for SAR geometries (Gutjahr et al., 2014).In either case the generation of epipolar image pairs is based on the underlying sensor geometries and relies on accurate sensor models in order to achieve pairs with strictly corresponding image lines.• Image matching in order to find correspondences of reference image pixels in the search image.In case of epipolar images highly sophisticated image matching algorithms can be applied which utilize the fact that the search is restricted to one dimension and thus local and also global optimization methods can achieve highly accurate results at fast runtimes.A good overview and benchmark of such algorithms is given in (Scharstein and Szeliski, 2002).For this work a custom-tailored version of the semi-global matching (SGM) algorithm, which was introduced in (Hirschmüller, 2008) was used.Disparity predictions can be calculated based on the sensor models and coarse elevation data (e.g. the SRTM DSM), such that the search range of the SGM algorithm can be limited.
The cost function compares two image patches and is defined as the Hamming distance of the two Census transforms within 9x9 pixel windows on the epipolar images.Within the semi-global optimization an adaptive penalty is used to preserve 3D breaklines.The image matching procedure yields two dense disparity maps pointing from the reference to the search image and vice versa.
• Spatial point intersection, in order to calculate ground coordinates from the corresponding image pixels retrieved from image matching by means of a least squares approach.To some extent, unreliable matching results can be identified and rejected in this step.This procedure results in a "cloud" of 3D points, irregularly distributed on ground.• DSM resampling, i.e. interpolation of a regular raster of height values from these 3D points.Remaining gaps may be filled using an appropriate interpolation mechanism.• DSM fusion to combine the information of forward and backward matching.Here the method described in (Rumpler et al., 2013) is applied, which takes all height measurements within a 3x3 neighbourhood and extracts the mode of this height distribution.Thus, local height errors can rather be eliminated than by a straightforward median based approach.

Utilization of Multiple Matching Results
A known trade-off in stereo data processing is that large intersection angles should be inherent to the stereo pair in order to imply higher geometric robustness and geo-location accuracy, but that on the other hand the images should be similar, i.e. acquired under small stereo intersection angles, in order to assure successful matching.Utilization of multiple stereo pairs and fusion of multiple matching results can be applied to consider both requirements and thus increase the 3D reconstruction quality.
In this study each stereo pair was processed independently, yielding two DSMs according to the workflow described in Section 5.1.Fusion is then also done according to this workflow, but integrates multiple stereo-derived DSMs retrieved from individual stereo pairs of the triplet.Thus, for the Trento test site all possible stereo pairs within a triplet, namely IMG1-IMG2, IMG2-IMG3 and IMG1-IMG3, were processed in both feasible stereo constellations, finally yielding six DSMs, which are then fused to one final DSM.
One advantage of this workflow is that each stereo pair can be processed individually on multiple computers.This is of particular importance with respect to image matching, which represents the workflow bottleneck in terms of runtime.Best feasible efficiency thus is achieved, because only the final DSM fusion step has to simultaneously access the multi-stereo processing results.It should be noted that this workflow can be used to process any other kinds of images as well, like in particular aerial images, e.g.UltraCam images (Leberl et al., 2003).

Digital Surface Models
For both test sites dense DSMs were extracted using the Pléiades data and employing the workflow presented above.While for the Innsbruck test site only a stereo pair was available, a tri-stereo data set could be processed for the Trento test site.For visual comparison and analysis, subsets of both test sites are presented below.
Figure 3 shows a subset of the Trento test site, which covers a hospital and its surrounding area of 360 x 360 m 2 .A Pléiades ortho-image, the corresponding LiDAR DSM, the stereoderived DSMs as well as the triplet-derived DSM are illustrated.
The terrain heights of the DSMs are scaled between 240 (black) and 300 m (white) a.s.l.Infrastructural changes due to ongoing construction activities since the LiDAR acquisition can be seen.E.g. there is a new helipad (center-left) and a parking lot (center-down), which previously used to be a park.The stereoderived DSMs in general are affected by occlusion areas, which use to increase with increasing stereo intersection angle, and which almost vanish in the triplet-derived DSM.The DSM derived from the IMG1-IMG2 Pléiades stereo pair looks different in comparison to the others and seems to be more similar w.r.t. the LiDAR model.This is due to the small stereo intersection angle, which implies higher image similarity and thus a higher performance of the stereo matching with reduced occlusion areas.However, due to this weak geometric disposition the height accuracy of this DSM is worse than for the other DSMs.
Exemplary quantitative analysis was made through comparison of LiDAR and stereo/triplet-derived DSMs.Due to the large temporal gap between LiDAR and Pléiades data acquisition only selected areas were analysed which are not affected by temporal change due to construction, vegetation growth or cloud cover.The results are summarized in Table 7.The standard deviations of the height differences confirm the interdependence of stereo intersection angles and 3D mapping accuracy, with worst accuracy achieved for the stereo pair  Adequate to the Trento test site, a quantitative accuracy assessment was made for a selected area of this test site, which could be supposed to be free of temporal changes.The results are included in

CONCLUSION AND OUTLOOK
The presented assessment of the mapping potential of Pléiades stereo and tri-stereo data revealed that the 2D geo-location accuracy of the initial sensor model is within the range of 8.5 m CE90 as stated by Astrium.For accurate 3D mapping procedures the sensor models have to be optimized employing ground control points.After such an optimization the 2D geolocation accuracy is in the range of 1 pixel and the 3D geolocation accuracy is in the range of 0.5 m in planimetry and in the range of 1 m in height.Such accuracy can be seen to be feasible at the best with respect to 3D models extracted from Pléiades data.
In the study also a tri-stereo Pléiades data set was incorporated within a 3D mapping workflow.It could be in general shown that the increased over-determination resulting from multiple stereo matching results yields higher quality surface models.
Due to the unfavourable intersection angle of this particular tristereo set, however, the potential of the triplet-based DSM generation could not be fully exploited.In future, tri-stereo sets should be acquired, which comprise a nadir image as well as forward and backward images, which are acquired at about the same forward and backward look angles.
When thinking of processing large amounts of Pléiades triplets for 3D mapping the current bottleneck is the GCP measurement, if to be done manually, which is required to optimize the Pléiades sensor models.Thus, integration of automatic procedures for GCP acquisition would be the next logical step in order to realize a fully automatic pipeline.Our idea is to define a GCP-chip database such that control points could be retrieved therefrom automatically via area based image matching.

Figure 3 .
Figure 3. Detail view of a hospital in test site Trento.

Figure 4 Figure 4 .
Figure 4 shows a subset of the Innsbruck test site, which covers an urban region of 1950 x 1000 m 2 .It illustrates an orthoimage, a LiDAR DSM, a stereo DSM generated by the commercial software Geomatica 2013 by PCI Geomatics 1 (denoted PCI-DSM), and the stereo-derived DSM as generated by the workflow sketched in Section 5.1 (JR-DSM).The DSM height values are scaled within 620 to 680 m a.s.l.When comparing LiDAR and Pléiades DSMs, temporal changes again become visible, like a new residential area in the upper-left image area.It is also obvious that the JR workflow preserves 3D breaklines better than the PCI workflow.For instance large buildings are decently reconstructed while they are missing in the DSM retrieved by PCI.These differences very likely result from the different stereo matching algorithms being used.For a partly forested region of test site Innsbruck with 550 x 550 m 2 Figure5shows a Pléiades ortho-image as well as a the height differences between the LiDAR and the Pléiades DSM, scaled from -25m to +25 m.Thus, bright areas indicate forest

Table 1 .
visualizes the densely reconstructed DSMs derived by our proposed method.It is clearly visible that both test sites cover difficult terrain from urban areas to steep mountainous zones.Acquisition parameters of the Pléiades data.

Table 2 .
A-priori geo-location errors of Pléiades image data.

Table 3
. RMS residuals after parameter adjustment.The decision to adjust the linear nominator coefficients is based upon previous experiments involving different adjustment settings, which have shown, that sub-pixel accuracy thereby can be achieved.It also corresponds to findings made by

Table 4 .
, recommending to apply only a shift to the RPCs and thus to reduce the need for GCPs even more.Related tests were made for the image IMG1 of the Innsbruck test data, using a different number of GCPs in a comparative assessment of shift versus nominator coefficient adjustment.Remainder points were used as independent check points.The results of this analysis, including mean and standard deviation values of check point residuals, are summarized inThe table shows, that utilization of an absolute minimum of GCPs, i.e. 1 for shift and 4 for nominator coefficient optimization, yields systematic errors for both scenarios, as expressed by the mean residual values.Appropriate overdetermination, e.g.utilization of 10 GCPs in this assessment, may reduce such systematic errors significantly.The numbers also confirm that removal of only a shift then might be sufficient to achieve reliable accuracy.

Table 5
. 3D RMS residuals achieved from initial as well as optimized sensor models (given in meters).

Table 7 .
Table 7 and show a slightly better accuracy than the Trento data sets.Pléiades ortho-rectified LiDAR minus Pléiades DSM Figure 5. Detail view of a partly forested area of test site Innsbruck.Statistics of height differences between LiDAR and Pléiades DSMs on non-forest regions.