Alenia Spazio has developed, on internal funds, a complete set of processing tools to produce three-dimensional images of the observed area using two complex images coming from SAR sensors like ERS-1 and SIR-C/X-SAR. This paper describes the whole processing of the complex SAR data to obtain the 3-D representation of the image.

In the frame of interferometric processing some steps have been pointed out: the evaluation of range and azimuth shifts between the two images; the filtering of the complex data in order to improve the coherence between the images; the extraction and filtering of the interferogram; the phase unwrapping; the phase-to-height conversion.

The whole processing has been tested and applied to images coming from SIR-C/X-SAR (and ERS-1) on Monte Etna in Sicily, in order to evaluate the eventual differences between the two kinds of sensors in terms of radar wavelength, baseline, and repetition time of the images.

2. Shift evaluation and images registration

The registration of the two SLC images is done using the following procedures, after the removal of the coarse misregistration carried out by estimating the pixel shift with the cross-correlation on the amplitude images. For the fine (sub-pixel) registration we use the technique proposed in Lin, et al., 1991.

(1)

The fringe fluctuation function (1) has been improved with respect to the original algorithm by introducing pixel by pixel a weight factor derived from the correlation matrix between the two images. Fig. 1 shows the results obtained with the fringe fluctuation algorithm applied to X-SAR images coming from Monte Etna in Sicily, during the SRL-2 mission.

3. Complex images filtering

The objective is to improve the correlation coefficient and so to reduce the phase and height noise observing that one of the causes of phase noise is the relative shift of the reflectivity spectra of the two images. This is principally due to baseline decorrelation.

3.1 Baseline decorrelation

It is well known (Rodriguez and Martin, 1991) that the spectra of the received signals represent different terrain reflectivity bands. This is due to the different angle of view of each pixel on the ground with respect to the two points of observation. Calling q1 the view angle of a pixel for the first antenna (platform pass) and q2 the view angle of a pixel for the other antenna (platform pass) it has been noted that the same spectral components of the signal coming from one acquisition is shifted by a frequency fW with respect to the other, being:

(2)

with , f0 the radar frequency, and a the local mean terrain slope. Observing that , the spectral shift becomes

(3)

The behavior of frequency shift versus local slope is plotted in Fig. 2 for X-SAR images.

One can observe that f increases with increasing terrain local slope. Let us call "blind angles" those angles which cause an infinite shift. Moreover if the terrain slope exceeds the off-nadir angle (layover) the absolute value of f decreases and it becomes negative.

As seen this spectral shift depends on baseline and local slope and it causes that part of each image spectrum, exactly f width, that is not common in the two images. When, extracting the interferogram, the spectra correlation is performed, the result is noisy because of the correlation of these not-common parts: hence the cross-correlation presents the following three contributions: a peak at the frequency of perfect alignment of the spectra, a sequence due to a common band in adifferent position, and a sequence due to the cross-correlation between the not- common bands.

Fig. 2. Spectral shift vs local slope (X-SAR)

It is important to remember (Prati and Rocca, 1993) that this spectral shift is the cause of the fringe generation. The interferogram, hence is described by

(4)

with Ry range resolution. So we do not have to remove this shift, but only to filter the not common parts of the spectra.

To do this, because of the dependency from terrain slope, it is necessary to estimate the shift in a sufficiently small portion of image in which a constant slope is supposed. In the following the processing of one subimage is described.

3.2 Range Filtering

It is necessary to design two filters with the same bandwidth and different central frequency. These will be obtained by realizing a real high-pass filter H(f) with and shifting the filter once on and once on .

After having estimated the spectral shift the project of filters is performed: the high pass FIR filter is realized with the Parks McClellan technique.

To determine the optimum filter the Golden Section algorithm has been used: it seeks the maximum of a function to find, varying the stopband from zero to the stopband at -31 dB, the optimum stopband that maximizes the correlation coefficient . It is interesting to note that if the scene is supposed to be white and uniform, the correlation coefficient between the two images correspond to the correlation coefficient between the transmitted signal and the f shifted one.

