UPGRADE OF FOSS DATE PLUG-IN: IMPLEMENTATION OF A NEW RADARGRAMMETRIC DSM GENERATION CAPABILITY

Synthetic Aperture Radar (SAR) satellite systems may give important contribution in terms of Digital Surface Models (DSMs) generation considering their complete independence from logistic constraints on the ground and weather conditions. In recent years, the new availability of very high resolution SAR data (up to 20 cm Ground Sample Distance) gave a new impulse to radargrammetry and allowed new applications and developments. Besides, to date, among the software aimed to radargrammetric applications only few show as free and open source. It is in this context that it has been decided to widen DATE (Digital Automatic Terrain Extractor) plug-in capabilities and additionally include the possibility to use SAR imagery for DSM stereo reconstruction (i.e. radargrammetry), besides to the optical workflow already developed. DATE is a Free and Open Source Software (FOSS) developed at the Geodesy and Geomatics Division, University of Rome “La Sapienza”, and conceived as an OSSIM (Open Source Software Image Map) plug-in. It has been developed starting from May 2014 in the framework of 2014 Google Summer of Code, having as early purpose a fully automatic DSMs generation from high resolution optical satellite imagery acquired by the most common sensors. Here, the results achieved through this new capability applied to two stacks (one ascending and one descending) of three TerraSAR-X images each, acquired over Trento (Northern Italy) testfield, are presented. Global accuracies achieved are around 6 metres. These first results are promising and further analysis are expected for a more complete assessment of DATE application to SAR imagery.


INTRODUCTION
Given their complete independence from logistic constraints on the ground and both weather and illumination conditions, Synthetic Aperture Radar (SAR) satellite systems may give an important contribution in terms of Digital Surface Models (DSMs) generation.Starting from SAR data, two different approaches to extract elevations can be considered: the phase-based interferometric technique and the intensity-based radargrammetric one.As deeply explained in (Gelautz et al. 2003), both interferometry and radargrammetry offer advantages and drawbacks.Among others, the main advantage of the interferometry technique rely on the fact that its accuracy depends on the signal wavelenght, which is smaller than the amplitude resolution and can therefore provide higher accuracy; on the other hand, to properly work, radargrammetry does not require coherence between images.At any rate, in order to perform a deeper analysis and obtain the best results both in terms of accuracy and completeness, it would be worthy considering as complementary the result obtained through the radargrammetric approach and the interferometric one, as described in (Crosetto et al. 1999;Gelautz et al. 2003;Sansosti 2004).This work is focused on the description of a new radargrammetric capability for DSMs generation, implemented in the Free and Open Source (FOSS) DATE plug-in, developed at the Geodesy and Geomatics Division of the University of Rome "La Sapienza".The radargrammetric technique use for DSM generation is not a new topic; as a matter of fact radargrammetry was first employed in the 1950s with ground and airborne radars, and then less and less exploited, due to the low amplitude resolution of SAR im- * Corresponding author agery, since the accuracy of radargrammetry-derived products directly depends on the spatial resolution of the SAR image.In recent years, mainly due to the availability of high-resolution SAR imagery collected by new satellite missions such as COSMO-SkyMed, TerraSAR-X , and RADARSAT-2, radargammetry has begun to deserve renewed attention.As a matter of fact, very high resolution SAR data reach up to 20 cm Ground Sample Distance (GSD) and have given new impulse to radargrammetry, having allowed new applications, developments and investigations (Toutin et al. 2009;Perko et al. 2011;Capaldo et al. 2012;Nascetti et al. 2015).
In the following, we illustrate the new radargrammetric capability implemented in DATE software and discuss the results obtained with two stacks (one ascending and one descending) of three TerraSAR-X images each, acquired over Trento (Northern Italy) testfield, which is characterized by very steep slopes and mixed morphology and land cover.
In details, in Section 2 a description of DATE characteristics and functionalities is provided, paying specific attention to the workflow implemented for SAR imagery processing; section 3 includes the data set description and the results presentation and discussion.In the end, in Section 4, some conclusions are drawn and future prospects are outlined.OSSIM is a Free and Open Source Software for Geospatial (FOS-S4G) (Moreno-Sanchez 2012), that appears as a set of libraries and applications for image, maps and terrain data processing; it belongs to the wide OSGeo family, has been actively developed since 1996 and offers a great amount of features and a very robust architecture (OSSIM wiki 2015).Amongst other benefits, OSSIM is able to handle a large amount of data as it manages to process different images in a sequential and totally automatic way.Having been developed as an OSSIM plug-in, also DATE appears as a FOSS4G, that is its source code is accessible, usable for any purpose, modifiable and redistributable.An Open Source Software (OSS) carries various benefits: if everyone can access the source code this leads to a rapid development of high class software, increasing the stability by skilled community review.Furthermore, in general an OSS appears with a high degree of modularity that enables an easy functionalities expansion and has the potential of breaking the project into small pieces and assigning them to different developers (Câmara et al. 2007).

