DSMS GENERATION FROM COSMO-SKYMED , RADARSAT-2 AND TERRASAR-X IMAGERY ON BEAUPORT ( CANADA ) TEST SITE : EVALUATION AND COMPARISON OF DIFFERENT RADARGRAMMETRIC APPROACHES

This work is focused on the analysis of potentialities of the radargrammetric DSMs generation using high resolution SAR imagery acquired by three different platforms (COSMO-SkyMed, TerraSAR-X and Radarsat-2) with particular attention to geometric orientation models. Two orientation models have been tested in this work: the rigorous Toutin’s model, developed at the Canada Center for Remote Sensing (CCRS) and implemented in the commercial software package PCI Geomatica, and the radargrammetric model developed at University of Rome La Sapienza and implemented in the scientific software SISAR. A full comparison and analysis has been carried out over Beauport test site (Quebec, Canada), where a LIDAR ground truth and a dense set of GNSS CPs (Check points) are available. Moreover, a preliminary comparison between the DSMs extracted, respectively with SISAR and PCI-Geomatica has been performed. The accuracy of the generated DSMs has been evaluated through the scientific software DEMANAL developed by Prof. K. Jacobsen of University of Hannover. As regards orientation models, the results shown that the Toutin’s model accuracy is slightly better than the SISAR one, even if it is important to underline that the SISAR model is computed without using a priori ground truth information. As concern DSMs assessment, the global DSMs accuracy in term of RMSE is around 4 meter and the two radargrammetric approaches show similar performances.


