ISPRS Journal of Photogrammetry and Remote Sensing 86 (2013) 11–20

Contents lists available at ScienceDirect

ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs

Radargrammetric registration of airborne multi-aspect SAR data of urban areas Michael Schmitt ⇑, Oliver Maksymiuk, Christophe Magnard, Uwe Stilla Photogrammetry and Remote Sensing, Technische Universitaet Muenchen (TUM), Arcisstr. 21, 80333 Munich, Germany Remote Sensing Laboratories, University of Zurich, Winterthurerstr. 190, 8057 Zurich, Switzerland

a r t i c l e

i n f o

Article history: Received 23 November 2012 Received in revised form 17 July 2013 Accepted 4 September 2013

Keywords: Synthetic Aperture Radar (SAR) Multi-aspect Airborne Registration Georeferencing Gauss-Helmert model

a b s t r a c t In this paper, the registration of decimeter-resolution airborne multi-aspect SAR (MASAR) data of inner city areas by application of the radargrammetric rager equations is investigated. The geometrical model is adapted to linear flight trajectories and zero-Doppler processed SAR data, whereas the observed trajectory parameters are adjusted using a strict Gauss–Helmert model and known ground control points. The significance of the estimated corrections is examined and the most suitable set of free parameters is determined. Finally, the methodology is applied to real test data of an airborne campaign over the city of Munich, Germany, and the feasability of the proposed radargrammetric registration method is shown. Ó 2013 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS) Published by Elsevier B.V. All rights reserved.

1. Introduction The coregistration of decimeter-resolution SAR imagery showing densely built-up inner city areas is still an unsolved problem due to the inherent side-looking imaging geometry (Stilla, 2007). While it is highly desirable to combine or fuse SAR images acquired from different aspects in order to replace zones of radar shadow with data taken from a different viewing angle (Schmitt et al., 2011), the corresponding layover effect that causes elevated objects to be displayed ‘‘collapsing’’ towards the sensor does not allow conventional image registration methods that are typically designed for orthographically projected imagery or imagery acquired from only slightly differing viewing angles (Zitova and Flusser, 2003). Therefore, only some first approaches that either require preliminary georeferencing (Middelmann et al., 2003; Soergel et al., 2004), or imagery with lower resolutions in the meter-domain (Curatolo et al., 2003; Dell’Acqua et al., 2004) have been published. Gernhardt et al. (2012) even avoid the whole image registration problem by fusing only geocoded persistant

⇑ Corresponding author at: Photogrammetry and Remote Sensing, Technische Universitaet Muenchen (TUM), Arcisstr. 21, 80333 Munich, Germany. Tel.: +49 89 289 22672; fax: +49 89 289 23202. E-mail addresses: [email protected] (M. Schmitt), [email protected] (O. Maksymiuk), [email protected] (C. Magnard), [email protected] (U. Stilla). URL: http://www.pf.bv.tum.de (U. Stilla).

scatterer point clouds derived from multi-temporal InSAR data stacks. A field related to multi-aspect SAR (MASAR) image registration is the combination of SAR and optical image data as both kinds of data differ in their imaging geometry due to their different imaging principles. Galland et al. (2005) provide an interesting solution for this situation based on both the photogrammetric and radargrammetric sensing equations. Building upon the insight that a fusion of decimeter-resolution multi-aspect SAR data of urban areas is only reasonably possible in object space, this paper presents a method that combines the wellknown calibration of the SAR flight geometry using one or more known ground control points (Wu and Lin, 2000; Sohn et al., 2005; Leberl, 1990) with a least squares approach formulated as a strict Gauss–Helmert model (Neitzel and Petrovic, 2008; Neitzel, 2010; Lenzmann and Lenzmann, 2004) and the simultaneous adjustment of an arbitrary number of acquisitions of the same urban scene taken from arbitrarily oriented flight tracks. The proposed approach does not use any approximations or assumptions with respect to the projection of the imagery to be registered; also it does not require a preliminary orthoprojection or similar preprocessing step. Although we have designed the mathematical model for linear flight tracks and zero-Doppler processed data, the algorithm can straight-forwardly be extended to more general imaging geometries. Therefore, the contribution of this article can be summarized as follows:

0924-2716/$ - see front matter Ó 2013 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS) Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.isprsjprs.2013.09.003

12

M. Schmitt et al. / ISPRS Journal of Photogrammetry and Remote Sensing 86 (2013) 11–20

 A new solution of the well-known radargrammetric rangeDoppler calibration problem utilizing a least-squares adjustment of observations only is introduced. Therefore, only one ground control point is necessary.  This adjustment is formulated in a strict Gauss–Helmert model, thus providing the correct solution even in the case of large residuals.  The model is described for the simultaneous calibration of an arbitrary number of different acquisitions of the scene, which enables the precise radargrammetric registration of multiple aspect SAR data of urban areas. After an introduction to the geometrical model and the parameter estimation framework a thorough analysis of all involved observations is carried out in order to determine which flight parameters are supposed to be considered correctable and which parameters should be considered constant. Also, a proper weighting of the observed variables by employment of a reasonable stochastical model is discussed. The findings of these preliminary investigations are then applied to a simultaneous adjustment of multiple SAR aspects. Finally, the applicability of the proposed approach is qualitatively and quantitatively discussed.