DATE WORKFLOW DESCRIPTION
DATE is a tool for fully automated DSM generation from High Resolution Satellite Imagery (HRSI).The developed tool is based on a hybrid approach, whereby photogrammetric and computer vision algorithms are mixed in order to take the best from both, with the aim to develop an efficient, precise and accurate tool.The synergy between computer vision and photogrammetry has taken hold over a decade ago and it has already returned outstanding results (Pierrot-Deseilligny et al. 2011).
An innovative approach based on a coarse-to-fine pyramidal scheme has been adopted to take advantage of iterative solutions at gradually increasing resolution in order to obtain the "quasi-epipolar" images from the original stereopairs.The procedure is iteratively applied, increasing the resolution when passing from a pyramidal level to the next one, until the desired resolution is achieved.It is quite evident that in this way at lower resolution it is possible to detect larger structures (low frequency) whereas at higher resolutions small details are progressively added to the already obtained coarser DSM.
Overall, the DSM generation processing chain consists of the following main steps hereafter illustrated: • images orthorectification

• dense image matching and disparity map generation
• DSM geocoding Raw satellite images are projected in a ground geometry using a coarse DSM (e.g.SRTM, ASTER GDEM) or other freely available DSMs for the area of interest; this is done using the image metadata Rational Polynomial Coefficients (RPCs) through the OSSIM built-in RPCs-based orientation model (OSSIM wiki 2015), available for the most common sensors and easily upgradable to the upcoming ones due to the modular structure of this tool.The orthorectified images generated in this way can be considered as "ground quasi-epipolar" at each pyramidal level, since the lower resolution adopted allows to approximately satisfy the epipolar constraint, with maximum errors on the transversal parallax of 2 pixels, without making use of the rigorous epipolar geometry, that is not straightforward with imagery acquired by pushbroom sensors differently from standard imagery with one projection center only.
At each pyramidal level a disparity map for the stereopair is computed using the Semi Global Matching (SGM) algorithm by Hirschmüller (Hirschmüller 2010) as implemented in the OpenCV library (Dröppelmann et al. 2010); this is actually the Semi Global Block Matching (SGBM) algorithm, that differs from the original one for considering only 5 research directions instead of 8, for matching not individual pixels but blocks of pixels (even if changing the default parameter it is possible to reduce the blocks to single pixels) and for having implemented a simpler Birchfield-Tomasi sub-pixel metric (Birchfield et al. 1998) instead of a mutual information cost function.In Figure 1 an example of a disparity map computed from a TerraSAR-X stereopair on the Trento testfield is reported.The obtained local pixel disparity values are then converted into height differences with respect to the used approximated DSM through the computation of a conversion factor based on RPCs model.In particular a constant conversion factor is estimated starting from the metadata RPC orientation parameters on a 9 x 9 ground regular grid.Specifically at each point two heights are computed, a minimum height (equal to the lowest height in the image) and a maximum one (equal to the highest height in the image) in order to be sure to analyse the complete height range of the considered area.Through the RPCs the corresponding minimum and maximum image points in master and slave images are computed from these minimum and maximum ground points, then the difference between master and slave (along i direction) between these minimum and maximum image points is calculated.The value thus obtained is related to the height difference for the minimum and maximum ground points.The conversion factor used to obtain metric disparity is the mean of these values achieved at each point.It has been evaluated that the conversion factor computed in this way causes a maximum error in the correction of the initial reference DSM in the region of decimetre, consistent with the desidered final accuracy.Note that no Ground Control Points (GCPs) were used, so that the DSMs georeferencing was based on the metadata RPCs only.
As last step of the deployed workflow, the metric disparity values obtained starting from local pixel disparity through the usage of the conversion factor, are summed up to the height of the input DSM at each pyramidal level iteration.The geocoded DSM is then saved in a standard geotiff format, exploiting OSSIM memory handling funcionalities.