So, the maximum detection algorithm can determine the optimum stopband that maximizes the correlation between the transmitted signal filtered with the centre on f₁ filter and the same signal filtered with the centre on f₂ and shifted by f filter.

3.3 Azimuth Filtering

The results obtained applying the processing described before to X-SAR and ERS-1 images lead to some observations.

The azimuth spectra of the X-SAR present two different Doppler centroids. In this case there is a common and a not-common part of the azimuth spectra and so in X-SAR images also an azimuth filtering is necessary.

In this case the Doppler centroids have been estimated and with the Parks McClellan method a high pass filter has been projected with N=127, , and , .

4. Phase unwrapping

4.1 Bidimensional analytical approach

The phase unwrapping is based on the well known weighted least-mean square estimation of the phase from the phase differences. It has been implemented on a general purpose computer for testing (maps of up to 4096x4096 pixels have been processed and unwrapped), and as a basis for a next implementation on a parallel computer.

Another very promising technique is based on direct integration of the Poisson equation coming from the derivation (second derivative) of the complex logarithm of the interferogram. With the second derivative operation we have reduced the problem to a well known one. Now, we have to solve the Neumann problem with non-homogeneous boundary conditions: the first derivative along the boundary has been previously calculated. Techniques and s/w programs to solve the above problem are widely described in the literature. We use the Fourier method to solve the homogeneous problem and then add the particular solution for the boundary conditions.

Two major practical points are: i) the calculation of derivatives; ii) the presence of noise (i.e. corrupted samples). In our implementation we use the Savitzky-Golay (S-G) smoothing filters for both computing the derivatives and smoothing noise. The characteristic of the S-G filter is to use c_n coefficients that preserve higher moments: for each point f_i, the least square fitted polynomial to all n_L+n_R+1 points is computed, and g_i is the value of that polynomial at position i..

The S-G filter allows us to compute directly the filtered derivative, simply taking the n^th derivative of the fitted polynomial and dividing g_i by ⁿ ( is the stepsize).

4.2 Global classification approach

Why do some algorithms find a lot of difficulties in doing what the human brain-eye system does in a very easy way? A try to understand the interpretation mechanism of human visual inspection of interferometric fringes leads to a new method for phase unwrapping.

The fringes, after their generation, are examined by a continuity tester which performs a segmentation of the 2-D fringe map into homogeneous areas. These bidimensional segments are then grouped into classes, the result of which is defined as connected sets of pixels having the same phase ambiguity.

Finally the discrete ambiguity problem is solved against classes which are mapped into a ambiguity matrix.

This algorithm has been implemented on a conventional HP Apollo RISC machine and results are 10 times faster than others if applied to L-band interferometric data (SIR-C). The only drawback is that non-contiguous regions (e.g. islands within the sea) cannot be uniquely assigned to an ambiguity map, but this is a common problem in the phase unwrapping matter.

5. Phase to height conversion

Most of the methods suggested in the literature need the knowledge of at least three points on the ground of which height, position in image, and absolute phase is known. In order to demonstrate the feasibility of a global interferometric mission to produce DEM of the whole Earth surface and especially of those areas about which the topographic data are unknown, a method that reduces the a priori information is developed. Fig.3 shows the typical geometry, where h(i,j) is the height of the pixel in the object, (i,j) its slant range distance, q(i,j) the off-nadir angle from which it is observed, H the platform height and Re the Earth radius at the observed point. The idea is to refer the analysis to the interferometer axis: to this the off-nadir angle is equal to the tilt angle, that is the baseline attitude with respect to cross-track direction. Along this axis the interferometric phase should be zero: the interferogram phases are referred to this zero point in terms of variation with respect to different off-nadir angles. So, the interferometric absolute phase can be written as

. (5)

with .

Deriving (5) with respect to off-nadir angle

	(6)
	(7)
	(8)

