Calculation of Interferometric Heights

surfaces and building shadow show low values in all three feature layers. Finally, for each InSAR image pair the classification into vegetated and urban areas is accomplished by reasonable adjustment of the fuzzification step.

As high vegetation shows the same layover effect as buildings, adjacent objects can be occluded by vegetation signatures. Hence, a fusion of multi-aspect classification results is carried out to work against this occlusion effect. First of all, this requires an individual geocoding of the classification results to a common coordinate system followed by a fusion step calculating the intersection be-tween the classified vegetated areas. Only those areas are considered as final vegetated areas that belong to this class in all InSAR data. First results on this multi-aspect classification were presented in [123]. The final results as well as the assessment of the classification results by comparison with electro-optical data are presented in Subsection 5.2.3. In the following steps of building reconstruc-tion, the classification result is used to support feature extraction (Section 4.6) and the generation of building hypotheses (Section 4.8).

The next processing step contains the subtraction of ϕ_flat to obtain the so-called flat earth corrected phase. This correction can easily be accomplished using commercial software (e.g. ENVI or ER-DAS IMAGINE) as long as the appropriate sensor model is implemented. The quality of the pro-vided configuration parameters (e.g. rapid or science orbit) determines whether or not post-processing is necessary. If post-post-processing is required, for example, ground control points (GCPs) are utilised. Since InSAR data of the new generation of satellites are provided with high precision orbit information, no additional processing has to be considered (see Fig. 4.6a, step “flat earth cor-rection”). This is different for airborne data, because these platforms show higher variations in track geometry than satellites on orbits. Furthermore, commercial software tools do not support the processing of airborne data in general because of their limited data accessibility and individual specifications. Considering that, for the interferometric processing of airborne data, in-house solu-tions are used.

The flat earth correction of airborne data underlies some assumptions: first, the length of the base-line is constant as the antennas are rigidly fixed. Second, the flight height over ground is steady and third, the flight geometry is stable during recording time. For a scene dominated by flat topography close to the mean height chosen for flat earth correction, a phase distribution with a significant narrow maximum is expected. The result of flat earth correction achieved with the original parame-ter set is given in Fig. 4.6b. The phase distribution (left) shows a wide curve rather than a signifi-cant narrow maximum. The reason for that is a linear phase trend observable in a corresponding range profile (right). This can be caused by imprecise configuration parameters or a natural slope of topography. Taking into account additional information about the topography of the recorded scene, a flat terrain can be assumed. Therefore, the goal of an improved flat earth correction is the reduction of this linear phase trend. This improvement of the flat earth correction contains a re-finement of the two angles ξ and θ, whereby the latter is specified by sensor height over ground

H and range distance r. In detail, ξ is corrected by varying the vertical B_V and horizontal B_H parts of the baseline, and θ by testing H, respectively. The estimation of optimal values is realised by a coarse-to-fine parameter search similar to the coregistration. The aim is to reach a narrow phase distribution with a minimum standard deviation. The reduction of the linear trend is visible in Fig. 4.6c. The numerical results are summarised in Subsection 5.2.2.

The resulting absolute phase values (c) are random and the Gaussian-like phase distribution can show a high σ. The latter can be caused by hilly topography or by insufficient motion compensa-tion during the SAR processing. To find out how strong this effect is, range histograms are calcu-lated from the InSAR phases considering phase values of several neighboured range lines. A linear drift of the distribution maximum becomes visible in the plot of all histograms (Fig. 4.6d left), in-dicating residual motion effects in the InSAR data. Since the fusion of multi-aspect InSAR data is one of the main goals of this approach, this aspect specific effect has to be taken into account. The subsequent phase centring can mitigate this motion effect in an easy way, and additionally a phase distribution with zero mean can be achieved. The centring is implemented as a simple rotation of the complex interferogram values into zero position. Thus, from the interferogram ^{S x y}

(

^,

)

^{with x}

and y as range and azimuth coordinates, the one-dimensional function S yy

( )

is calculated:

( ) ( )

all y ,

x

S y =

∑

S x y ^(4-6)

This azimuth profile is filtered by averaging values in the range of y'= −δ,...,+δ. In the following, this result is termed as S yy

( )

. The parameter δ has to be chosen large enough to avoid effects on the phase of objects. The final zero-centred interferogram results from:

( ) ( ) ( )

zero , , ^y

( )

y

S x y S x y S y S y

∗

= ⋅ (4-7)

where S yy

( )

^∗ denotes the complex conjugate of S yy

( )

. In Fig. 4.6e the centred phase distribution (left) is narrower than the one in (c – left). Also the linear trend in the plot of the range histogram is not longer visible (Fig. 4.6d – right). Of course, care has to be taken for hilly topography as this centring can eliminate topography portions.

Before converting the interferometric phases to height values, phase unwrapping has to take place.

Since many unwrapping algorithms fail in urban areas, a simple step of phase shifting is integrated.

That means phase values significantly below the phase average are shifted upwards by 2π. As the

slant range [pixel]

0 500 1000 1500 2000 -π

+π

phase

-π +π

phase -π0 500 1000 slant range [pixel] 1500 2000 +π

phase

-π +π

phase -π0 500 1000 slant range [pixel] 1500 2000 +π

phase

