Proceedings of International Geoscience and Remote Sensing Symposium (IGARSS'96). Lincoln, Nebraskq USA. May 27-31.1996
Algorithms for Calculation of Digital Surface Models from the Unwrapped Interferometric Phase Wolfgang Goblirsch, Paolo Pasquali Remote Sensing Laboratories, Department of Geography, University of Zurich Winterthurerstrasse 190, CH-8057 Zurich - Switzerland fax: +41 1 362 52 27, email:
[email protected]
ABSTRACT where ,441denotes the calibrated unwrapped phase difference and h denotes the wavelength. The range-Doppler equations are:
The focus of this paper is on interferometric DEM generation algorithms from imaging geometry analysis. The proper description of the interferometric imaging geometry is presented. As a result of solving the two sets of range-Doppler equations one also obtains the interferometric alongtrack baseline component. This fixes the imaging positions of each pixel in the case of known antenna tracks. Conversely, if the relative positions of the antenna tracks are not precisely known, one may adjust the interferometric baseline by comparing the azimuth-offsets from image correlation to the calculated alongtrack baseline. As an example, the angle of convergence is estimated using single-pass interferometric DOSAR data motion compensated to individual reference tracks. A new method of estimating the absolute phase difference is suggested which exploits the range-offsets obtained from correlating the complex image pair. Finally, it is shown that layover produces a residual tilt in interferometric DEM's if the crosstrack-coordinate is left unfiltered.
(4)
where:
1. INTRODUCTION Airborne interferometric S A R allows the generation of potentially high resolution Digital Surface Models (DSM) since the interferometric phase is a function of the mapped terrain surface. In order to map the interferometric phase to the surface of the terrain model, two different operations have to be executed.
= (x , y , z) , position of the imaged groundpoint
$ (t,)
= velocity ofantennap
(t,)
= velocity ofantenna q
d
= velocity of the groundpoint
a
= processed squint angle trackp
P
= processed squint angle track q
Geometrically, the two range equations represent two spheres centered at the imaging positions on the two tracks. The Doppler equations represent cones with their tips located at the imaging positions, and their axes pointing in the direction of the velocity vectors. For zero squint the Doppler cones degenerate into zeroDoppler planes.
The first step is the phase unwrapping operation which aims to guarantee the correct multiplicity to the phase difference between 0 and 2n. Due to the inherently noisy interferometric phase this operation is already subject to intense scrutiny by the remote sensing scientific community in order to develop better noise compensating algorithms. The second step transforms the unwrapped phase difference to the surface of the terrain model. Methods are presented for improving this second step.
3. A NEW PHASE CALIBRATION APPROACH Now the final step in interferometric DSM generation is the calculation of the groundpoint coordinates for each matched image pixel on the basis of the antenna tracks and the unwrapped phase difference. In order to do this we first have to find the correct phase difference offset A& to be added to the original unwrapped phase difference Aq0 so that equation (1) will result in the true range difference. This operation is termed phase calibration and if no special calibrakion setup is implemented in the interferometric S A R system usuall!y requires a ground control point (gcp). In certain cases, however, the use of ground control might be avoided using the rangeoffsets Arc obtained from image correlation.
This transformation follows both from two sets of range-Doppler equations representing the imaging geometry, one for each antenna, and from the phase difference relation involving the slantrange difference of the antennas. We will term this set of five equations the Interferometric Range-Doppler Equations (IRDE). In this paper we first focus on the description of the IRDE. Then we suggest a new phase calibration method. From the exact solution of the IRDE follows a baseline adjustment approach. Finally, we demonstrate how separate filtering of the vertical- and crosstrackDSM-coordinates improve the precision of the DSM.
In the case when the range-offsets Arc from the correlation of the compl1:x image pair represent the absolute range difference Ar in equation (1) we are able to calculate Aqk with the help of equation (6):
4zArc A$k = --
2. THE INTERFEROMETRIC RANGE-DOPPLER EQUATIONS
h
A$O
(6)
where Aqk is averaged over all matched pixels with a sufficiently high correlation coefficient. We use the same window-size for filtering A+o used in correlating the complex image pair.
The IRDE describe the imaging geometry of a processed interferometric image pair based on the two antenna flight tracks denoted by p (t) = ( p , ( t ) , p2 ( t ) , pg (t) ) and 9 ( t ) where t is the along-track, or azimuth time, referenced to one track (master track). This time is found with the pixel number j, in azimuth as well as the azimuth pixel spacing 8, as follows: t = (j - 1) 6,,/v , where vp P P P is the averaged velocity on track p. The imaging time tq of the corresponding pixel on the slave track is found with an analogous relation. In a similar fashion the range rp from the master track is found from the slantrange pixel spacing &, the range pixel number i,, and the range rP0 to the first image pixel with: r = rp = rpo + ip8rp.
14s an example, given a range pixel size of 0.624 m and an offset precision of 11100 of a pixel we would get an acceptable error in the phase-offset of about 8' by averaging over 100 pixels. Unfortunately, in airborne interferometric SAR-systems both antennm usually have different internal delay times so there is an
unkown constant range-offset which has to be found by careful rangecalibration. With data from two identical antennas in one-pass interferometry or from one antenna in two-pass interferomertry, the new phase calibration approach should work since the offsets from complex image correlation should offer sufficient precision and the number of pixels used can be quite large.
The range rq to the corresponding slave pixel is obtained from the interferometric phase difference relation:
0-7803-3068-4/96$5.00@1996IEEE
h
656
4. SOLUTION O F THE INTERFEROMETRIC RANGEDOPPLER EQUATIONS
i
Once the absolute range-offset has been established with the help of the calibrated phase difference, we are able to solve the 4 rangeDoppler equations for the groundpoint coordinates d and the nominal azimuth-offset tq - 5. The nominal azimuth-offset disappears for parallel tracks and zero-Doppler imaging geometry in which case we only have to solve the two range equations and one zero-Doppler plane equation. In all other cases we have to solve 4 different equations, two of which are quadratic and the other two linear for the groundpoint coordinates. In the general case of nonlinear tracks, an analytic expression of the solution cannot be given. However, we can always solve two range equations and one Doppler equation in an analytic fashion, starting with an initial azimuth-offset, and then testing the solution with the remaining Doppler equation. This establishes an iterative approach for solving the IRDE with arbitrary precision [3].
0.0
-0.5
o o 200
e00
600
PO0
ni,mutn
Psxel
[O 8
~ 2 0 0
,000
ml
Figure 1: Roll angle variation
For the case of linear tracks a noniterative solution can be given which results from the real roots of a fourth order polynomial in the vertical DSM-coordinate z or the crosstrack DSM-coordinate y. The DSM x-coordinate is conveniently selected in the groundtrackdirection of the master track and is thus given by: x = (pl+rsina)cosO
(7)
g = groundrange to the groundpoint
where 8, denotes the pitch angle of the master track velocity vector. The variable p, = j,S, denotes the alongtrack coordinate of the scene, where jp denotes the alongtrack pixel number and 6 , is the alongtrack pixel spacing. The transformation between the vertical- and the crosstrackcoordinate is found with:
Ag = groundrange difference p, = azimuth-coordinate
Po = projection of slave track squint angle p a, = tangent plane projection of a
I \\
Typically in airborne S A R ,one motion compensates with respect to a fixed sensor altitude for which the azimuth-offset or alongtrack interferometric baseline component B, can be expressed in simpler terms: = e(tany)'-ytany+rsinaX
( r + A r ) s i n p 1 + (tany)
\
1 + (tany)
groundpoint
where:
Figure 2: Imaging geometry in the tangent plane for horizontal linear tracks
y = angle of convergence do = groundtrack separation at the first image pixel
5 . ESTIMATING THE ANGLE OF CONVERGENCE FROM THE CORRELATION AZIMUTH-OFFSETS
e = - (ddtany + pl) = distance to groundtrack intersection By = (B, - e)tany = crosstrack baseline
We exploit equation (9) for an alternative method in estimating y. For airbome two-antennae interferometric SAR systems we usually have 11.1 i0.01" and a = p. In this case equation (9) simplifies to:
The projection of the imaging geometry to the tangent plane is pictured in Figure 2 for right pointing antennae.
B x G - (y tany + Arsina) (11) which shows that the absolute value of the azimuth-offset increases with range in the crossed-tracks case and is also terrain dependent.
In the case of horizontal tracks we require y and do for the calculation of the vertical, or the crosstrack, DSM-coordinate. The variable do can be given as a function of the absolute platform attitude angles at the start of the scene for an airborne single-pass interferometric S A R system. One has to use a gcp for adjusting do if these angles are not precise. If the platform roll angle variation is available we might find the angle of convergence with the following equation:
The azimuth-offset Ax from correlating the complex image pair usually is the along-track baseline component B, plus a constant A,. Thus in order to calculate y we must use the difference of the correlation azimuth-offset obtained at different ranges where y (r2) # y ( r l ) : ( A x ( r l ) -Ax('')) tany =
where B is the physical baseline length, L, denotes the alongtrack length of the processed scene, and Ap is the difference of the platform roll angles at the beginning and at the end of the processed scene segment. This works especially well when the roll variation is nearly linear as demonstrated with IMU data pertaining to the interferometric data set gathered by DOSAR during the 1994 campaign in Switzerland [l], [2]. Figure 1 shows the roll angle variation for a particular scene segment which was used to generate a DSM with the help of the IRDE for crossed tracks. The crossed-tracks geometry was caused by motion compensating with respect to individual reference track-directions.
+ s i n a ( A r ( r l ) -Ar(r2))
Y ('2) - Y ( r l )
(12)
Of course y is still a function of y so we would have to use an iterative approach to find the angle of convergence. But if we are only interested in the order of magnitude of y we can use an approximated crosstrack-coordinate obtained, for example, from a low resolution reference DEM. We correlated the complex image pair processed from 200 MHz DOSAR data [I] with a 16x16 window. The pixel-size of the used data set was 0.58~0.62m2 (azimuthxrange). The offset precision was set to pixel-sizd20. The offsets were updated at 9.3 m intervals in azimuth, and at 10 m intervals in range. When we inserted the corresponding offsets into equation (12) for all '2 rl = 10 m no trend could be detected. Averaging y over all samples resulted in y = O.oOIo. We have yet to confirm this method with additional data sets.
In the given example we obtained 0.0008° for y by inserting the values of B = 1.06 m, Ap = 0.75' and L, = 960 m into equation (IO). The value of y was confirmed by using gcp's for adjustment of the imaging geometry during DSM generation.
-
657
6. LAYOVER AND NOISE INFLUENCE IN THE DSM CROSSTRACK-COORDINATE It is not adviseable to filter the phase difference heavily before DSM calculation because this will reduce singular DSM-elevation values unnecessarily. On the other hand, leaving the phase difference essentially unfiltered will cause DSM deformation in the vertical and crosstrack DSM-coordinates. If we want to measure the elevations of individual objects more precisely we might want to accept a certain degree of noise in the vertical DSM-coordinate z. However, the same degree of noise is unacceptable in the crosstrack-coordinate y since it will turn the slantrange into a less monotonic function of y [l]. But even if we have strongly reduced the noise influence by filtering the interferometric phase will we measure some crosstrack displacements of steep surface features. This can be attributed to the layover effect resulting from the SAR slantrange geometry:
sulted in more or less residual tilt which a crosstrack filtering window size of of 3x3 demonstrates in Figure 5. The filtering window size has to be adjusted to the expected overlaid slantrange intervals. On the other hand, one should not filter too strongly in order to preserve crosstrack topograiphic information as much as possible. As a rule of thumb the maximum window size n for crosstrack smoothing may be calculated with the following equation: n = [ (r2-rl)/6r] (13) where: rl = ,/(yo r2 =
Az and Ay are the vertical and illuminated crosstrack dimensions of the overlaid object, yo is the minimum groundrange, and zrefis medium elevation of the scene.
When the slope facing the radar is greater than the off-nadir angle we have more than one ground location corresponding to the same slantrange; a situation whch is equivalent to a nonmonotonic slantrange as a function of groundrange. As a result, signatures from locations whose ranges correspond to the same range sample mix and we obtain only one elevation value per slantrange pixel, thereby changing the form of the imaged features upon slantrange-tocrosstrack transformation with equation (8).
560 540
Et 520
L A
c:
.E!Cl 500 tl
2; 480 460
Their correlation coefficient reaches values in excess of 0.87 where correct elevations are obtained in non-overlaid regions.
4401.. . I . . . A . . . ...I .. . . 2800 3000 3200 3400 3600 3800 4000
~~~~
\\"
480 460
"
'
'
Figure 5: Crosstrack elevation profile filtered with 3x3 window sue
560 540
'
440
2820
2840
7
440 2800
3000
3200
7 3400 3600 3800 Crosstrack [ m ]
4000
4200
4
'
'
4200
Crosstrack [m]
Without signal mixing, a symmetric peak in groundrange pictured in Figure 3 would be tilted towards the radar in slantrange-coordinates as in Figure 4.The tilt is such that it would be exactly offset when going from slantrange to groundrange, thereby restoring the original form again. But in reality the peak in slantrange will be more symmetric because one can obtain only one elevation for each range pixel. Transforming the more symmetric peak back to groundrange, however, results in a net tilt away from nadir. 5 8 0 1 ' '
5
r_(
These layover areas are not masked with DOSAR data for two reasons:
Due to the high range resolution of the 200 MHz range-chirp data, very small surface objects are resolved which tend to have radar facing slopes greater than the off-nadir angle. Therefore, we also want to preserve as well as possible the form o f these objects in the surface model.
+ Ay) 2 + ( p j - ( zref + Az) ) 2
2860
2880
2900
2920
Groundronge [m]
Figure 3: Peak in groundrange-coordinates
Figure 6: Crosstrack elevation profile filtered with 27x27 window size
7.CONCLUSION It was found that one has to solve the Interferometric RangeDoppler Equations for a correct description of the interferometric imaging geometry. Doing so one also obtains a model o f the alongtrack interferometric baseline which may be adjusted with offsets from complex image correlation in the case of crossed-tracks geometry. Using the range-offsets in combination with the suggested phase calibration approach one would be less dependent on ground control points; a situation especially convenient in spaceborne interferometry for large area investigations. High resolution interferometric DSM generation requires separate filtering of the DSM crosstrack-coordinate. The window size has to be adjusted to the local ambiguities which occur as a result of layover.
\
460
4 4 0 .
3760
. , . . . 3780 3800 I
.
REFERENCES .
.
.
3820
3840
3860
3880
slontronge [m]
[ 11
\N.Goblirsch et al: Accuracy of Interferometric Elevation Models Generated from DOSAR Airborne Data, IGARSS'95 vol. 1, pp. 770 -
[2]
N. Faller, E. Meier: First Results with the Airbome Single-Pass DCC; AR Interferometer, IEEE Transactions on Geoscience and Remote Sensing, vol. 33, no. 5, 1995, pp. 1230-1237
[3]
\K. Goblirsch: Error Sources in Optimization of Parameters of Airborne Interferometric Synthetic Aperture Radar. Thesis, Department of Geography, University of Zurich 1996, in press
Figure 4: Peak in slantrangecoordinates This happens with all surface features having sufficiently steep radar facing slopes which is the case for mostly small-scale objects. Thus when we straighten out the crosstrack-coordinate relative to slantrange so much that local ambiguities disappear we will obtain a monotonic slantrange which removes the tilt in the DSM. This has been done with 27x27 moving average filter of y in Figure 6. Smaller window sizes re-
774