Radargrammetric capability
Due to SAR imagery relevance in many acquisition conditions (see Section 1), a new radargrammetric functionality has been added to DATE capabilities.The main idea is the same as for optical imagery: a preliminary image orthorectification, a subsequent dense image matching and disparity map generation to finally obtain a geocoded DSM.
In order to use the same processing chain developed for the optical images, it is necessary to somehow retrieve RPCs for SAR images, since, in general, most of them are not supplied with RPCs.To do this, SISAR (Software per le Immagini Satellitari ad Alta Risoluzione), a scientific software developed at the Geodesy and Geomatics Division of the University of Rome "La Sapienza" (Capaldo et al. 2011) has been adopted.As a matter of fact, it contains a tool for RPCs generation based on the radargrammetric rigorous orientation model.In principle, RPCs can be generated by a terrain-depended scenario without using the physical sensor model or by a terrain-independent scenario.Nevertheless, the first approach is not recommended for two relevant reasons (Capaldo et al. 2014).SISAR follows the second approach and is able to generate RPCs using a known physical sensor model.The procedure to generate RPCs implemented within SISAR includes three main steps: 1) at first the image is orientated through the already established radargrammetric rigorous orientation model; 2) further, 3D object grid with several layers slicing the entire terrain elevation range is generated; the horizontal coordinates of a point of the 3D object grid are calculated from the image corner coordinates with a regular posting, the corresponding (i, j) of the image grid are calculated using the computed orientation model; 3) finally, the RPCs are estimated in a least squares solution, having as input the coordinates of the 2D image grid points and of the 3D object grid points (Capaldo et al. 2014).The RPCs generated in this way are used with other proper imagery metadata files for SAR imagery processing with DATE software.

DSMs merging
In order to exploit information both from ascending and descending SAR images stack (often available for SAR images), a strategy for DSMs merging has been developed.After the single DSM geocoding phase, each generated digital model of the stack is fused with the others: a variability map is extracted for each pixel, and the final pixel values is obtained with a weighted computation of the contribution of each DSMs, according to Equation 1.In such a way a global, more accurate and complete DSM is generated exploiting all the information available from the stack.
where hi,j = pixel merged height value σasc = height standard deviation for the asc.stack σ desc = height standard deviation for the desc.stack µasc = height median value for the asc.stack µ desc = height median value for the desc.stack

EVALUATION
The above described workflow for automatic DSM generation from SAR imagery, has been tested with TerraSAR-X images acquired over Trento (Northern Italy) area.The DSMs have been generated with a cell grid at 8.0 x 8.0 m resolution.