INTRODUCTION
The satellite remote sensing technology offers the opportunity to have continuous observation of Earth's surface for territorial application, with short acquisitions and revisit times, satisfying the demand for monitoring rapid changes in the ground and anthropic activities.In particular, the SAR (Synthetic Aperture Radar) high resolution satellite imagery could offer night-andday and all-weather functionality (clouds, haze and rain penetration), that represents an important advantage for time-series analysis and for rapid mapping.Starting from SAR data, one of the most important derived products are the Digital Surface Models (DSMs/DEMs).These kind of products allow to have a synoptic knowledge of the land morphology, with different level of accuracy and details, depending on the characteristics of the sensor.In details, there are two main techniques to generate DSMs from SAR data: the well-known interferometric approach and the less used but promising radargrammetric one (Hanssen, 2001).
Radargrammetry was first used in the 1950s and it represents an alternative solution able to avoid the classical decorrelation problem affecting the interferometric technique especially over areas with vegetation/forest, and, at present, the importance of the radargrammetric approach is rapidly growing due to the new high resolution imagery (up to 1 m of ground resolution) which can be acquired by new sensors.One of the main advantages of radargrammetry, compared with the interferometry, is the less influence by atmospheric effects.Basically, atmospheric effects on the SAR imagery are same in the StereoSAR or in the InSAR.However, radargrammetry uses the magnitude (intensity) value while InSAR uses the phase difference in SAR imagery.Considerably, magnitude is less affected than phase by atmospheric heterogeneous.The atmospheric disturbance is undesirable for the InSAR processing but not much of a concern for the radargrammetry processing (Yu et al., 2010).
Formerly, the radargrammetric approach was less and less used, due to the quite low resolution in amplitude of radar image, if compared to their high resolution in phase (Leberl, 1990).Actually, the interest for the radargrammetric approach to Digital Surface Models (DSMs) generation has been growing thanks to the availability of very high resolution imagery acquired by new SAR (Synthetic Aperture Radar) sensors, as COSMO-SkyMed, Radarsat-2 and TerraSAR-X, which are able to supply imagery up to 1 m GSD (Ground Sample Distance) (Toutin and Gray, 2000, Capaldo et al., 2011b, Perko et al., 2011).In particular, the Italian COSMO-SkyMed mission is a four-satellite constellation (the last was launched in 2010), each equipped with an X-band SAR sensor; the Canadian Radarsat-2 platform (launched in 2007) follows on from Radarsat-1 and operates at C-band radar; the German TerraSAR-X (launched in 2007) and TanDEM-X are a twin satellite with a high-resolution X-band SAR sensor payload.Moreover, the incoming ESA (European Space Agency) mission, called Sentinel-1, is equipped with a C-band imaging radar and will become a valid free source of this kind of data.
As concern the radargrammetric processing, SAR imagery are characterized by proper deformations due to the different acquisition geometries and processes, which should be duly taken into account during the two fundamental steps for DSMs generation (image orientation and image matching) if their potentialities have to be fully exploited (Fayard et al., 2007).
Like for optical high resolution satellite data, SAR imagery orientation may be performed using rigorous or Rational Polynomial Functions (RPFs) models with related coefficients (Rational Polynomial Coefficients -RPCs).The former approach tries to International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-1/W1, ISPRS Hannover Workshop 2013, 21 -24 May 2013, Hannover, Germany model the physical image acquisition processes, the latter is a purely analytical approach based on the RPCs supplied by the image providers.
As regards stereo geometry, the optimum configuration for the radargrammetric application is when the target is observed in opposite-side view; however it causes large geometric and radiometric disparities hindering the image matching procedure.A good compromise is to use a same-side configuration stereo pair with a base to height ratio ranging from 0.25 to 2 to increase the efficiency in the correlation image process (Meric et al., 2009).
In particular, it is well known that different image matching approaches have been developed within the photogrammetry and computer vision research fields.In all matching algorithm, there are two fundamental aspect that must be take into account (i) the definition of a primitive model and consequently of an identification criterion (ii) the choice of a strategy for the search of homologous points on a couple of images.Nevertheless, for SAR imagery, apart from the proper deformations already mentioned, also the typical "salt and pepper" aspect due to speckle noise has to be duly considered to carry out a successful matching.
The goal of this paper is to test two SAR orientation models, the rigorous Toutin's model and the SISAR radargrammetric model, and two matching algorithms, implemented in PCI Geomatica and in SISAR, using data acquired over the test site of Beauport (Quebec, Canada) by three different high resolution satellite, COSMO-SkyMed, TerraSAR-X and Radarsat-2.The former orientation model has been developed at the Canada Center for Remote Sensing (CCRS) and implemented in the commercial software package PCI Geomatica, the latter and the SISAR matching have been developed by the research team of Geodesy and Geomatics division of the University of Rome La Sapienza and implemented in the scientific software SISAR (Software Immagini Satellitari ad Alta Risoluzione).
The two radargrammetric procedures follow different approaches both for the orientation model and matching algorithms.The main difference between two orientation models is the use of GCPs: Toutin's model uses at last 8 GCPs to orientates a stereopair while SISAR orientates the SAR images without using a priori ground truth.In addition, the matching algorithm implemented in SISAR is based on a coarse-to-fine hierarchical solution with an effective combination of geometrical constrains and an ABM algorithm, with an original procedure and implementation.The PCI-Geomatica matching method uses quasi-epipolar geometry to reduces possibility of incorrect matches and a crosscorrelation criteria.It is important to underline that, regarding the DSMs extraction only preliminary results are here reported and discussed.
After a first description of the several images and the reference data, main SISAR features and a brief introduction to PCI Geomatica approach are reported.Finally two different kind of comparison and the relative results are shown: on the orientation models performances and on the DSMs generation and assessment, the latter performed using the software DEMANAL developed by Prof. K. Jacobsen at Leibniz University Hannover.

DATA SET
The available data for the experiments on Beauport (Canada) test site were COSMO-SkyMed (CSK), TerraSAR-X (TSX) and RA-DARSAT-2 (R2) imagery (Tab.1).In detail, the following stereo pairs have been considered: a same side CSK in Spotlight Mode, supplied by e-Geos; a TSX in Spotlight Mode supplied by CCRS; a R2 in ultrafine Mode supplied by CCRS, through German Aerospace Center (DLR).The R2 images cover a 20×25 Km 2 area with a slow resolution at level of 3×3 m 2 , the CSK and TSX images cover a smaller area 10×10 Km 2 with different ground resolution, respectivley 1×1 m 2 and 2×2 m 2 (Tab.1).The overlap area of the three stereo pairs is shown in Fig. 1.

SISAR PACKAGE
The radargrammetric processing chain implemented in SISAR is outlined hereafter with particular attention to the orientation model and to the image matching strategy adopted.

