PIXEL-BASED APPROACH FOR BUILDING HEIGHTS DETERMINATION BY SAR RADARGRAMMETRY

Numerous advances have been made recently in photogrammetry, laser scanning, and remote sensing for the creation of 3D city models. More and more cities are interested in getting 3D city models, be it for urban planning purposes or for supporting public utility companies. In areas often affected by natural disaster, rapid updating of the 3D information may also be useful for helping rescue forces. The high resolutions that can be achieved by the new spaceborne SAR sensor generation enables the analysis of city areas at building level and make those sensors attractive for the extraction of 3D information. Moreover, they present the advantage of weather and sunlight independency, which make them more practicable than optical data, in particular for tasks where rapid response is required. Furthermore, their short revisit time and the possibility of multi-sensor constellation enable providing several acquisitions within a few hours. This opens up the floor for new applications, especially radargrammetric applications, which consider acquisitions taken under different incidence angles. In this paper, we present a new approach for determining building heights, relying only on the radargrammetric analysis of building layover. By taking into account same-side acquisitions, we present the workflow of building height determination. Focus is set on some geometric considerations, pixel-based approach for disparity map calculation, and analysis of the building layover signature for different configurations in order to determine building height. * Corresponding author. This is useful to know for communication with the appropriate person in cases with more than one author.


INTRODUCTION 1.1 Motivation
Due to the high resolution they can achieve, space borne SAR sensors are gaining importance for applications in urban areas.Several SAR approaches for building height determination have been implemented in the last decade, the more successful relying on mono-aspect or multi-aspect InSAR Data or on Persistent Scatterer Interferometry (PSI).For more details, please refer to (Soergel et al. 2010).Moreover, due to their weather and sunlight independency, such sensors are very attractive, especially in cases of natural or technological disaster.Often, the infrastructures are affected, and a rapid change detection analysis is required in order to help emergency troops rescuing possible victims or finding suitable places for transitional shelters.In such cases, 3D information is useful.The interferometric methods mentioned above present the drawback of using several images up to bigger stacks of acquisitions, taken under the same incidence angle.Due to the time needed for acquiring such data, they are not suitable for rapid application.However, operational SAR sensors like TerraSAR-X, TanDEM-X, and COSMO-Skymed, allow a fast revisit time of the same area.This time does not exceed two days considering each satellite alone and can be reduced to a few hours considering multi-sensor constellations.Thus, it is possible to analyse urban areas with several acquisitions taken under different incidence angles within a short time span.In this paper, we exploit this opportunity, showing the possibility of rapid building height determination by radargrammetry, using data from TerraSAR-X.

State-of-the-Art
The wide field of radargrammetry was first investigated by (Leberl et al., 1990).The first methods developments served the 3D mapping of the surface of planet Venus by processing SAR images acquired by the Magellan Mission (NASA).In (Leberl et al., 1994), the authors compare the quality of several intensity-based matching methods used for the Venus processing.Since then, several approaches have been implemented on higher resolution SAR imagery, principally for determining DEM of mountainous area (e.g., Fayard et al., 2007), canopy heights (e.g., Perko et al., 2011), or DEMs of glacier regions (Toutin et al., 2013).These approaches all rely on pixel-based normalised cross-correlation calculation for matching.Moreover, they use images pyramids in order to reduce computing time and obtain robust disparity calculation.Although considering steep slopes, where foreshortening effects occur, these approaches do not consider layover areas.Although well developed for mapping of large rural areas, radargrammetric processing in urban area with high-resolution data is still at an early stage.Existing approaches on this topic can be separated into two groups following different strategies on the matching method.The first consists of using featurebased matching.For example, (Simonetto et al., 2005) classifies bright lines in same-side images and search for crossing points in the detected binary images.A discrete dynamic processing is then used for matching the bright crosses.(Soergel et al., 2009) detects silent lines and points in orthogonal-side images before merging them within a production system.Finally, (Goel et al., 2012), makes use of Bayesian inference for estimating the absolute height of single point scatterers, using at least three same-side images taken under different angles.These approaches show good results for height estimation, but do not consider all the information contained in the building signature ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-3/W3, 2013CMRT13 -City Models, Roads and Traffic 2013, 12 -13 November 2013, Antalya, Turkey This contribution has been peer-reviewed.The double-blind peer-review was conducted on the basis of the full paper.19 of the intensity image.The second group consist of using pixelbased matching, at the example of (Oriot et al., 2003).There, several acquisitions taken under the same incidence angle, but under a circular trajectory (i.e different azimuth angles), are considered.Normalized cross-correlation calculation between overlapping images, followed by geocoding of all resulting height maps, permits to retrieve the DEM.This approach shows good results as well.However, due to the acquisitions configurations, the differences between images being matched are small.Moreover, the different layover lengths and shapes that occur in images taken under varying incidence angles are not considered.
In our approach, we show the applicability of pixel-based matching for building height determination by radargrammetry, considering only the disparity of layover areas.The novelty of our work is given by the exploitation of all information contained in the layover areas.In the following, we mainly focus on same-side stereo-radargrammetric configurations.First, we present the overall workflow of the radargrammetric processing, giving some more details on radiometric and geometric considerations.Second, we present the used matching method (Section 2).Then, we analyse the resulting building signature in the disparity map (Section 3).In Section 4, we show first results of building height determination.As conclusion, we give a discussion about future improvements.