Deriving (7) with respect to , , and h with H'=H+Re and substituting (6) we obtain every image pixel (i,j) slope behavior toward the interferometric phase and slant range increasing distance (as function of off-nadir angle), where (j is the range coordinate and f_r is the range sampling rate) and (i_m,j_m) is the only ground reference point about which we should know both its position in the image toward the interferometer system
( (i_m,j_m) , (i_m,j_m) ) and the absolute height (h (i_m,j_m) ). In this way it is possible to reconstruct the height in the whole image knowing information of just one ground control point and the unwrapped phases, with the simple implementation of a linear equation in which the constant term is given by the attitude of the interferometric system and by the reference point height. Moreover implementing the same procedure to the wrapped phase one can obtain the phase behavior without a flat terrain contribution, so that it is possible to know the level contour lines on the image. Fig. 4 shows the level contour lines referred to in Fig. 1. Fig. 5 shows the corresponding 3-D reconstruction.

6. Multifrequency (and Multibaseline) Interferometry.

The multifrequency capability of the SIR-C/X-SAR system has been used to improve the achievable performances of the final interferometric product. Two techniques have been used: the first one to help the processing of the low correlated X-band data; the second to improve the final sensitivity by the fusion of the L- and X-band DEM's.

The unwrapping of the X-band interferogram is difficult due to: i) the low level of correlation affecting large areas of the examined scene; ii) the quite large baseline which results in very narrow fringes where the slope is steeper. To help the unwrapping process the L-band interferogram is first unwrapped, then it is scaled (according to the wavelength ratio), registered, and subtracted from the X-band fringes. The result (delta fringes) is an interferogran at X-band with a few wide fringes which are easily unwrapped; the result is then added to the previously subtracted interferogram giving the final X-band product. The processing can be improved by weighting the "delta fringes" with the correlation map before unwrapping; obviously the result is a merging of the L-band (where the correlation is low) and X-band unwrapped phases.

The second technique (only partially tested) uses the multiresolution analysis (via the wavelet decomposition) to perform the fusion between the L- and X-band. The X-band and the scaled L-band unwrapped phases are firstly decomposed using the wavelet decomposition, then a final image is composed using the HF components of the X-band and the LF components of the L-band.

7. Conclusions

All the processing has been developed in order to realize a massive parallel technique for a following implementation on a parallel computer. The results are very interesting especially as far as the coherence improvement and phase unwrapping are concerned. The trend is that of a greater use of the interferometric technique to classify radar images using low noise coherence images and to use the low noise fringes to develop a differential Interferometry to detect small changes in the earth's surface.

References

[1] Bombaci, O., F. Rubertone, and A. Torre "Alenia Spazio Activities on SAR Interferometry: application on scientific mission," IGARSS `95, Florence, 1995.

[2] Bombaci, O. and A. Torre, "Alenia Spazio Research Activities on SAR Interferometry," EUROPTO `95, Paris, Sept. 1995.

[3] Bombaci, O., F. Impagnatiello, and A. Torre, "New features for on-line processing in interferometric applications" EUSAR `96, March 1996. [4] Lin, Q., J. F.Vesecky, and H. A. Zebker, "New Approaches in Interferometric SAR Data Processing," IEEE Trans. on Geoscience and Renote Sensing, vol. 30, no. 3, May 1992.

[5] Rodriguez, E. and J. Martin, "Theory and design of interferometric SAR's," Proceedings of IEEE, 1991.

[6] Prati, C. and F. Rocca,"Improving slant-range resolution with multiple SAR surveys," IEEE Transaction on Aerospace and Electronic Systems, vol. 29, no. 1, January 1993.

Fig. 1. X-SAR fringes over Monte Etna, Sicily, obtained with fringe fluctuation algorithm

Fig. 3. Phase-to-Height projection geometry

Fig. 4. Mt. Etna Flat Terrain subtracted interferogram (SIR-C L-band)

Fig. 5. Mt. Etna DEM

Table of Contents

Converted to HTML by Alvin Wong, al.wong@jpl.nasa.gov

The Jet Propulsion Laboratory
4800 Oak Grove Drive
Pasadena, Cailfornia 91109