SISAR orientation model
The radargrammetric model implemented in SISAR is based on the equation of radar target acquisition and zero Doppler focalization.The radargrammetry technique performs a 3D reconstruction based on the determination of the sensor-object stereo model, in which the position of each point on the object is computed as the intersection of two radar rays coming from different positions and therefore with two different look angles.
Actually, these radar rays can be simply modeled as two segments of measured lengths centered in two different positions (along two different satellite orbits), so that the intersection generating each object point is one of the two possible intersections between two circumferences centered in the two different positions and laying into two planes orthogonal to the two different satellite orbits whose radii are equal to the segment measured lengths.In zero Doppler geometry the target is acquired on a heading that is perpendicular to the flying direction of satellite; the second equation is the slant range constrain (Capaldo et al., 2011a).The couple of equations in a ECEF system reads: where: • X P , Y P , Z P are the coordinates of the generic ground point P in the ECEF coordinate system • X S , Y S , Z S are the coordinates of the satellite in the ECEF coordinate system • u X S , u Y S , u Z S are the Cartesian components of the satellite velocity in the ECEF coordinate system • D s is the so-called "near range" • ∆r is the slant range resolution or column spacing • I is the column position of point P on the image Moreover, the time of acquisition of each lines can be related to line number J through a linear function, since the satellite angular velocity can be considered constant along the short orbital arc related to the image acquisition where start time is the time of start of acquisition, PRF is the Pulse Repetition Frequency, t is the time of line acquisition and J the corresponding line number.
The first step for the image orientation is the orbital estimation; the goal is to estimate the satellite position for each line number according to zero Doppler geometry.In the metadata file, available with SAR imagery, the ECEF position ad velocity of satellite related to the time are supplied at regular interval through state vectors, whose number depends on the considered SAR sensor.
The orbit interpolation has been performed by Lagrange polynomials: the polynomials degree depends on the state vectors number being one unit lower; in particular for COSMO-SkyMed data are available 15 state vectors, for TerraSAR-X 12 state vectors and for RADARSAT-2 5 state vectors.In this way it is possible to estimate the stereo orientation without GCPs.This involves the considerable advantage since it is not necessary to know a priori ground information and to select the points on the SAR imagery, which may be a difficult operation due to the speckle affecting the radar imagery.In fact, as is possible see in the Fig. 3 the homologous points individuation on radar images is much more difficult than in the optical one, so that the point positions may result affected by significant errors.

SISAR matching strategy
It was underlined that the development of a fully automatic, precise and reliable image matching method that adapts to different images and scene contents is a challenging problem.Dissimilarities between images due to occlusion, geometric distortions, radiometric differences and speckle noise must be take in account and this is one of the reasons why many different image matching approaches have been developed in recent years (Gruen et al., 2006).
Generally, a matching algorithm is composed of two essential parts: a primitive model to identify a correspondence between the pixels of the two or more images and a search strategy to find the matching candidates.An innovative image matching algorithm, presently under patenting by the University of Rome "La Sapienza", has been developed.It is based on area based primitive model and on an hierarchical solution with geometrical constrain.The correspondences are looked analysing the signalto-noise ratio (SNR) along two perpendicular search paths.The leading idea, that guided the development of this algorithm, has been to search the homologous primitives directly in the object space re-projecting and re-sampling the stereo images into a 3D ground grid.In this way it is possible not only to limit the research space, but also fetch the images in the same ground geometry, allowing a more easier, robust and reliable homologous points recognition.Moreover, experimental results have highlighted that an image enhancement should be consider in order increase the number of matched points; in this work different The model is based on precise metadata information and reflects the physical reality of the complete acquisition geometry taking in account all the distortions generated during the image formation.Toutin's model is based on simple and straightforward equations where some unknowns parameters should be estimated with few ground control points (at least 8, but usually 12) (Toutin, 2004).
OrthoEngine matching algorithm need epipolar or semiepipolar image geometry to check the homologous points.Three different parameters influence the final radargrammetric DSMs accuracy: Epipolar downsample, details and DTM sampling pixel.
The Epipolar downsample represents the number of image pixels and lines that will be used to calculate one epipolar image pixel; the details are the level of detail that you want in the extracted DEM; the DTM sampling pixel is the number of image pixels and lines (sampling frequency) that will be used to extract one DEM pixel.This is the only informations that is possible to declare since the PCI-Geomatica algorithms are under protection (PCI, 2012).

