| 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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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