2. Radargrammetric calibration in a strict Gauss–Helmert model Papers explaining the radargrammetric calibration of SAR imaging geometries using one or more ground control points have been published for decades (Schreier, 1993). While it is not uncommon that least-squares based approaches are employed, e.g. (Hellwich and Ebner, 2000; Raggam and Gutjahr, 2000; Yue et al., 2008), to our knowledge no solutions adjusting only observations in a strict Gauss–Helmert model were proposed to this day. In this article, we show the formulation of the well-known range-Doppler approach in the framework of this sophisticated model. The main advantage of this estimation method is that it does not depend on the introduction of any unknown parameters and just seeks to estimate corrections that adjust the observations in a least squares sense (Mikhail, 1976, pp. 137 ff.). One of the core advantages of this implementation is the fact that an arbitrary number of flight tracks can be calibrated simultaneously without any loss in redundancy: Usually, the parameters describing the flight track are considered as unknowns in the calibration, such that each track would introduce new unknowns and therefore reduce the available redundancy. Due to that, the number of necessary ground control points and the number of tie points that link the images one to another increases. In contrast to that, in our formulation for every flight track that is included in the estimation, the redundancy is even enhanced, because only additional conditional equations without any unknowns are introduced. An additional benefit is that in a simultaneous adjustment, all flight tracks will be linked via the ground control points that have jointly been used for the calibration. Thus, a precise combination of multi-aspect SAR and InSAR data of urban areas is enabled, an otherwise non-trivial task.

First, the range and azimuth coordinates of any ground control point measured in the scene can be derived from its pixel position using the linear relations

t ¼ t 0 þ dtðr  1Þ;

ð1Þ

and

R ¼ R0 þ dRðc  1Þ;

ð2Þ

respectively. Here, t0 denotes the time at which the first azimuth bin was imaged, whereas R0 is the near-range value of the system given in slant range geometry. dt and dR denote the pixel spacing in azimuth (time) and range direction, and r and c are the row (azimuth) and column (range) coordinates of the control point, counted from the center of the top left image pixel. Using (1), the three-dimensional sensor position during each azimuth bin can be derived:

sðtÞ ¼ s0 þ v  t:

ð3Þ

Analogous to t0, s0 describes the sensor position during the acquisition of the first azimuth bin, and v describes the sensor velocity. Eqs. (2) and (3) are then put into the range and Doppler equations

R  kp  sðtÞk ¼ 0;

ð4Þ

vðp  sðtÞÞ ¼ 0;

ð5Þ

where p is the 3D position of the ground control point in world coordinates. Again it is important to note that (5) is a simplified Doppler equation that only holds for zero-Doppler processed data. It can, however, easily be extended to other types of SAR processing without loss of generality for the reasoning of this article. If (4) and (5) are set up for every ground control point in every scene, a non-linear system of conditional equations is constructed. Since all included parameters (see Table 1) are known from the flight navigation control and the SAR focusing, respectively, we propose to model the correction of the sensor geometry by a Gauss–Helmert model without unknowns (Mikhail, 1976, pp. 137 ff.). 2.2. Parameter estimation in the strict Gauss–Helmert model As mentioned before the system of conditions consists of the range and Doppler Eqs. (4) and (5). They depend not only on the flight track, but also on the positions of the ground control points and their corresponding image positions. Therefore we have the objective

0

1 ^ f 1 ðbÞ C ^ ¼B B .. C ¼ 0; fðbÞ @ . A ^ f n ðbÞ

ð6Þ

^ ¼bþm ^ are the corwhere n is the number of flight tracks and b rected observations. The vector of conditional equations for a single flight track i 2 {1, . . . , n} is

1 ð1Þ ^ fi;r ðbÞ C B B f ð1Þ ðbÞ ^ C C B i;d C B B . C ^ f i ðbÞ ¼ B .. C; C B B ðmÞ ^ C B fi;r ðbÞ C A @ 0

2.1. Recapitulation of range-Doppler geometry As mentioned before, SAR imaging geometry calibration routines are usually based on the well-known range-Doppler equations (Leberl, 1990), which – in our case – are combined with the assumption of linear flight trajectories and zero-Doppler processed data. This is justified for airborne SAR data with limited swath length and an appropriate focusing of the raw radar data (Magnard et al., 2012).

ð7Þ

ðmÞ ^ fi;d ðbÞ

ðjÞ

ðjÞ

with m being the number of ground control points and fi;r and fi;d representing (4) and (5) for flight track i and ground control point

13

M. Schmitt et al. / ISPRS Journal of Photogrammetry and Remote Sensing 86 (2013) 11–20 Table 1 Observation types. Observation

Unit

Order of magnitude

Standard deviation

Sensor parameters Sensor starting position, s0 = [s0X, s0Y, s0Z]T Sensor velocity, v = [vX, vY, vZ]T Starting time, t0 Azimuth pixel spacing (time), dt Near range, R0 Range pixel spacing (distance), dR

m m=s s s m m

106 101 103 104 103 101

0.1 0.02 0 0 0.02 0

Ground control point measurements GCP position in WGS84, p = [pX, pY, pZ]T GCP image measurements, r, c

m pixel

106 104

0.03 0.5

j 2 {1, . . . , m}. Eq. (6) is the objective we want to achieve by calculating the correction vector m^. In general, for the application of a Gauss–Helmert model, the conditional equations have to be linearized, i.e. a Jacobian matrix containing the partial derivatives of the conditional equations with respect to all the corrections is needed:

B¼

@fðm Þ ; @m

ð8Þ