SISAR orientation model vs Toutin's model
To test the effectiveness of the rigorous models implemented in SISAR and in OrthoEngine, the stereo pairs have been orientated and the model accuracy has been evaluated computing the RMSE over CPs residuals (RMSE CPs).
One of the most important difference between the SISAR radargrammetric model and the Toutin's model is the use of GCPs.
As previously mentioned, Toutin's model needs at least 8 GCPs.Actually, in our these tests a particular set of 12 GCPs has been selected to orientate the three stereo pairs.This set of points has been chosen following a particular criterion.The GCPs should be homogeneously distributed both horizontally and in the height.
On the contrary, the SISAR model, based on metadata information only, is able to orientate the SAR imagery without using any GCPs.
Specifically, for the RADARSAT-2 data a set of 60 Ground Points (GPs) acquired with GNSS technology is available.In Figure 4  The horizontal accuracy of SISAR software is at level of 2-4 m for RADARSAT-2 data, 4-6 m for TerraSAR stereo pair and 2-3 m for COSMO-SkyMed imagery.In particular, the vertical accuracy is better than 3 m for the whole sensors and it is important to underline that the scientific software presents a large average both in up and planimetry.
As regards OrthoEngine software, the orientation accuracy of RADARSAT-2 images is at level of 2 and 3 m in horizontal and vertical directions respectively; for TerraSAR-X products, is around 3-4 m both in horizontal and in vertical; the COSMO-SkyMed images were orientated with an accuracy better than 2 m in up and around 3 m in planimetry.The results (Tab.2) shown that the Toutin's model accuracy is better than the SISAR model one, overall for COSMO-SkyMed and Radarsat2 data and the difference is about 1 m in up.For TerraSAR-X product, SISAR give better result in up but in East the accuracy is poor, at level of 6 m.
This comparison should consider that the intrinsic high accuracy of GNSS and LiDAR points is compromised by the collimation errors, and this phenomena could affect the statistical evaluation.
In fact, it has to be pointed out that the identification of GPs on the SAR imagery is usually much more difficult than in the case of optical imagery, so that an average error of 2-3 pixels (if not larger) should be considered.

DSMs generation and assessment
Afterwards, having studied the different behaviour of the orientation models, the next step has been analysed the DSMs extracted using the two different radargrammetric approaches.Preliminary tests have been carried out in order to characterize the radargrammetric mapping potential of the different sensors using the three available stereo pairs (see table 1) over the Beauport test site.Several DSMs have been generated over the overlap area (see figure 1); the tile is characterized by a flat forested area and a small build up zones near the St. Charles lake.
Subsequently, the DSMs have been assessed through the scientific software DEMANAL developed by the Prof. K. Jacobsen of University Leibniz of Hannover (Jacobsen, 2005).The height differences has been computed re-sampling the generated DSMs over the LiDAR one using a bilinear interpolation method and the accuracy statistics (RMSE, bias, standard deviation) has been evaluated at the 95% of probability level.
As regards the OrthoEngine DSMs extraction, it is important to underline that several tests has computed.This tests were carried out under the supervision of Prof. Thierry Toutin at CCRS in order to select the optimal parameters for radargrammetric DSMs generation.In particular, the imagery has been preprocessed using an Enhanced Lee adaptive filter (embedded in OrthoEngine software) with the aim of reduce the speckle noise and make easier the matching process; COSMO-SkyMed and RADARSAT-2 have been filtered using a template windows of 13x13 pixels; TerraSAR-X with a template window of 11x11 pixels.Moreover, the epipolar imagery generated has been downsampled of a two factor for all the tests.
OrthoEngine software does not allow to extract point clouds; the generated DSMs are sampled over a regular ground grid with a horizontal posting of 4x4 meters.
In table 3 the assessment results are shown; the accuracy level in term of RMSE is about 4 m for all the three platforms; the only differences are the slightly higher values of bias obtained with TerraSAR-X stereo-pair and the high value of blunder, about 250 meters, reported for Cosmo-SkyMed DSM.
As concern the SISAR processing, the OrthoEngine filtered imagery have been used in order to make a comparison using the same preprocessing algorithm.Some preliminary DSMs have been generated only with COSMO-SkyMed and TerraSAR-X data.
Starting from the native point clouds generated by SISAR, the corresponding regular DSMs have been sampled over a regular ground grid (4x4 meters) using a Kriging interpolation algorithm.Both, point clouds and regular DSMs have been assessed through DEMANAL software.
As shown in Table 3 the accuracy level of point clouds is about 3 meters for TerraSAR-X and 4 meters for COSMO-SkyMed stereo-pairs.The RMSE values grow up to 4-4.5 meters in the regular DSMs and they are practically the same with respect to those obtained with OrthoEngine software.