Overview
In this section, we briefly present our overall workflow for building height determination, before explaining in more details two particular steps.This workflow is part of an overall approach for change detection in urban area that we presented in (Dubois et al., 2013).In this approach, pre-event interferommetric data are fused with post-event radargrammetric data in order to determine building changes.In both datasets, relevant building features are extracted (i.e.corner lines and building heights) and compared.For more details about the overall concept, our test area, and the acquired data, please refer to this paper.In the following, we focus on the radargrammetric part of this approach, at building level.Here, we determine building heights and some characteristic features based on the analysis of the radargrammetric disparity map at building location.
Figure 1 shows the workflow for radargrammetric processing.
After calibration (see Section 2.2) and resampling, the slave image is coregistered with the master image, in slant range geometry.Then, the disparities between both images are calculated using pixel-based approach (see Section 2.4).A previous analysis of the acquisition configurations allows reducing the search area for matching (see Section 2.3).Using the preliminary extracted building corner line, the obtained disparity map is then filtered, and characteristic features are extracted (see Section 3).Finally, the building height is determined (see Section 4).

Image calibration
As the images are taken under different incidence angles, calibration is a mandatory step in order to minimize the radiometric differences between the images and make them comparable.In our approach, we used the radar brightness β 0 , which represents the radar reflectivity per unit area in slant range (Fritz et al., 2007): Here, k s is the calibration factor and A the magnitude of the considered pixel.After calibration, a resampling of the slave image occurs, so that both images have the same sampling.

Image Coregistration
As explained in (Toutin et al., 2000) and more specifically for building areas in (Dubois et al., 2013), the quality of the disparity calculation depends highly of the acquisition configuration, i.e. if it is same-or opposite-side stereo, if the incidence angles are steep or shallow, and if the convergence angle (intersection angle between both acquisitions) is large or small.In the following, we focus on same-side configurations, considering several incidence and convergence angle configurations.As explained in (Méric et al., 2009), a transformation of the images in epipolar geometry is mandatory in order to limit the search for matches, and so decrease the computation time.In fact, considerations of acquisition configurations can reduce the search along both range and azimuth direction.Nevertheless, considering layover in urban areas, the search for matches is not so trivial.As several scatterers contribute to the intensity of one single image cell, matching of homologous points is a matter of compromises.As we already explained these effects thoroughly in (Dubois et al., 2013), we simply remind here the principal considerations and conclusions for the self-containment of this paper.Figure 2a & 2b show a schematic representation of the effects occurring in layover areas after coregistration, in ground and slant range, respectively.Façade point A is imaged in A' in image m and in A'' in image s.The distance d between both points is the disparity we want to determine.It can be split into two parts: dr, due to the difference of incidence angles of both images, and da, due to the difference of heading angles ζ of both images.Estimating dr and da allows reducing the search area for matches along the range and the azimuth direction, respectively.However, in the layover, matching of the two façade contributions contained in A' and A'' involves matching of ground and roof contributions too, although they do not represent the same scatterers.Determining dr and da gives thus an idea about the matching error induced to ground and roof points.Figure 2c shows da for different façade scatterer heights, for the radargrammetric configuration of incidence angles 21°/52°.This configuration shows the highest difference of heading angles.We can conclude that for increasing façade This contribution has been peer-reviewed.The double-blind peer-review was conducted on the basis of the full paper.20 point heights, i.e. building heights, da increases.In our test area, the mean building height is less than 40m, so that da stays among 1 pixel, which is negligible.Figure 2d shows a summary table of the predicted dr and da for all the radargrammetric same-side configurations we have, considering a building height of 40m.We can observe that da varies between 1 and 5 pixel, depending on the chosen configuration.Furthermore, dr can be very long, for configurations with a very large convergence angle.With these considerations, we can then reduce the search area for matches, as explained is the next section.Moreover, it is also clear that for a point situated on the roof edge, dr corresponds to the difference of layover lengths.

