ESA Earth Home Missions Data Products Resources Applications
    22-Jul-2014
EO Data Access
How to Apply
How to Access
FRINGE '96 Workshop: ERS SAR Interferometry, 30 September - 2 October 1996
Services
Site Map
Frequently asked questions
Glossary
Credits
Terms of use
Contact us
Search


 
 
 

FRINGE 96

The Developing of a Wide Area Interferometric Processor

Daniel Carrasco, Susanna Sanz, Ricardo Sousa, Antoni Broquetas

Universitat Politècnica de Catalunya. Signal Theory and Communications Department. Electromagnetic Engineering and Photonics Group (D3-EEF).Gran Capitán s.n. Campus Nord UPC - D3. 08034 Barcelona SPAIN.
E-mail: diaz@voltor.upc.es
http://gaig.upc.es

Abstract

ERS-1 SAR interferometry has been assessed over small areas. However, the development of an operational interferometric processor capable to process a quarter SLC scene running on a low-end workstation is a challenging task which requires new adaptive algorithms and blocking strategies in order to deal with a large amount of data. This paper summarizes the present status of an interferometric processor under development. The processor covers all stages from SLC images to DEM. Special effort is dedicated to make it as automatic as possible. Its main features are as follow: two different algorithms are combined in order to achieve a robust registration, the spectral shift filtering and the interferogram filtering are performed by means of a short time Fourier Transform which allows adaptive processing and two phase unwrapping algorithms are used: either LMS or weighted LMS. LMS performs well when dealing with not too difficult relief. W-LMS has a higher computational load but it allows to isolate low coherence or hard to unwrap areas using masks. Some results are shown using ERS-1 data.
Keywords: SAR Interferometry processor, STFT, LMS phase unwrapping, DEM.

Introduction

The development of an interferometric processor which can obtain a Digital Elevation Model from a pair of SAR images over a full (100x100 km) or quadrant-size (50x50 Km) SLC image is a challenging task. Besides the huge amount of data (every quarter scene occupies 150 MB) which will require image blocking processing, the algorithms to be used must be adaptive and robust in order to deal with image oddities.

At the D3-EEF Department such an interferometric processor is currently being developed and some promising results have been obtained. This article summarizes the present status of the processor development, along with the key algorithms and the latest results.

The processor is divided in the following blocks: image registration, "common band" filtering, interferogram generation, ellipsoidal Earth correction, interferogram fringe filtering, phase unwrapping and geocoding. When processing a quarter SLC scene, 1GB of disk space is required besides the space used by the starting SLCs. Every stage is done by a different program but all of them are controlled from a single shell script and use only one parameter file which makes it user friendly. From the interferogram generation up to the final geocoding, the data flow can be conducted through pipes rather than through intermediate disk files in order to increase the throughput. If intermediate results are required at these final stages, the information should be redirected to the hard drive with a significant time penalty. Dynamic blocking strategies are adopted depending on the computer maximum available RAM. Next sections will show the most remarkable facts about the algorithms used in the processor.

Interferogram generation

The first step is the image registration, which is performed in two stages: coarse and fine registration. Coarse registration is based in an amplitude cross correlation in a window at the center of the scene. Fine registration must compensate the pixel misregistration and its variation along the image. The fine registration is calculated at a regular grid along the scene (for example 3 by 3 points) or at some user selected points, by maximizing the image coherence. An amplitude cross correlation is also performed in order to detect aberrant offsets between the images caused when the test point is located over the sea or a high incoherence area. The offsets obtained between the two images at these points are fitted to a fine registration model (a first or second order surface) which will allow to calculate the pixel misregistration at any point of the scene. Finally, a variable index interpolation along the image is performed based on ideal interpolation.