f denotes the vector of conditional equations, m the vector of corrections. A visualization of the structure of the Jacobian matrix for a typical configuration is illustrated in Fig. 1. We seek to estimate these corrections m from the linearized set of conditional equations

Bm^ þ w ¼ 0;

ð9Þ

based on the least squares principle. In most of the common literature on least squares estimation, the linearization is only carried out at the position m0 = 0, leading to the vector of contradictions

w ¼ fðm 0 Þ:

ð10Þ

We, however, intend to follow the approach of Lenzmann and Lenzmann (2004) and thus to renew the linearization for every update of m0, which yields the actual vector of contradictions

w ¼ Bm 0 þ fðm 0 Þ;

ð11Þ

that has to be recalculated in every iteration of the minimization of the objective function. If (11) is used instead of (10), this is called ‘‘strict’’ Gauss–Helmert model in the rare literature about the topic and was shown to equal common Total Least Squares solutions (Neitzel and Petrovic, 2008; Neitzel, 2010). The afore-mentioned objective function is in both cases defined as

^T

X¼m

m  2k^ T ðBm^ þ wÞ;

^ Q 1 bb

ð12Þ

where Qbb is the cofactor matrix of the observations that can be derived from their covariance matrix and is used for the weighting of the observations. If Qbb is the unity matrix, an unweighted least ^ is the vector of correlates and calculated squares solution is found. k by

^ ¼ ðBQ BT Þ1 w: k bb

ð13Þ

^ this finally leads to the If (12) is minimized with respect to m^ and k, least squares estimator for the vector of corrections:

^ m^ ¼ Q bb BT k:

ð14Þ 6

^ Þ is fulfilled to a satisfying precision (e.g. 10 ), the iteration If fðm ^ can be aborted and the final solution is found. Otherwise, m 0 ¼ m is used for the next iteration and a new set-up of (8) and (11). ^ , the adjusted obserWith the estimated vector of corrections m vation vector

^ ¼ b þ m^; b

ð15Þ

can be derived easily. The main benefit of the strict formulation of the Gauss–Helmert model is that the conditional equations will be fullfilled to almost arbitrary precision after the adjustment. A table illustrating the benefit numerically is presented in Section 4.4. 2.3. Implementation of a step size parameter If the conditional equations (9) are fulfilled, the solution of the objective function (12) is given by T

^ Xmin ¼ m^T Q 1 m^ ¼ k w:

ð16Þ

Obviously the solution depends on the vector of contradictions. From our experience, however, with the applied mathematical model and the described optimization, chances are high that the solution is oscillating. Therefore, we propose the utilization of a suitable step size parameter such that m0 in (11) is replaced by av0 constrained to a P 0. The objective function to determine the optimal step size parameter a for updating the contradictions in every iteration is given by

kwk22 ! min :

ð17Þ

This objective enables to construct a solution in which the norm of the iteratively updated contradiction vectors is a monotonically decreasing sequence in the number of iterations. Objective (17) is a simple one-dimensional optimization problem with the analytical solution T

a¼

f p ; pT p

ð18Þ

Fig. 1. Exemplary structure of the Jacobian matrix for the case of just a single flight track with two ground control points (p1, p2) imaged. The magnitude of the partial derivatives is shown from white (no derivative) to black (very large derivative).

14

M. Schmitt et al. / ISPRS Journal of Photogrammetry and Remote Sensing 86 (2013) 11–20

where p = Bm0 given that the Jacobian B and the conditional equations f are evaluated at m0 and assumed to be independent of m and therefore constant during one iteration. The latter is only fulfilled if B and f are linear in m but in practice this procedure works very well leading to a solution in only a few iterations. The reason for the introduction of the step size parameter a constraining the iterative updates is that the linear system (11) shows a significant sensitivity according to small perturbations in B and m. As stated in (Golub and Loan, 1996, pp. 80 ff.) the condition number j(B) = kmax/kmin, where k are the singular values of B, is a measure for this kind of sensitivity. It indicates the amount of amplification of w for a small perturbation of m. In our case it is especially the derivation of the Doppler equation with respect to the azimuth (time) pixel spacing dt which introduces significant imbalances in the singular values of the Jacobian matrix (cf. Fig. 1) leading to j(B)  106. 2.4. Covariance matrix of observations As mentioned in Section 2.2 already, it is possible to accomplish a weighted least squares solution using the cofactor matrix Qbb derived from the covariance matrix Cbb of the observations by

Q bb ¼

1

r20

Cbb :

ð19Þ

A priori, the so-called variance factor r20 is usually set to 1. Its actual value is estimated after the adjustment (see Section 2.5). In general, a weighting of the observations can always be considered reasonable, especially if they are very heterogeneous in their nature. In the functional model employed for this study, observations of different units and different orders of magnitude are combined, e.g. pixel coordinates, world coordinates, time measurements, and velocities (see Table 1). Furthermore, the stochastical model can also be used to exclude certain observations from the least squares adjustment if they are to be considered constant and their variance is set to 0 (or rather a small number  > 0 for computational reasons). Different approaches for the stochastic model are compared in the experimental Section 4. 2.5. Posterior statistical evaluation Based on variance–covariance propagation, it is possible to deduce the covariance matrix of the corrections and the adjusted observations after the parameter estimation. The covariance matrix of the corrections is found to be

 1 ^ 20  Q bb BT BQ bb BT Cm^m^ ¼ r BQ bb ;

ð20Þ

