Sapienza
Università di Roma
Facoltà di
Ingegneria
Strumenti e Metodiche:
Metodiche: il futuro
3D modeling and Digital Surface Models generation from radargrammetry
an additional resource for the Italian SAR constellation COSMOCOSMO-SkyMed
Paola Capaldo, Francesca Fratarcangeli, Andrea Nascetti,
Francesca Pieralice, Martina Porfiri and Mattia Crespi
DICEA - Area di Geodesia e Geomatica
Università di Roma "La Sapienza”
via Eudossiana, 18 - 00184 Rome, Italy
[email protected]
Dal Sole e dalle stelle ai satelliti per
osservare, misurare, comunicare, viaggiare
ASSOCIAZIONE
NAZIONALE
UFFICIALI
TECNICI EI
Società Geografica Italiana
Villa Celimontana, Roma, 29 aprile 2014
DICEA -Area di Geodesia e Geomatica
Sapienza
Università di Roma
Facoltà di
Ingegneria
Aim of the work
Development a complete radargrammetric approach for the
3D modeling and the generation of Digital Surface Models (DSMs)
ISPRS TWG VII-2 & EARSeL SIGs 3D+ Radar Remote sensing
Illustrate the actual potentialities of DSMs generation from high
resolution satellite Synthetic Aperture Radar (SAR) imagery with
a radargrammetric stereo-mapping approach
• The model has been implemented in the scientific software SISAR (Software
per Immagini Satellitari ad Alta Risoluzione), by the research group of Geodesy
and Geomatic Division - University of Rome “La Sapienza”
• Tests on COSMO-SkyMed SpotLight and Stripmap imagery
DICEA -Area di Geodesia e Geomatica
1
Sapienza
Università di Roma
Outline
Facoltà di
Ingegneria
Radargrammetric technique
• Introduction
• The orientation model: the geometric aspect
• The image matching: the radiometric aspect
DSM Assessment strategy
Como Case Study
• Dataset : SpotLight imagery
• DSM analysis and assessment
Conclusions and future work
DICEA -Area di Geodesia e Geomatica
2
Sapienza
Università di Roma
Facoltà di
Ingegneria
Radargrammetric technique
Radargrammetry:
technique
similar
to
photogrammetry but based on distances instead of
angles; allows a 3D reconstruction starting from a
stereo model
Optimal geometric configurations:
• B/H ratio ranging from 0.25 to 0.7 is
recommended to have a good stereo
geometry
• opposite-side view cause large geometric
and radiometric disparities, hindering the
image matching process
• a good compromise is to use a Same-side
to enable an easier image matching
DICEA -Area di Geodesia e Geomatica
3
Sapienza
Università di Roma
Facoltà di
Ingegneria
Advantages of the Radargrammetry
The advantage of this approach are:
• work with just a couple of images
• short time to collect the data (half day to quite few days) thanks to the
independence of satellite radar acquisition from weather (clouds), daylight and
logistic constraints (as for airborne data collection)
• use the amplitude (not the phase) of the SAR imagery, it does not require
coherence between images, as for the most known and used interferometric
approach (InSAR)
• new high resolution imagery (up to 1 m GSD), which can be acquired by
COSMO-SkyMed, TerraSAR-X and RADARSAT-2 sensors
In order to demonstrate the radargrammetric mapping potentialities of high
resolution SAR data, two test sites was established in the area of Como and
Merano (Northern Italy), characterized by mixed morphology and land cover
DICEA -Area di Geodesia e Geomatica
4
Sapienza
Università di Roma
Facoltà di
Ingegneria
Radargrammetric DSM generation
Geometric aspect
Radiometric aspect
Images
orientation/model
parameters estimation
Images
matching/homologous
points detection
DSM
The two main steps for DSMs generation from SAR imagery according to
the radargrammetric approach are:
•the stereo pair orientation
•the image matching
DICEA -Area di Geodesia e Geomatica
5
Sapienza
Università di Roma
Facoltà di
Ingegneria
Radargrammetric model for SAR imagery
r r r
vS ⋅ ( P − S ) = 0
zero-Doppler constrain
r r
 P − S = DS + CS ⋅ I slant range constrain