slant range [pixel]

0 500 1000 1500 2000 -π

+π

phase

-π +π

phase

-π +π

phase

b

-π +π

phase

-π +π

phase

slant range [pixel]

0 500 1000 1500 2000

-π +π

phase

-π +π

phase

slant range [pixel]

0 500 1000 1500 2000

-π +π

phase

-π +π

phase

-π +π

phase

slant range [pixel]

0 500 1000 1500 2000 -π

+π

phase

slant range [pixel]

0 500 1000 1500 2000

c

azimuth azimuth

-π +π

-π

azimuth +πazimuth

-π +π

-π +π

d

-π phase +π

-π phase +π

-π phase +π

e

coregistration SLC data (SSSS₂) SLC data (SSSS1)

interferometric phases

calculation of heights

interferogram calculation coherence calculation

normalised interferometric heights post -processing of phases

phase centring flat earth

correction

phase shifting flat earth correction with iterative adaptation

coherence image interferometric heights

mask generation DTM calculation normalisation of heights post-processing of heights

coregistration SLC data (SSSS₂) SLC data (SSSSSLC data (SSSS1) ₂) SLC data (SSSS1)

interferometric phases

calculation of heights

interferogram calculation coherence calculation interferogram calculation coherence calculation

normalised interferometric heights post -processing of phases

phase centring flat earth

correction

phase shifting flat earth correction with iterative adaptation post -processing of phases

phase centring flat earth

correction

phase shifting flat earth correction with iterative adaptation

phase centring flat earth

correction

phase shifting flat earth correction with iterative adaptation

coherence image interferometric heights coherence image interferometric heights

mask generation DTM calculation normalisation of heights post-processing of heights

a

76.5 height [m] 104.5

76.5 height [m] 104.5

f

Figure 4.6: Workflow of interferometric height calculation (a); histograms of InSAR data and correspond-ing slant range profile – result of “flat earth correction” (b), result of “flat earth correction with iterative adaptation” (c), histogram plot of azimuth drift (d), histogram of “phase cen-tring” and image patch (e), result of “calculation of heights” with “phase shifting” (f)

mean of the phase distribution corresponds to the local terrain height, this step leads to an enlargement of the unambiguous range over ground. Hence, phase jumps (i.e., the transition be-tween +π and −π) at elevated objects, for example buildings, are reduced at the cost of possible exaggeration of sinks in the terrain. Results of the InSAR height estimation without and with con-sidering the phase shifting step are shown in Fig. 4.6e (right) and f, respectively. In the given ex-ample, the 2π unambiguous range corresponds to approximately 30 m and therefore phase ambiguity issues arise at buildings taller than 15 m. With an offset value of ³4π, the unambiguous range of the height is shifted from [66.1 m , 94.0 m] to [76.5 m , 104.5 m]. The phase shifting effect is visible at the U-shaped building (Fig. 4.6e,f right) that is featured by roof heights of 16 m up to 23 m. Obviously, less discontinuities (i.e., transitions from black to white) are observable in the roof area of this building. If no abrupt steep descents are present in the scene, this step of phase shifting compensates phase jumps at building roofs in a proper way.

Considering larger areas with slight topography, the global phase shifting can fail due to varying local terrain heights. For that reason, normalised InSAR heights are generated, which are essential for the subsequent filtering of building primitives (Subsection 4.6.2). The processing of these heights contains the following steps, whose intermediate results are shown in Fig. 4.7.

First, a binary filter mask (mask_i) is computed to define image pixels used for the Digital Terrain Model (DTM) generation. Only pixels of coherence value ( )γi above a minimum coherence threshold (γ_min)and an InSAR height value ( )h_i close to the global mean terrain height ( )h are considered, as given in (4-8). Thus, this mask contains the value 1 for objects like fallow land and rough man-made areas showing high coherence and mean terrain height. Elevated objects such as trees and buildings as well as areas of low coherence (e.g. shadow areas) are assigned to the value 0.

The definition of the parameter γ_min depends on the temporal and spatial length of the InSAR con-figuration. For the given example γmin=0 5. was chosen. Additionally, the parameter σ_h can be es-timated directly from the InSAR heights. In case of rising topography, σ_h will also increase possibly affecting low building signatures. A solution would be the determination of local σ_h.

a b c d e

Figure 4.7: Post-processing of an extended region: InSAR heights (a), coherence (b), mask (c), DTM (d), and normalised InSAR heights (e)

 ≥ − ≤

=



1, if and

0, else

min

i i h

i

h h

mask γ γ σ

(4-8) Second, based on the filtering of the InSAR heights (a) with the binary mask (Fig. 4.7c), the DTM (d) is calculated. The pixel values of this DTM result from a weighted averaging of height values over an area of 50m 50m× in ground range resolution. Such a large area is mandatory for smooth-ing local structures probably included due to a high σ_h.

Finally, the height differences, the normalised InSAR heights (e) using the InSAR heights (a) and the DTM (d), are calculated. In the following, these heights are investigated to filter building primi-tives (Subsection 4.6.2) and to extract building heights (Subsection 4.8.2).

Im Dokument 3D Building Reconstruction From High Resolution Multi-Aspect Interferometric Synthetic Aperture RADAR Data (Seite 87-91)