Disparity Map Calculation
In our approach, we decided to show the applicability of pixelbased approaches for matching and calculating disparity of radargrammetric acquisitions, in layover areas.In Figure 3 a schematic representation of the applied method and used notations are given.The grey and blue parallelograms represent the layover areas.Here, we used the normalised crosscorrelation, whereby the local maximum of correlation between a template and the search window corresponds to the matching location.By using the conclusions of previous section, we define an appropriate search area, in order to reduce matching errors and computation time.With this method, matched points belong to the façade, as the main contribution of the layover is due to the façade backscattering.Depending on the configuration, we define the width w search and length l search of the search window as follow:

l l dr
Values of min dr calc can be found in the last row of Figure 2c.It represents the maximum disparity value that can be determined using l search and l temp.But, in order to retrieve the correct building height, min dr calc should have at least the values of Figure 2c.Namely, the acquisition configuration induces the offset dr that has to be recovered in order to determine the correct height.Thus, the search window length has to be adapted, depending of the template window size.For now, the search window is symmetric, centered on the pixel we want to match.However, knowing the direction of displacement, we privilege disparities leading to this direction, if their correlation value is almost as good as the maximum of correlation of the whole search window (at least 95%).
Figure 4 shows first results of disparity map calculation at building location for two different same-side configuration (ascending in blue and descending in grey).Although very noisy, some pattern can be recognized.In the next section, we analyse more in detail these patterns.