v
P
r
S
position of a ground point P
r
vS
satellite position corresponding to the point P
DS
near range
satellite velocity
CS slant range resolution (column spacing)
I
column coordinate of the point P on the image
DICEA -Area di Geodesia e Geomatica
6
Sapienza
Università di Roma
Facoltà di
Ingegneria
Model: Orbit Interpolation details
The orbital arc related to the image acquisition in Spotlight mode is relatively
short (about 10 Km) and fifteen orbital state vectors are available in the
metadata
State vectors
A Lagrange Polynomial Interpolation is
used in order to retrieve the satellite
position
and
velocity
at
the
corresponding row
Lagrange Interpolation
Link between rows and
satellite position&velocity
Model metadata parameters
•start-time, PRF: linear function that relates the
time of acquisition of each GP (Ground Point)
to its line number J (1)
t ( GP ) = StartTime +
•DS: near range (2)
DS (2)
1
⋅J
PRF
(1)
DICEA -Area di Geodesia e Geomatica
7
Sapienza
Università di Roma
Facoltà di
Ingegneria
Radargrammetric DSM generation
Geometric aspect
Radiometric aspect
Images
orientation/model
parameters estimation
Images
matching/homologous
points detection
DSM
The two main steps for DSMs generation from SAR imagery according to
the radargrammetric approach are:
•the stereo pair orientation
•the image matching
DICEA -Area di Geodesia e Geomatica
8
Sapienza
Università di Roma
SAR Radiometric Information
Facoltà di
Ingegneria
The identification of details and features on the SAR images is usually much more
difficult than in the case of optical imagery
This kind of images are affected by a lots of distortions that have to be considered in
order to develop an image matching algorithm
Spotlight, SAR image
Google Earth, aerial image
DICEA -Area di Geodesia e Geomatica
9
Sapienza
Università di Roma
Facoltà di
Ingegneria
Characteristics distortion of the SAR observation system:
FORESHORTENING
Foreshortening effect is caused by the SAR imaging principle: measuring signal
travel time and not angles as optical systems do.
Points a, b and c are equally spaced when vertically projected on the ground.
However, the distance between a’ and b' is considerably shortened compared to
b'-c', because the top of the mountain is relatively close to the SAR sensor
Example of a SAR image affected by foreshortening
DICEA -Area di Geodesia e Geomatica
10
Sapienza
Università di Roma
Facoltà di
Ingegneria
Characteristics distortion of the SAR observation system:
LAYOVER
Layover occurs in the case of a very steep slope (i.e. high buildings, mountain),
targets in the valley have a larger slant range than related mountain tops.
The ordering of surface elements on the radar image is the reverse of the ordering
on the ground. Generally, these layover zones, facing radar illumination, appear
as bright features on the image.
Example of a SAR image affected by layover
DICEA -Area di Geodesia e Geomatica
11
Sapienza
Università di Roma
Facoltà di
Ingegneria
Characteristics distortion of the SAR observation system:
SHADOWING
A slope away from the radar illumination with an angle that is steeper than the
sensor depression angle provokes radar shadows
Shadow regions appear as dark (zero signal) with any changes due solely to
system noise, sidelobes, and other effects normally of small importance
Example of a SAR image affected by shadowing
DICEA -Area di Geodesia e Geomatica
12
Sapienza
Università di Roma
Facoltà di
Ingegneria
Characteristics distortion of the SAR observation system:
SPECKLE NOISE
• For a single surface type (like an
homogeneous terrain) important
grey level variations may occur
between adjacent resolution cells
• These variations generate a grainy
texture, characteristic of radar
images
• This creates a "salt and pepper"
appearance that is called Speckle
SAR raw image
DICEA -Area di Geodesia e Geomatica
13
Sapienza
Università di Roma
Speckle noise filtering
Facoltà di
Ingegneria
Speckle compromised the image matching process and must be reduced:
• SAR image multi-look processing: Independent measurements of the same
target can be averaged in order to smooth out the speckle
• Filtering techniques: Moving window filters are used. Different algorithms have
been proposed to properly shape the impulse response of the filter within the
window (i.e. Lee, Kuan, Gamma-Map filters)
SAR raw image
SAR speckle Image processed with Lee filter (windows 7x7 pixel)
DICEA -Area di Geodesia e Geomatica
14
Sapienza
Università di Roma
Image matching
Facoltà di
Ingegneria
Automatic detection of homologous points
Homologous points: couple of image points related to the same object on the ground
How can we develop an image matching algorithm?
What parameters should we must consider?
DICEA -Area di Geodesia e Geomatica
15
Sapienza
Università di Roma
Image matching
Facoltà di
Ingegneria
Define a matching primitive
The first step of image matching process is to define the matching entity, that is a primitive
(in the master image) to be compared with a portion of other (slave) images, in order to
identify correspondences among different images.
Area Based Matching (ABM): a small image window, composed of grey values, represents the
matching primitive and the principal methods to assess similarity are cross-correlation and
Least Squares Matching (LSM)
Feature Based Matching (FBM): basic features, that are typically the easily distinguishable
primitives in the input images, like corners, edges, lines, are used as main class of matching
Example of two windows primitives
Left: extracted points with Harris operator.
Right: Extracted edges with Canny operator.
DICEA -Area di Geodesia e Geomatica
16
Sapienza
Università di Roma
Image matching
Facoltà di
Ingegneria
Define a search criteria
The second step in developing a matching algorithm is to choose your search criteria, in
recent years many techniques have been developed, among those most used are:
• Epipolar Geometry: images are projected in an epipolar geometry in order to limit the
research space from two dimension to one dimension; the corresponding points are located
along the same epipolar line
• Region Growing Algorithm: the search for correspondence is carried out starting from a
small number of known homologous points and then expands to the rest of the image
• Semi-Global Matching: (SGM): proposed by Hirschmuller (2005 and 2008), successfully
combines concepts of global and local stereo methods for accurate pixel-wise matching at low
runtime
Region growing processing step
Epipolar geometry
DICEA -Area di Geodesia e Geomatica
17
Sapienza
Università di Roma
SISAR matching strategy
Facoltà di
Ingegneria
An original algorithm, which is presently under patent procedure, has been
developed for SISAR :
• based on a coarse-to-fine hierarchical solution with an effective combination of
geometrical constrains and an Area Based Matching (ABM) algorithm
• homologous points are looked for by cross-correlation and signal to noise ratio
thresholds
DICEA -Area di Geodesia e Geomatica
18
Sapienza
Università di Roma
DSM assessment strategy
Facoltà di
Ingegneria
For all the tests performed has been
used the same validation procedure:
• The regular DSMs or the points
clouds have been compared with the
reference
DSM\DTM
through
DEMANAL software, developed by
Prof. K. Jacobsen - Leibniz University
Hannover
• The accuracy, in terms of Root Mean
Square Error (RMSE) was computed
at the 95% probability level, so that
the LE95 was evaluated
Demanal window screenshot
DICEA -Area di Geodesia e Geomatica
19
Sapienza
Università di Roma
Facoltà di
Ingegneria
Data Set COSMO-SkyMed – Como
Area
Como
Acquisition
data
Coverage
2
(Km )
Mean incidence
angles (degrees)
Orbit
Look
side
24/6/201
10 x 10
27.8
Descending
Right
28/6/2011
10 x 10
55.4
Descending
Right
17/6/2011
10 x 10
50.8
Ascending
Right
7/8/2011
10 x 10
28.9
Ascending
Right
B/H
0.8
0.6
Spotlight imagery - zero-Doppler/slant range geometry
These data were acquired within an
ASI COSMO project:
Title: Exploitation and Validation of
COSMO-SKyMed Interferometric SAR
data for Digital Terrain Modelling and
Surface Deformation Analysis in
Extensive Urban Areas
DICEA -Area di Geodesia e Geomatica
20
Sapienza
Università di Roma
Facoltà di
Ingegneria
Reference DSM
• The DSM reference
was acquired with
LiDAR technology
• Horizontal Spacing:
1.0 x 1.0 m
• Vertical Accuracy:
0.25 m
• made available by the
“Regione Lombardia”
LiDAR Reference DSM
DICEA -Area di Geodesia e Geomatica
21
Sapienza
Università di Roma
Facoltà di
Ingegneria
Como DSM accuracy assessment
TILE 1
TILE 2
• Images have been preprocessed using Lee filter with a window size 7×7 pixels
• A descending and an ascending stereo pairs were available, so two different digital
models have been generated
• Two tiles have been selected for the analysis in the Como urban area and a basic
classification of the soil coverage has been performed, wooded (green area) and urban
(blue area), in order to investigate the different behavior of the algorithm
DICEA -Area di Geodesia e Geomatica
22
Sapienza
Università di Roma
Facoltà di
Ingegneria
Como DSM accuracy assessment: Points Clouds
Como CSM Absolute Error [m] on matched points - Tile 1
DSM
BIAS
ST.DEV.
RMSE LE95
LE95
# Points
Total
Ascending
-2.02
7.56
7.83
19.94
224970
Descending
-0.34
7.85
23.57
127665
Ascending
-1.82
7.84
Wooded
6.72
6.96
16.79
97120
Descending
-0.70
Ascending
Descending
6.47
17.31
61520
-2.17
6.43
Urban
8.25
8.53
21.13
127782
0.08
9.35
9.35
25.93
65972
Como CSM Absolute Error [m] on matched points - Tile 2
DSM
BIAS
ST.DEV.
RMSE LE95
LE95
# Points
Total
Ascending
-1.92
8.01
8.24
21.32
242741
Descending
-0.60
10.54
25.24
87287
Ascending
-4.87
10.53
Wooded
7.54
8.97
20.22
25638
Descending
-0.48
10.88
23.99
8340
Ascending
-1.60
10.87
Urban
7.99
8.15
21.32
217192
Descending
-0.66
10.52
10.54
25.05
78374
The two stereopairs have been processed
separately and the relative points clouds
have been assessed, the results remarks:
• in Tile 1 the RMSE of ascending and
descending points clouds is around 8 m
• a better accuracy has been reached in
the wooded area, about 7 m, whereas in
the more complex morphologies of the
urban area the RMSE grew up to 9 m
• in tile 2 the different behavior in
wooded and urban area is not highlighted
and a similar level of accuracy is
detected
• in the descending pair a lower number
of the matched points is detected
DICEA -Area di Geodesia e Geomatica
23
Sapienza
Università di Roma
Facoltà di
Ingegneria
Como DSM accuracy assessment: Tile 1
Starting from the points clouds, the DSMs have been generated, estimating the heights on a 5 m x 5 m grid by a
linear interpolation, after a Delaunay triangulation
A third product has been created by the merging the two opposite side points clouds that have been previously
filtered through the removing of matched points with lower correlation
Merging
DSM Descending
DSM Ascending
DICEA -Area di Geodesia e Geomatica
Merged DSM
24
Sapienza
Università di Roma
Facoltà di
Ingegneria
Como DSM accuracy assessment: Tile 01
Como DSM Absolute Error [m] - Tile 1
DSM
BIAS
ST.DEV.
The results remarks:
RMSE LE95
LE95
# Points
Total
Ascending
-1.07
7.79
7.86
21.87
80849
Descending
1.53
10.24
10.35
33.14
80849
Merged
-1.10
6.94
7.02
18.14
80849
Wooded
Ascending
-0.69
7.10
7.14
18.10
36835
Descending
1.32
8.53
8.63
27.06
36835
Merged
-0.88
6.07
6.14
15.55
36835
• the accuracy on the whole area is
around 8 m and 10 m for the ascending
and the descending DSMs respectively,
and around 7 m in the merged product
• the wooded area is more accurate with
respect to the urban area, and also in
this case the merging produces an
improved accuracy (about 1 m)
Urban
Ascending
-1.45
8.32
8.45
22.01
44107
Descending
1.73
11.86
11.98
38.10
44107
Merged
-1.40
7.59
7.72
20.18
44107
DICEA -Area di Geodesia e Geomatica
25
Sapienza
Università di Roma
Facoltà di
Ingegneria
Como DSM accuracy assessment: Tile 01
A deeper analysis has been carried out in order to understand the difference accuracy
value between ascending and descending stereo pairs :
• the accuracy of ascending DSMs is better than the descending ones, since more
homologous points have been detected on the ascending stereo images
• this is probably due to lower quality of radiometric information in one of the
descending images that present some small zones affect by artifact distortions
Radiometric artefact distortions of descending image
DICEA -Area di Geodesia e Geomatica
26
Sapienza
Università di Roma
Facoltà di
Ingegneria
Como DSM accuracy assessment: Tile 2
Starting from the points clouds, the DSMs have been generated, estimating the heights on a 5 m x 5 m grid by a
linear interpolation, after a Delaunay triangulation
A third product has been created by the merging the two opposite side points clouds that have been previously
filtered through the removing of matched points with lower correlation
Merging
DSM Ascending
DSM Descending
Merged DSM
DICEA -Area di Geodesia e Geomatica
27
Sapienza
Università di Roma
Facoltà di
Ingegneria
Como DSM accuracy assessment: Tile 2
Como DSM Absolute Error [m] - Tile 2
DSM
BIAS
ST.DEV.
The results remarks:
RMSE LE95
LE95
# Points
Total
Ascending
-0.84
8.68
8.72
24.12
74620
Descending
1.06
11.59
11.64
34.62
74620
Merged
-0.99
8.28
8.34
21.49
74620
Wooded
Ascending
-4.80
8.04
9.37
26.59
7972
Descending
1.92
8.25
8.47
19.06
7972
Merged
-4.33
7.61
8.75
20.42
7972
• the accuracy on the whole area is
around 9 m and 12 m for the ascending
and the descending DSMs respectively
• the accuracy is around 8 m in the
merged product
• the wooded area is small and it is not
more accurate with respect to the urban
area
Urban
Ascending
-0.42
8.60
8.61
24.16
66196
Descending
0.97
12.11
12.15
36.73
66196
Merged
-0.61
8.30
8.32
21.39
66196
DICEA -Area di Geodesia e Geomatica
28
Sapienza
Università di Roma
Conclusions
Facoltà di
Ingegneria
• The image matching and the DSM generation have been performed in SISAR
software and the accuracy is strictly related to the terrain morphology: over a
urban area the RMSE is around 7-8 m, the vegetated areas are well recognized
• The use of two same-side stereo pairs acquired from different look side seem to
be a solution to reconstruct the 3D geometry in the presence of foreshortening,
shadows and layover
• Accuracy results are comparable with other research groups results (Toutin,
Raggam et al.), they are at level of COSMO-SkyMed hand-list and interferometric
DSMs, even better (recent call by ASI)
• Radargrammetry is likely to became an effective complement/alternative to
InSAR, since it may work even with a couple of images with good performances
over forested areas too
• This results show that radargrammetry and the InSAR techniques should be
integrated in order to exploit at best SAR data, in particular this method could be
a resource to fill the gaps due to the lack of coherence in interferometric DSMs
DICEA -Area di Geodesia e Geomatica
29
Sapienza
Università di Roma
Facoltà di
Ingegneria
Future work (already ongoing)
• Tests to improve the potentialities of the automatic matching
procedure for DSMs generation in urban areas or in with more
complex morphologies
• More research on filtering techniques (i.e SAR speckle wavelet
filtering)
• Make
direct
comparisons
radargrammetrc techniques
between
interferometric
and
DICEA -Area di Geodesia e Geomatica
30
Sapienza
Università di Roma
Facoltà di
Ingegneria
Thank you for your kind attention
DICEA -Area di Geodesia e Geomatica
31