the covariance matrix of the adjusted observations to be

   1 ^ 20  Q bb  Q bb BT BQ bb BT Cb^b^ ¼ r BQ bb ;

ð21Þ

with the posterior variance factor

r^ 20 ¼

^ m^T Q 1 bb m r

;

ð22Þ

3. MEMPHIS test data For a real-data based assessment of the methodology proposed in this paper, a multi-aspect SAR dataset acquired by the German MEMPHIS sensor over the city of Munich, Germany, in 2011 has been used. 3.1. System description The MEMPHIS sensor was first introduced by the Institute of High Frequency Physics and Radar Technology of the FGAN (Forschungsgesellschaft für angewandte Naturwissenschaften, now part of Fraunhofer Society) in 1998 (Schimpf et al., 2002). The acronym stands for Millimeterwave Experimental Multifrequency Polarimetric High-resolution Interferometric System. It operates at 35 GHz and 94 GHz (Ka band and W band), respectively. Being an experimental, modular and removable system, it is typically mounted on a C-160 Transall airplane of the German Armed Forces. Using various antenna shapes and configurations, many different SAR modes are possible: single-pass multi-baseline cross-track interferometry with four receiving antennas, dual-pol circular or linear polarimetry and monopulse for MTI applications. The radargrammetric registration presented in this paper can for example be considered as a pre-processing step for urban surface model or building reconstruction as proposed in (Schmitt and Stilla, 2011) or (Maksymiuk and Stilla, 2012). 3.2. Raw data processing For the 2011 campaign, the navigation data were acquired with a differential GPS (dGPS) system and a precise inertial navigation system (INS). The GPS, INS and SAR systems were synchronized through event markers, second markers and an NMEA link. After processing with the commercial software Inertial Explorer, the navigation data were finally smoothed with a Kalman filter to avoid small variations in the millimeter range, which would introduce artifacts in the focused SAR data. The lever arms between the dGPS antenna, the INS and the SAR antennas fixed in operating position were measured using a terrestrial surveying method with a few centimeters accuracy. Since MEMPHIS is a stepped-frequency radar system, it successively transmits 8 chirps of 200 MHz bandwidth with a 100 MHz carrier frequency shift between each other, thus building together a full bandwidth of 900 MHz. The raw data from each chirp are first focused in range using chirp replica with the conventional matched filtering technique. The range compressed data are then put together in the frequency domain through an algorithm based on (Lord, 2000). The azimuth compression is performed with the Extended Omega-K algorithm (Reigber et al., 2006), resulting in a zero-Doppler slant range geometry. The block diagram in Fig. 2 shows the processing chain of the algorithm. The motion compensation is a critical step for reaching high geolocation accuracy. The first order motion compensation (moco) is achieved as follows:  The navigation data are upsampled to the PRF1 rate.  A linearized track is defined using a least squares method on the position data in global Cartesian coordinates. A constant linearized velocity is defined.  For each echo, the projection of the vector linking the real and the linearized track onto the mid-range line of sight is evaluated. The phase and position of the range-compressed (RC) data are corrected according to this value.

r is defined as the redundancy of the equation system, which equals the number of conditional equations for the Gauss–Helmert model without unknowns. It has to be noted that the variance factor also contains an assertion about the a priori stochastical model that has been employed by the relation given already in (19): Using (22), the ‘‘true’’ covariance matrix of original observations can be calculated by

b bb ¼ r ^ 20  Q bb : C

ð23Þ

^ 20 is to 1, the better were the stochastic This means: The closer r assumptions about the original observations.

1

PRF: pulse repetition frequency.

M. Schmitt et al. / ISPRS Journal of Photogrammetry and Remote Sensing 86 (2013) 11–20

15

Fig. 2. Block diagram of the Extended Omega-K algorithm (moco: motion compensation, (I) FFT: (Inverse) Fast Fourier Transform, SLC: Single Look Complex).

Fig. 3. Quadruple corner reflector consisting of four trihedrals.

Table 2 MEMPHIS test data. System parameters Wavelength Off-nadir angle H Range pixel size dR Azimuth pixel size dAz

8.5 mm (Ka band) 60° (Track 14), 65° (Tracks 15–18) 0.17 m 0.05 m

Track 14 Heading angle Flying altitude

340° 768 m

Track 15 Heading angle Flying altitude

70° 714 m

Track 16 Heading angle Flying altitude

250° 709 m

Track 17 Heading angle Flying altitude

160° 714 m

Track 18 Heading angle Flying altitude

340° 712 m

 The RC data are interpolated in azimuth direction according to the constant linearized velocity, getting a regular (constant) spatial azimuth sampling. The second order motion compensation consists of the following process:

Fig. 4. Flight track configuration and overlapping scenes of the 2011 MEMPHIS MASAR campaign (optical image of downtown Munich Ó 2012 Google).

The motion compensation handles objects lying at the altitude of the reference surface without error. In case of large deviations between the real and linearized flight paths, geolocation errors and focusing degradation get worse with growing distance of the imaged objects from the reference height. As the test site is a rather flat area, this does only concern very high buildings. 3.3. Available test data The multi-aspect dataset that has been used for the experiments in this paper is comprised of four aspects from anti-parallel and orthogonal trajectories, and one aspect acquired from a different flying altitude. A summary of the relevant scene parameters can be found in Table 2, a sketch of the track configuration in Fig. 4. Since corner reflectors usually have a limited bandwidth and therefore can only be seen from a specific flight direction, we have employed two quadruple corner reflectors, each consisting of four adjacent trihedral corners (see Fig. 3). 4. Calibration results 4.1. Calibration of a single aspect

 The vector linking the real and linearized track is projected onto the line of sight of each range sample (using a flat reference surface). The difference between this projection and the one applied in the first order motion compensation is used for correcting the phase and position of the data.