CONCLUSIONS AND FUTURE WORK
The purpose of this work was the comparison of two different radargrammetric approaches implemented, respectively, in the scientific software SISAR and in the OrthoEngine tool of the commercial software PCI Geomatica.Considering the two main steps of the radargrammetric processing, in this paper an orientation model analysis and an assessment of generated DSMs have been illustrated.
At first, in order to compare and evaluate the level accuracy of 3D reconstruction obtained from the two orientation models, several analyses have been carried out on three stereo pairs acquired on Beuport (Canada) test site by three different platforms, i.e.COSMO-SkyMed, TerraSAR-X and RADARSAT-2.The main difference between the two orientation models is that the SISAR model used only the metadata information, whereas the Toutin's model needs at least 8 GCPs to refine some unknown parameters.), respectively in the East, North and Up direction.Particular, Toutin's model were estimated using 12 GCPs selected from the available set of CPs for all three stereo pairs.In order to evaluate the accuracy using the same set of check points, for both model the residuals have been computed over the remaining CPs.
The results show that the RMSE value ranges about from 2 to 4 meters in Up direction for all the tests and the accuracy level appears consistent with the imagery resolution.Particularly, the Toutin's model accuracy is slightly better than the SISAR one; this comparison should consider that the intrinsic high accuracy of GNSS points is compromised by the collimation errors, and this phenomena could affect the statistical evaluation.
Moreover, to define the real effectiveness of radargrammetric techniques the imagery acquired on Beauport test site have been totally processed and several DSMs have been generated.The DSMs accuracy has been evaluated through a reference LiDAR ground truth using the scientific software DEMANAL developed by Prof. K. Jacobsen at Leibniz University Hannover.
A tile characterized by a flat forested area and a small build up zones has been selected considering the overlap area between the imagery.The DSMs global accuracy reached in term of RMSE is about 4 meter for all the used sensors and, considering the smooth terrain morphology, the accuracy level is comparable to the orientation results.As regards the comparison between the two radargrammetric approaches, the results highlighted similar performances, even if higher standard deviation values, instead of lower bias ones, have been observed in SISAR DSMs compared to OrthoEngine DSMs.These slightly discrepancies are probably due to the different interpolation algorithms used in the two processing chain.
This comparison should be widened to larger areas possibly characterized by more complex morphology and unfavourable land cover, like mountainous forested areas.Moreover, additional tests should be performed both on Stripmap images and on data acquired by coming SAR satellite sensor, as Sentinel-1, an ESA C band imaging radar mission to provide imagery for GMES (Global Monitoring for Environment and Security) user services.

Figure 1 :
Figure 1: Beauport overlap area On the Beauport area a set of 60 Control Points (CPs), acquired by GNSS survey, are available and an other set of 20 points was obtained from LiDAR DSM.In fact, the ground truth of the area test is a DSM acquired by LiDAR technology (Fig. 2).

Figure 3 :
Figure 3: Example of situations faced with during point selection on radar imagery

International
Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-1/W1, ISPRS Hannover Workshop 2013, 21 -24 May 2013, Hannover, Germany speckle filter methods (Lee, Kuan, GammaMap) have been investigated and embedded in the radargrammetric software.A specific speckle dynamic filtering technique has been designed and implemented into the radargrammetric processing chain.4 PCI-GEOMATICA PACKAGE PCI-Geomatica v. 2012 is a commercial software that provides support for common GIS and image processing tasks.Using the implemented tool OrthoEngine, it is possible to extract Digital Surface Model by optical and SAR imagery.The Ortho-Engine radargrammetric processing is composed by an orientation model, developed by Toutin, and a image matching algorithm.

Table 3 :
International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-1/W1, ISPRS Hannover Workshop 2013, 21 -24 May 2013, Hannover, Germany Beauport overlap area -CSK, TSX, R2 DSMs accuracy comparison Starting from the available set of GNSS CPs, the models accuracy has been assessed in term of Root Mean Squared Error (RMSE