BUILDING SIGNATURE ANALYSIS
Figure 5 shows the expected building signature in the disparity map, by using the described method on simulated data.Here, the difference of heading angle ζ has not been considered, as it is almost negligible for our area.Both corresponding master and slave images are also represented.The images shown here consist on binary images that contain uniform random noise.In future work, this step will be improved by using more sophisticated SAR simulation tools.As a general observation, we recognize two parallelograms of homogeneous values, P 1 (green) around the building corner line, and P 2 (blue) around the end of the layover of the master image.These observations can also be done in Figure 4 on real data.However, on real data, the parallelogram situated around the corner line is less noisy than the other one, and seems larger.The higher intensity values of the corner line probably make the correlation calculation more reliable.Moreover, for the configuration represented in Figure 4g & 4h, the lower part of P 2 is influenced by horizontal bright patterns corresponding to parking lots.A closer analysis of the signature of the simulated data though indicates that both parallelograms have the same size.In fact, for template window length l temp smaller than the layover length of the master image l m , the widths w 1 and w 2 of parallelograms P 1 and P 2 can be expressed as follow: Furthermore, both building corner lines and layover borderlines are situated exactly in the middle of their respective parallelogram.Thus, determining the disparities on the middle line of P 2 (layover border line of master) and converting them into height values enables the retrieval of the building height.Figure 4 underlines this conclusion for both configurations: the larger the template window, the larger the parallelograms, and the better both can be recognized.On Figure 5, we can also recognize another homogeneous area behind the parallelogram P 2 (red).Nevertheless, this area is hardly discernible in real data (Figure 4), thus we do not further use it the following.The rest of the building signature in the disparity map of Figure 5 is characterized by noise, the noisy area between both parallelograms having the width: ( ) For different parameterizations, i.e. for a different ratio of l temp and l m , the analysis of the building signature and thus the retrieval of the building height are not so trivial.We will consider it in future work.In the following, we only consider configurations for which (l temp -1) ≤ l m .
Another characteristic of the building signature are the disparities values in each parallelogram.Both parallelograms show homogeneous values, corresponding to the disparity of the line they represent.Thus, P 1 contains disparity values around zero, as the building corner line is located at the same position in both master and slave images, after coregistration.Whereas P 2 shows disparity values equal to the difference of the layover lengths (l' s -l m ).Determining the corresponding heights, thus the mean height value of parallelogram P 2 , lead to the building height.

HEIGHT DETERMINATION
According to the considerations made in Section 2.3 (see also Figure 2a), we can express the disparity d as follow, using basic trigonometric relations: Considering the conclusions made in (Dubois et al., 2013), we can deduce the formula for the height determination: Figure 6 visualises the last steps of the workflow shown in Figure 1, from the disparity map to the height calculation.Here, only the result of the ascending configuration (Figure 4d) is represented.
As the disparity map is very noisy, a filtering of the disparity map is mandatory.First, the building corner line is extracted from the master image automatically.For this, the line detector introduced by (Tupin et al., 1998) is applied, and a hough transform is performed on the detected image.The corner line is then used as input for the filtering (see Figure 6a) that is adapted to the previous version used for interferometric phases (Dubois et al., 2012).In detail, it consists of dynamic masks along the building orientation whereby a coherence based weighting is applied on the phase values.Here, as the data statistic of disparity maps is quite different, we first make use of a simple mean calculation within each filter window.The result of this filtering is shown in Figure 6d.The shape of parallelogram P 2 is clearly recognizable and shows almost homogeneous values.However, the estimated heights are quite under-estimated.Thus, in a second step the correlation values obtained by evaluating the disparities (see Figure 6e) are taken into account.Only the disparity values showing correlation higher than 0.85 are considered as weights in the filtering, (see Figure 6e).Figure 6f shows the height map resulting from this new weighted filtering.In this, P 2 is still clearly identifiable.So far, the parallelogram is defined manually, but its automatic recognition is planned for the near future.An evaluation of the mean height of P 2 leads to a building height of about 30 m.It is still lower than the reference height of 34 m provided by our ground truth data (IGN BDTopo©), but shows an improvement in comparison to the simple mean filtering.The underestimation results from the values at the upper part of the layover (red in Figure 6f), which are due to neighbouring effects.In future work, we will further investigate these effects, in order to distinguish with exactitude those that are due to the neighbouring objects from those that are due to changes in the building structure at the building boundaries.

DISCUSSION AND CONCLUSION
In this paper, we presented a pixel-based approach for building height estimation by radargrammetry relying only on the analysis of the building layover in the disparity map.A ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-3/W3, 2013CMRT13 -City Models, Roads and Traffic 2013, 12 -13 November 2013, Antalya, Turkey This contribution has been peer-reviewed.The double-blind peer-review was conducted on the basis of the full paper.23 weighting of the disparity values during the filtering enables obtaining first estimations of building heights.Our future work will focus on improving this estimation, by enhancing first the disparity map calculation.For this, a hierarchical strategy as described in (Perko et al., 2011) as well as the use of another matching criterion as described in (Tupin et al., 2002) will be investigated.Furthermore, more tests on different buildings will be performed in order to assess this new approach.Moreover, we presented results for same-side configurations.In the near future, we want to investigate the potential of oppositeside configurations, and the combination of both, because the buildings of the test area show high similarity on opposite-sides.Additionally, such configurations permit smaller convergence angle, and thus better similarity.Figure 7 shows opposite façade of the same building and corresponding intensity images.
For the far future, we will use the extracted features and building heights to detected changes between pre-and postevent data.

Figure 1 :
Figure 1: Workflow of the radargrammetric processing

Figure 2 :Figure 3 :
Figure 2: Influence of the heading angle ζ on the coregistration; a) 3D representation of the facade point A mapped in ground geometry; b) schematic representation in top view in slant-range geometry; c) result for configuration 21°/52°; c) overview of the pixel disparities for h=40m, for same-side configurations.

Figure 5 :Figure 4 :
Figure 5: Building signature in the disparity map of simulated data and corresponding notations and ζ the convergence angle.

Figure 6 :
Figure 6: Result of height determination; a) corner line extraction on master image (asc.47°); b) slave image (asc.36°); c) disparity map calculation with l temp =39 and l search =109; d) result of filtering with mean method and observable parallelogram P 2 ; e) correlation map and binary image corresponding to correlation values higher than 0.85; f) height map after weighted filtering; g) corresponding building façade.