The first series of experiments conducted aimed at the determination of the observed parameters that are supposed to be corrected in opposition to the ones that are to be considered error-free. Therefore, the calibration routine proposed in Section 2

16

M. Schmitt et al. / ISPRS Journal of Photogrammetry and Remote Sensing 86 (2013) 11–20

is first applied to all available aspects separately. After the adjustment, the relation between correction and the respective standard deviation, the so-called t-value

  m ^ t ¼  ;

rm^

ð24Þ

is determined for every parameter with the standard deviation having been extracted from the covariance matrix of corrections. The higher the t-value is, the more significant is the corresponding correction. Therefore, this relation is intended to deliver a feeling for the importance of each correction, i.e. it can be used to determine which of the observed parameters are prone to correction and which can be considered constant. Besides the t-values, however, also a priori knowledge about observations and their physical meaning have to be considered. In general, they can be classified in three different groups: trajectory parameters (s0 and v), system parameters (t0, dt, R0, and dR) and ground control point measurements (p and (r, c)). The final goal of each radargrammetric calibration is always a refinement of the trajectory parameters based on ground control points. Therefore, both these groups are intended for correction by the adjustment by principle. Therefore, the investigations of this chapter are mainly meant for an evaluation of the correctability of the system parameters. Among them, t0 is the time stamp linking the initial position of the aircraft s0 and the acquisition of the first azimuth bin. It can purposefully be considered constant if a local time synchronization is used in the sense that both image acquisition and navigation processing start at t0 = 0. 4.1.1. Calibration with unweighted observations In Fig. 5, the averaged t-value results of a separate calibration of all available aspects are shown for the case of an unweighted least squares estimation. It can nicely be seen that most parameters show a comparable significance. The main exceptions are found to be the pixel size in azimuth direction dt as well as the pixel size in range direction dR, where the corrections were estimated close to 0. In addition to them, also t0 and the GCP pixel position in azimuth (r) are corrected by just about half the standard deviation of the respective correction. The next test used the same configuration (separate calibration of all five available aspects) but with all parameters related to the ground control points being fixed. The intention behind this experimental set-up was to determine the distribution of the corrections if all errors are supposed to occur in the flight trajectory parameters and not the control measurements. Although this is not a valid assumption in reality, we expected to get additional information and a better feeling for the significance of the involved parameters. Fig. 6 shows the related t-values. Again, dt and dR are the least

Fig. 6. Median t-values for the separate unweighted calibration of 5 aspects with fixed ground control related parameters. The black error bars indicate the standard deviations of the t-values.

significant parameters (in fact they were almost not corrected at all). Also, again t0 is found to be only corrected to a fraction of the correction’s standard deviation. 4.1.2. Calibration with weighted observations From the findings of these first two experiments, a reasonable weighting of the observations by introduction of a stochastical model as described in Section 2.4 has been added to the estimation. The standard deviations used for the stochastical model are listed in Table 1. Besides the system parameters dt and dR that have been found not to be prone to correction in the previous tests, also the system parameter t0 has been considered a fixed value in this model due to the reason explained in Section 4.1. All other parameters have been weighted based on considerations about the accuracies of the navigation systems and the flight processing. Therefore, the sensor position was defined as the least accurately known parameter, whereas the sensor velocity was considered quite precise. The image coordinate of the ground control point were just measured manually, yielding an accuracy of 0.5 pixel at best. The GCP world coordinates have been measured using a differential GPS and the German SAPOS2 reference station system, so that an accuracy of about 3 cm can be assumed. Fig. 7 shows the resulting t-values calculated from the averaged results of the separate weighted adjustments of all available aspects. It seems interesting that all GCP related observations are comparably significant, while for the system related parameters R0 and s0Z are found to be the most important. It should be noted, however, that all calculations were carried out in the WGS84 reference system, so that none of the coordinates can directly be related to an interpretable dimension like flying altitude. The comparably large standard deviations of the t-values are also noteworthy. They indicate that the respective parameter is not equally significant in all aspects. From experience and detailed analysis of the data, we found that the high-altitude aspect 14 only needs very little corrections, whereas aspect 15 seems to contain more a priori mis-calibration errors, which can also be figured from Table 3. 4.2. Joint calibration of multiple aspects

Fig. 5. Median t-values for the unweighted separate calibration of 5 datasets acquired from different viewing angles. The black error bars indicate the standard deviations of the t-values.

After looking at the separate calibration of the available datasets, the algorithm is applied to all aspects simultaneously for a detailed assessment of the correction properties in operational MASAR cases. The main advantage of such a joint calibration is that a kind of network between all the parameters of the different flight tracks and the ground control points is established, such that any inconsistencies between them are efficiently eliminated by the adjustment. 2

Satellite Positioning Service of the German State Survey.

17

M. Schmitt et al. / ISPRS Journal of Photogrammetry and Remote Sensing 86 (2013) 11–20

Fig. 7. Median t-values for the weighted separate adjustment of the 5 aspects. The black error bars indicate the standard deviations of the t-values.