Dataset
The available data set consists of two stacks (one ascending and one descending) of three TerraSAR-X images each over Trento (Northern Italy) area.The images were acquired on January 2011, all with a right look side (Table 1).(Agugiaro et al. 2012).
In order to assess the results obtained, a 2.5D comparison with respect to a reference DSM was performed.The reference dataset used to assess the generated DSMs in this work is a LiDAR DSM with grid posting 1.0×1.0m and a mean elevation accuracy of 0.15 m, freely available on the website of "Provincia Autonoma di Trento" (http://www.territorio.provincia.tn.it).

Results analysis and DSMs validation
In order to be able to process SAR images with DATE, RPCs have been generated with SISAR (see Section 2.2).A 3D grid of 10×10×12, to discretize the volume to be investigated, was used.
To be sure that no significant error is introduced with the RPCs use with respect to the zero Doppler SAR model, an assessment using the Tie Points supplied by the vendors and available in the metadata files is carried out for each image.A standard deviation in the order of 10 −3 pixel is obtained; this means that, for all the images, the zero Doppler model is well repeated by the generated RPCs.
The height differences between the DSMs generated with DATE and the LiDAR DSM were computed together with their statistics (standard deviation, root mean square error (RMSE), bias and LE95).The indices used for the DSMs assessment are the statistic parameters usually adopted for the assessment of the global DSMs (as SRTM, ASTER, TANDEM-X DSMs) (Rodriguez et al. 2006;Gesch et al. 2014).It is important to recall that the DSMs were generated on the basis of the metadata RPCs only, without any orientation correction based on GCPs.
Overall, 6 DSMs were generated, 3 from the ascending stack and 3 from the descending one, using all the possible images combinations and always adopting as master the image with the lower incidence angle.The results obtained from the ascending stack are reported in It is evident that the results obtained with the descending stack are quite better than the statistical parameters achieved with the ascending one, both in terms of bias and standard deviation.Besides, a high bias for the DSMs generated from the ascending stack is present.This bias is not compliant with the absolute geolocation accuracy of the zero Doppler SAR model (Capaldo et al. 2014) and, in the near future, the reason of its presence has to be investigated thoroughly.
The statistical results achieved through the merging procedure are presented in  Note that the error map is computed with respect to the LiDAR DSM and its value is positive when the reference LiDAR height is above the extracted DSM.It is quite evident that in the merged DSM error map most of the biggest errors are removed achieving a better result than the single individual DSM, even if in some areas the unfavorable contribution of the ascendent stack lead to higher values in the merged error map than those pertinent to the descendent stack error map.In order to avoid this unwanted errors, an automatic mask for layover and foreshortening areas could be generated for each DSM (Vassileva et al. 2015), leading to an optimal merging of all the available information.

CONCLUSIONS
Due to SAR imagery relevance in many acquisition conditions, a new radargrammetric capability for the FOSS DATE software has been developed and here is presented.The main workflow is the same identified for optical imagery: a preliminary image orthorectification, a subsequent dense image matching and disparity map generation to finally obtain a geocoded DSM.Since SAR imagery. in general, are not supplied with RPCs, SISAR, a scientific software developed at the Geodesy and Geomatics Division of the University of Rome "La Sapienza", is used to generate them.Furthermore, in order to exploit information both from ascending and descending SAR images stack, a strategy for DSMs merging has been developed.
This new functionality has been tested on two stacks (one ascending and one descending) of three TerraSAR-X images each, acquired over Trento testfield; for the merged DSM, global accuracies around 6 metres have been achieved.Better results could be obtained masking each generated DSMs in layover and foreshortening areas, in order to eliminate these difficult areas and obtain a more accurate merged DSM.
Finally, it would be worthy considering as complementary the result obtained through the radargrammetric approach and the interferometric one, but unfortunately, at the moment, an interferometric pair was not available to make feasible this comparison, which could be carried out in the future.

Figure 3 :
Figure 3: (a) Error map of the median DSM of the ascending stack, (b) Error map of the median DSM of the descending stack and (c) Error map of the merged DSM Located in the NorthEast part of Italy, Trento testfield lies in a valley at 200 m height, surrounded by the Alps with peaks at

Table 4 .
It can be seen that a standard deviation less than 6 metres is achieved exploiting all the available DSMs.The merged DSM is shown in Figure2.