After the image registration, the "common band" filtering is required (Gatelli' 94) . This will improve the coherence between the images. However, it has been found that for baselines around 150 m, this coherence improvement in the interferogram is not significant because the total amount of phase residues after the interferogram filtering scarcely depends on whether the "common band" filtering was performed or not. Anyway, the "common band" filtering will be used in order to improve the quality of the coherence image which will be used at the later stages of the processing.

Next step is the interferogram formation by multiplying one image by the complex conjugate of the other one. The order of the multiplication is determined in order to assure the correct sign of the phase.

The "flat Earth" correction uses orbital information from the CEOS headers. The orbit is interpolated for every range line in the scene. By combining the range distance equation (a sphere), with the Earth model (WGS 84 ellipsoid) and the Doppler equation (a plane) (Schreier'93) and solving the non linear system by iterative methods, a "flat Earth" interferogram can be simulated and substracted from the real one so that the remaining fringes are due only to the relief.

Interferogram filtering with the Short Time Fourier Transform

The fringes in the Earth corrected interferogram are too noisy and need some filtering before the phase unwrapping process. It is well known that the Fourier Transform of the complex interferogram will show (ideally) a peak at the frequency determined by its fringe frequency. For low fringe frequencies (nearly flat terrain) the interferogram two-dimensional spectrum will be low-pass, i.e. the energy will be concentrated in the center. On the other hand, for high fringe frequencies (steep slopes) the position of the peak of maximum energy in the Fourier transformed interferogram will move away from the center towards high frequency areas. It follows up that the interferogram must be band-pass filtered around the frequency determined by its fringe frequency in order to reject all the noise out of the information band.

However, since the fringe frequency changes along the interferogram depending on the scene relief so does the position of the maximum in the transformed domain. Therefore, an adaptive filtering is needed. Should a low-pass or a band-pass fixed filter be used instead, the high frequency fringes will destroyed. The Short Time Fourier Transform (STFT), widely used in image processing, is a suited technique for this purpose (Lim'80) .

STFT technique decomposes the interferogram in blocks which partially overlap in both directions. An overlapping factor of 50% was selected. Each block is windowed with a triangular (in both directions) window and Fourier transformed. The peak of the interferogram is detected and a two-dimensional band-pass filter is applied around it. This filter is designed with the window method using a Hanning Window. Once the block has been filtered, it is inverse-Fourier transformed back to the time domain. Since the blocks overlap and every sample in the interferogram belongs to four blocks (due to the 50% block overlap), the four contributions from each block must be added together to generate the filtered interferogram.

The dimensions of the blocks in which the interferogram is decomposed are selected so that the fringe frequency is uniform inside them. The aspect ratio of each block is 4:1 (azimuth-range) in order it to be square in the scene.

STFT not only allows adaptive fringe filtering but provides a straight forward blocking strategy for processing big interferograms. The residue reduction obtained reaches the 50% if compared with a non adaptive filtering processing (Carrasco'96). STFT is also used at the "common band" filtering stage.

The filtered interferogram will go to the phase unwrapping stage, but it can also be used to obtain a slope compensated coherence estimation. First of all, the round Earth correction must be undone in order to obtain the full interferometric phase (relief and flat Earth component). If the resultant phase is substracted from one of the original images (after registration and range filtering) the coherence calculation will not be byased by the interferometric phase and bigger patches can be used improving the estimation quality. This method outperforms the results obtained by only substracting the flat Earth component by means of a FFT spectral estimation since it also cancels the relief term.

Phase unwrapping

Two phase unwrapping algorithms have been tested: the LMS and the Weighted-LMS (Ghiglia'94). The LMS is a well known method as a counterpart to traditional branch-cut/ghost-lines or residue tying strategies (Goldstein'88). By minimizing the difference (in a least minimum squares sense) between the partial derivates of the wrapped and the unwrapped phase, which turns into the Poisson equation, it always provides a solution. However, since it considers the scene as a whole and it does not allow greater than pi phase jumps in the unwrapped phase, the solution may suffer from serious averaging problems around the phase discontinuities. This is because LMS provides only the irrotational component of the phase field. The rotational part, not present in the LMS solution, contains the residues and therefore the phase discontinuities greater than pi. LMS can be implemented in a very efficient way by using FFTs (Pritt'94) and performs well when dealing with not too rough relief.

The weak point of LMS, the inadequate treatment of discontinuities and areas where the phase is not reliable such as water surfaces, can be overcome by using the WLMS algorithm. WLMS adds a weight between 0 (not reliable) and 1 (reliable) to each sample of the phase map. This is equivalent to weighting each equation in the Poisson problem. The linear system can no longer be solved by using FFTs. The new weighted linear sparse system is solved iteratively with a preconditioned conjugated gradient algorithm. At every iteration step, the LMS algorithm is applied and the precondition used is the solution of the unweighted LMS problem. Convergence is guaranteed in N steps for NxN problem, but it is usually reached much faster, in a few iterations, depending on the steep of the discontinuities.

Several weighting procedures have been tested. First of all, the coherence map was used. It was found that the WLMS solution was better if the coherence map was quantized to just two levels, 0 and 1. This kind of mask allows the isolation of water covered surfaces and other areas where due to the lack of coherence phase unwrapping is not possible. The second masking strategy uses residues by surrounding them with zeros (morphological dilatation). This approach performs well when dealing with areas with short discontinuities (short ghost lines) which are masked by the dilated residues. It must be said that if an area is weighted with zeros, the WLMS algorithm will calculate the phases inside of it by interpolating from the values in the border. This will allow eventual phase discontinuities greater than pi along the zero-weighted area as though it was a ghost line. Should any knowledge of the ghost lines be available it could be used this way to allow phase discontinuities in the solution leaving all the noise residues to be solved by the LMS minimization.

In the figures 1, 2 and 3 an example of WLMS phase unwrapping is shown. figure 1 shows a 8x8 Km interferogram. Figure 2 shows the mask used from the dilatation of the residues. On the right side it can be seen the phase discontinuities caused by the river Ebro. The upper left side of the side shows residues due to de rough terrain. Finally, figure 3 shows the result from the WLMS.

Figure 1: 8x8 Km interferogram showing river Ebro in Tarragona, Spain.

Figure 2: Mask obtained from the residue morphological dilatation.

Figure 3: WLMS unwrapped phase

WLMS computational load and memory requirements are higher than LMS, in fact, as it was stated before, every WLMS iteration requires an LMS execution. Both algorithms will have to be modified in order not to keep all the image in memory but use a blocking strategy.

Some results

Figure 4 shows a 50x50 Km (quarter scene) filtered interferogram. It corresponds to an area crossed by river Ebro in Tarragona, Spain. Hilly areas can bee seen at both left and right sides whereas at the upper and the lower the relief is much smoother.

Figure 4: 50x50 Km filtered interferogram (thumbnail)

Next figures show a complete interferometric process from the SLCs up to the DEM. The perpendicular baseline is 160 m and the time between the two acquisitions was three days. The scene is a 16x16 km area in Tarragona (city) and Salou Cape - Port Aventura. The relief has no remarkable discontinuities, but all the lower side of the image corresponds to the Mediterranean sea. The phase unwrapping algorithm was the WLMS which allowed the isolation of the sea. The heights range from 0 to 175 m and the relief has been magnified in the final DEM. Geocoding was done using the orbit information to project the slant range image over the ground.

Figure 5: Amplitude detected image of the City of Tarragona (notice the harbour) and Salou Cape in Spain.

Figure 6: Coherence image

Figure 7: Earth Corrected interferogram

Figure 8: Filtered interferogram

Figure 9: DEM.

Conclusions

A SAR interferometry processor and its main algorithms along with some results has been presented. It is capable to process quadrant size images although phase unwrapping keeps being the critical step depending on the scene relief. Further work will be oriented to adopt a blocking strategy in the LMS/WLMS phase unwrapping in order it not to be limited by the RAM available and test new phase unwrapping algorithms capable to deal with rough relief such as the multibaseline approach. Geocoding is being improved by using the same geometric approach than in the "flat Earth" removal in order to generate a UTM positioned DEM.

Acknowledgments

This work is supported by INDRA ESPACIO and the Comisión Interministerial de Ciencia y Tecnología (CICYT, TIC 96-0879). The authors acknowledge the ESA ERS Fringe Working group for providing the ERS-1 images.

References

Carrasco, Sanz, Sousa, Broquetas, 1996:
Interferometría SAR con datos ERS-1. Actas del XI Symposium Nacional de la URSI, Madrid, septiembre 1996, vol. 2, pp.45-48.
Gatelli et al., 1994:
The wavenumber shift in SAR interferometry. IEEE Trans. on Geoscience and Remote Sensing, vol. 32, nº4, pp. 855-864, July, 1994.
Ghiglia, Romero, 1994:
Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods. Journal of the Optical Society of America - A, vol. 11, nº1, pp. 107-117, Jan. 1994
Goldstein, Zebker, Werner, 1988:
Satellite radar interferometry: two-dimensional phase unwrapping, Radio Science, vol. 23, nº4, pp. 713-720, Jul.-Aug. 1988.
Lim, 1980:
Image restoration by short-space spectral substraction. IEEE Trans. on Acoustics, Speech and Signal Processing, vol. 29, pp. 191-197, 1980.
Pritt, Jerome, Shipman, 1994:
Least-Squares two-dimensional phase unwrapping using FFT´s, IEEE Trans. on Geoscience and Remote Sensing vol. 32, nº3, pp. 706-708, May 1994
Schreier (ed.), 1993:
SAR Geocoding: data and systems. Wichmann, 1993.

Keywords: ESA European Space Agency - Agence spatiale europeenne, observation de la terre, earth observation, satellite remote sensing, teledetection, geophysique, altimetrie, radar, chimique atmospherique, geophysics, altimetry, radar, atmospheric chemistry