4.2.1. Calibration with unweighted observations The first experiment for the simultaneous calibration of all MASAR datasets was carried out in the vein of the single aspect experiments of Section 4.1 – with unweighted observations. As Fig. 8 shows, dt and dR remain almost uncorrected also in this case. Interestingly, the error bars showing relatively high standard deviations of the t-values, especially for the sensor position and velocity parameters, indicate a large spread around the depicted medians already for the unweighted case. Again, this means that large differences between the t-values in the involved tracks occur: Looking at s0X, for example, reveals that the estimated corrections are 1.5 times their standard deviation on average. However, there must also be single cases where the correction is less significant – as the large standard deviations of the t-value indicate. 4.2.2. Calibration with weighted observations In order to finalize the experiments, we have carried out the weighted adjustment as described in Section 4.1.2, but this time for all aspects simultaneously. The resulting t-values are shown in Fig. 9. Again, the most significant corrections are found to be among the sensor positions (particularly s0X and s0Z, the near range value R0 and the ground control point parameters (particularly r and pX). 4.3. Choice of adjustable parameters It has been shown that the system-inherent parameters linked to the pixel size are known to a very high precision so that almost no correction is estimated for them even for unweighted observations. Therefore, we have decided to consider them constant and leave them out of the least squares adjustment in the weighted case. In addition to that, we canceled also the t0-parameter. On the one hand, this choice was justified by relatively low significance values, on the other hand t0 is by its nature defined to be the time stamp during which the first azimuth bin is being acquired. Since the synchronization between the navigation system and the MEMPHIS sensors is very precise with an error below

Fig. 8. Median t-values for the unweighted joint calibration of all 5 aspects. The black error bars indicate the standard deviations of the t-values.

Fig. 9. Median t-values for the weighted joint calibration of all 5 aspects. The black error bars indicate the standard deviations of the t-values.

¼ 8:3  105 s , t0 can validly be defined to be 0. The remaining variances have been set according to reasonable approximations based on our knowledge of the available navigation system quality as well as the raw SAR data focusing. In the end, we came to the conclusion that the relevant system parameters to be corrected have to be the sensor position, the sensor velocity, and the near range offset. In addition to that, all ground control point related parameters are clearly observations in a narrow sense, and thus have to be coupled with realistic variances as well.

1 8PRF

4.4. Quality of calibration The quality of the calibration can best be assessed if the remaining error between the image of a corner reflector and its projection into the image is measured. While the available test data contained errors up to some pixels in non-calibrated geometry, the calibration method proposed in this article causes this projection error to vanish completely (see Table 3 and Fig. 12). More evidence of the calibration quality can also be found in the following section, where it is evaluated for the registration of multiple aspects.

Table 3 Projection of a corner reflector into the five images, before the calibration (left) and after the calibration (right). Note that after the calibration also the image measurements of the corner have been corrected to subpixel precision. The small remaining error in aspect 14 is mainly due to transformations between UTM32 and WGS84 coordinate systems. Track

14 15 16 17 18

Before calibration

After calibration

Measured position [r, c]

Projected position [r, c]

Measured position [r, c]

Projected position [r, c]

15461 18042 12825 18453 12502

15460.58 18040.93 12826.43 18455.09 12502.34

15460.65 18042.32 12825.11 18452.64 12502.14

15460.66 18042.32 12825.11 18452.63 12502.14

593 1527 2182 3053 980

592.85 1524.82 2180.44 3053.03 977.51

593.34 1526.97 2182.24 3053.32 979.57

593.34 1526.97 2182.24 3053.32 979.57

18

M. Schmitt et al. / ISPRS Journal of Photogrammetry and Remote Sensing 86 (2013) 11–20

Table 4 Averaged L1-norm of contradictions of all aspects for different adjustments and methods. It can be seen that the strict Gauss–Helmert model completely removes all contradictions. Adjustment

Unweighted separate calibration Unweighted separate calibration with fixed GCPs Weighted separate calibration Unweighted joint calibration Weighted joint calibration

Gauss–Helmert Classic

Strict

0.00945749 0.02112841

0.00000000 0.00000000

0.00348838 0.02431621 0.01216487

0.00000000 0.00000000 0.00000000

Another proof for the accuracy of the calibration can be taken from the fulfillment of the conditional Eqs. (6) and (7), respectively, as any remaining contradiction indicates that the underlying model does not fully fit the observations. Table 4 shows the L1-norm of

the final contradiction vectors (10) resulting from the five experimental setups of Sections 4.1 and 4.2 for the conventional and the strict formulation of the Gauss–Helmert model. It can clearly be assessed that only the strict Gauss–Helmert model is able to remove all contradictions and therefore to fulfill the conditional equations to a satisfying precision. 5. Multi-aspect image registration Since the final goal of this paper is the radargrammetric registration of MASAR data via range-Doppler geometry and the object space represented by the locally flat ellipsoid surface, both the calibration quality as well as the feasability of the proposed approach for MASAR registration is best assessed using real MASAR datasets. As has been shown, the proposed method enables the simultaneous calibration of multiple airborne SAR images acquired from arbitrary flight paths. After the joint calibration, the images can conveniently be projected onto a reference plane on the locally flat

Fig. 10. MASAR data (in slant range geometry, pixels approximately squared) acquired by the MEMPHIS sensor. The images show the area around Technische Universität München and Alte Pinakothek in Munich, Germany. The red rectangle depicts the scene which is used to explain the calibration and registration in detail later. Note the different viewing directions indicated by the arrow in the upper left corner. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. Images obtained from fusing all five aspects shown in Fig. 10. Left is the uncalibrated case and on the right is the result of our proposed method. As can be seen the details are much better represented due to the calibration.

M. Schmitt et al. / ISPRS Journal of Photogrammetry and Remote Sensing 86 (2013) 11–20

19

Fig. 12. Comparison of the reflector representation before (a) and after (b) calibration. Obviously the reflector appears much more focused in the calibrated case. Image (c) shows a part of a road with some vehicles and a hoarding along the left side of the road. It can be seen that the level of detail is increased significantly in the calibrated case (d) making it easier to recognize the different objects.

Fig. 13. Maximum intensity fusion of anti-parallel tracks 14 and 17. The image shows the layover extent of the State Museum of Egyptian Art in Munich, Germany, from two opposing viewing directions. Note that both the front and the back facades are clearly visible.

ellipsoid surface such that they are registered with respect to all structures lying on this plane. For the projection, the range-Doppler equations (4) and (5) are combined with the ellipsoid equation

p2X þ p2Y ða þ hÞ

2

þ

p2Z 2

ðb þ hÞ

 1 ¼ 0;

ð25Þ

where a and b are the semi-major and semi-minor axes, respectively, and h is the height of the corresponding pixel above the reference ellipsoid (Schwaebisch, 1998). This system of three equations can be solved for p = [pX, pY, pZ]T using an iterative optimization algorithm as for example the trust-region dogleg method (Powell, 1970). If the images are resampled to a common grid afterwards, they can all be displayed in the georeferenced world coordinate system. For visualization purposes, we have decided to display the average intensity over all georeferenced aspects. Fig. 10 shows a dataset consisting of five aspects with a red rectangle depicting the area used for an exemplary image registration. This area is of special interest because of its flat topography, since the correct projection critically depends on the height h of each pixel. This means that for flat but built-up urban terrain only structures located on the reference plane are displayed properly, e.g. streets, paths, lot boundaries or suitably placed corner reflectors. In Fig. 11, a before-and-after comparison is shown for some ground-level point scatterers including one of the corner reflectors that were used as input for the calibration. It can clearly be seen that the calibration leads to a severely better focused overlay image of the involved aspects. The calibration leads to a proper fusion of image pixels resulting in a finer representation of objects. This is a prerequisite for 3D reconstruction but is also advantageous in automatic object detection and recognition. An even closer look on the details is shown in Fig. 12. As mentioned above, the reflector was also used for calibration and with Fig. 12(a) and (b) one can easily verify that the calibration successfully attains a proper image registration. Beside the reflector there are other interesting structures in the scene, namely the road and the vehicles which are approximately at the same height as the referenced reflector. Fig. 12(c) and (d) shows a comparison on this closeup and the increase in the level of detail of the structures. It should be mentioned that at the time of image aquisition there was a hoarding along the left side of the street. In addition to that, Fig. 13 depicts the layover characteristics for a building with respect to different aspects acquired from flight tracks with different heading angles. Although this might seem

20

M. Schmitt et al. / ISPRS Journal of Photogrammetry and Remote Sensing 86 (2013) 11–20

disturbing at first glance, the different displacements can be used as valuable input information for radargrammetric image analysis: For example some papers have already been published aiming at the reconstruction of building and bridge heights by exploiting their layover extent measured in ground range imagery (Soergel et al., 2009; Wegner and Soergel, 2008). If, however, a detailed urban surface model of the scene were available, it could be used to create a true orthoprojected version of the available multi-aspect SAR images, where shadow areas of one aspect could be filled with data from a complementary aspect (Schmitt and Stilla, 2010).

6. Conclusion In this paper, a method for the simultaneous calibration of multi-aspect SAR data of urban scenes using a least squares approach formulated as strict Gauss–Helmert model is proposed. After the flight parameters that are supposed to be corrected are determined and distinguished from the constant system parameters, the algorithm is successfully applied for the radargrammetric registration of MASAR data. It has been shown that the projection of jointly calibrated SAR images can be fused, e.g. by simple averaging, after they have been projected on the locally flat Earth ellipsoid surface. Finally, the presented approach both provides an interesting perspective for the generation of true SAR ortho-images and is also a promising pre-processing step for radargrammetric image analysis. Acknowledgements The authors would like to thank Thorsten Brehm and Dr. Stephan Stanko of the Fraunhofer Institute for High Frequency Physics and Radar Technology FHR for providing the MEMPHIS test data. References Curatolo, F., Dell’Acqua, F., Gamba, P., Lisini, G., 2003. Multiangle high resolution SAR image co-registration: a feature-based approach. In: Proceedings of Tyrrenian International Workshop on Remote Sensing, pp. 667–673. Dell’Acqua, F., Gamba, P., Lisini, G., 2004. Coregistration of multiangle fine spatial resolution SAR images. IEEE Geoscience and Remote Sensing Letters 1 (4), 237– 241. Galland, F., Tupin, F., Nicolas, J.-M., Roux, M., 2005. Registering of synthetic aperture radar and optical data. In: Proceedings of IEEE International Geoscience and Remote Sensing Symposium, pp. 3513–3516. Gernhardt, S., Cong, X., Eineder, M., Hinz, S., Bamler, R., 2012. Geometrical fusion of multitrack PS point clouds. IEEE Geoscience and Remote Sensing Letters 9 (1), 38–42. Golub, G., Loan, C., 1996. Matrix Computations. John Hopkins University Press, Baltimore. Hellwich, O., Ebner, H., 2000. Geocoding SAR interferograms by least squares adjustment. ISPRS Journal of Photogrammetry and Remote Sensing 55 (4), 277– 288. Leberl, F.W., 1990. Radargrammetric Image Processing. Artech House, Boston. Lenzmann, L., Lenzmann, E., 2004. Strenge Auswertung des nichtlinearen Gauss– Helmert-Modells. Allgemeine Vermessungsnachrichten (2), 68–73.

Lord, R., 2000. Aspects of Stepped-frequency Processing for Low-frequency SAR Systems. Phd Thesis, Department of Electrical Engineering, University of Cape Town. Magnard, C., Brehm, T., Essen, H., Meier, E., 2012. High resolution MEMPHIS SAR data processing and applications. In: PIERS Proceedings, pp. 328–332. Maksymiuk, O., Stilla, U., 2012. A concept for building reconstruction from airborne multi-aspect SAR data. In: Proceedings of IEEE International Geoscience and Remote Sensing Symposium, pp. 7405–7408. Middelmann, W., Pepelka, V., Thoennessen, U., 2003. Registration of multiaspect InSAR images. In: Proceedings of SPIE, vol. 5095, pp. 98–109. Mikhail, E.M., 1976. Observations and Least Squares. IEP, New York. Neitzel, F., 2010. Generalization of total least-squares on example of unweighted and weighted 2D similarity transformation. Journal of Geodesy 84 (12), 751– 762. Neitzel, F., Petrovic, S., 2008. Total Least Squares (TLS) im Kontext der Ausgleichung nach kleinsten Quadraten am Beispiel der ausgleichenden Geraden. Zeitschrift fuer Geodaesie, Geoinformation und Landmanagement 133 (3), 141–148. Powell, M., 1970. Numerical Methods for Nonlinear Algebraic Equations. Gordon and Breach, pp. 115–161 (Ch. A Fortran subroutine for solving systems of nonlinear algebraic equations). Raggam, H., Gutjahr, K., 2000. InSAR block parameter adjustment. In: Proceedings of the 3rd European Conference on Synthetic Aperture Radar, pp. 493–496. Reigber, A., Alivizatos, E., Potsis, A., Moreira, A., 2006. Extended wavenumberdomain synthetic aperture radar focusing with integrated motion compensation. IEE Proceedings Radar, Sonar and Navigation 153 (3), 301–310. Schimpf, H., Essen, H., Boehmsdorff, S., Brehm, T., 2002. MEMPHIS – a fully polarimetric experimental radar. In: Proceedings of IEEE International Geoscience and Remote Sensing Symposium, pp. 1714–1716. Schmitt, M., Stilla, U., 2010. Utilization of airborne multi-aspect InSAR data for the generation of urban ortho-images. In: Proceedings of IEEE Geoscience and Remote Sensing Symposium, pp. 3937–3940. Schmitt, M., Stilla, U., 2011. Fusion of airborne multi-aspect InSAR data by simultaneous backward geocoding. In: Proceedings of Joint Urban Remote Sensing Event, pp. 53–56. Schmitt, M., Magnard, C., Brehm, T., Stilla, U., 2011. Towards airborne single pass decimeter resolution SAR interferometry over urban areas. In: Stilla, U., Rottensteiner, F., Mayer, H., Jutzi, B., Butenuth, M. (Eds.), Photogrammetric Image Analysis, Lecture Notes in Computer Science, vol. 6952. Springer, Berlin/ Heidelberg, pp. 197–208. Schreier, G. (Ed.), 1993. SAR Geocoding: Data and Systems. Wichmann. Schwaebisch, M., 1998. A fast and efficient technique for SAR interferogram geocoding. In: Proceedings of IEEE International Geoscience and Remote Sensing Symposium Proceedings, pp. 1100–1102. Soergel, U., Cadario, E., Thoennessen, U., 2004. Geocoding and fusion of airborne high-resolution multi-aspect SAR data. In: Proceedings of 5th European Conference on Synthetic Aperture Radar. Soergel, U., Michaelsen, E., Thiele, A., Cadario, E., Thoennessen, U., 2009. Stereo analysis of high-resolution SAR images for building height estimation in cases of orthogonal aspect directions. ISPRS Journal of Photogrammetry and Remote Sensing 64 (5), 490–500. Sohn, H.-G., Song, Y.-S., Kim, G.-H., 2005. Radargrammetry for DEM generation using minimal control points. In: Proceedings of IEEE International Geoscience and Remote Sensing Symposium, July 2005, vol. 2, pp. 1162–1164. Stilla, U., 2007. High resolution radar imaging of urban areas. In: Fritsch, D. (Ed.), Photogrammetric Week. Wichmann, pp. 149–158. Wegner, J.D., Soergel, U., 2008. Bridge height estimation from combined highresolution optical and SAR imagery. In: International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, vol. 37, pp. 1071–1076. Wu, J., Lin, D.-C., 2000. Radargrammetric parameter evaluation of an airborne SAR image. Photogrammetric Engineering and Remote Sensing 66 (1), 41–47. Yue, X.J., Huang, G.M., Zhang, Y., Zhao, Z., Pang, L., 2008. Multi-photo combined adjustment with airborne SAR images based on F, Leberl ortho-rectification model. In: International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, vol. 37, pp. 357–360. Zitova, B., Flusser, J., 2003. Image registration methods: a survey. Image and Vision Computing 21 (11), 977–1000.