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

A Phase Unwrapping Method Based on Network Programming

Mario Costantini ESA ESRIN
Via G. Galilei, c.p. 64
00044 Frascati, Roma, Italy
Phone/Fax: +39-6-94180516/280
E-mail: marioc@mail.esrin.esa.it

Abstract

Phase unwrapping is the reconstruction of a function on a grid given its values modulo . A new phase unwrapping method, which is a significant advance from existing techniques, is described and tested. The method starts from the fact that the neighboring pixel phase differences can be estimated with possibly an error which is an integer multiple of .This suggest the formulation of the phase unwrapping problem as a minimization problem with integer variables. In fact, it is possible to equate this to the problem of finding the minimum cost flow on a network, for the solution of which there exist very efficient techniques. The tests performed confirm the validity of this approach.
Keywords: Phase Unwrapping, Network Programming, SAR Interferometry

Introduction

Phase unwrapping is the reconstruction of a function on a grid given the value modulo of the function on the grid. We refer to these two quantities as the unwrapped and the wrapped phase respectively. In the last few years an increasing interest has been devoted to phase unwrapping, mainly due to the developing of SAR (Synthetic Aperture Radar) interferometry (Zebker and Godstein, 1986), although there are applications of phase unwrapping in several other fields.

Basically all the existing phase unwrapping techniques start from the fact that it is possible to estimate the neighboring pixel differences of the unwrapped phase when these differences are less than . From these, the unwrapped phase can be reconstructed up to an additive constant. The methods differ in the way they overcome the difficulty posed by the fact that this hypothesis may be somewhere false, which cause the estimated unwrapped phase differences to be inconsistent, that is their "integral" depends on the integration path.

Branch cuts methods (Goldstein et al.,1988; Prati et al., 1990) unwrap by integrating the estimated neighboring pixel differences of the unwrapped phase along paths avoiding the regions where these estimated differences are inconsistent. The problem of building cuts delimiting these regions is very difficult and the resulting phase unwrapping algorithm is very expensive computationally.

In least squares methods (Fried, 1977; Hudgin, 1977; Hunt, 1979), unwrapping is achieved by minimizing the mean square deviation between the estimated and the unknown neighboring pixel differences of the unwrapped phase. Least squares methods are very efficient computationally when they make use of fast Fourier transform techniques (Takajo and Takahashi, 1988; Ghiglia and Romero, 1994). But the resulting unwrapping is not very accurate, because least squares procedures tend to spread the errors that are instead concentrated on a limited set of points. To overcome this problem a weighting of the wrapped phase can be useful. However, the weighted least squares algorithms proposed (Ghiglia and Romero, 1994; Pritt, M. D., 1996) are iterative and not as efficient as the unweighted ones. Moreover, the accuracy of the results depends on the weighting mask used.

Here we propose a new method for phase unwrapping. We exploit the fact that the neighboring pixel differences of the unwrapped phase are estimated with possibly an error which is an integer multiple of . This leads to formulate the phase unwrapping model as the problem of minimizing the deviations between the estimated and the unknown neighboring pixel differences of the unwrapped phase with the constraint that the deviations must be integer multiple of . This constraint prevents the errors from spreading; so that a weighting of the data is not necessary. In any case, weighting is allowed in our formulation (without loss of efficiency), and can be useful when there are large regions of noisy data. Minimization problem with integer variables are usually very complex computationally. However, recognizing the network structure underlying our problem makes available very efficient strategies for its solution. In fact, our problem can be equated to the problem of finding the minimum cost flow on a network, for the solution of which there exist very efficient algorithms.

Finally, we present tests of our method performed on simulated data and on the real wrapped phase obtained from the interferometry of a pair of ERS-1 and ERS-2 Tandem SAR images. The results are really satisfactory: they demonstrate the consistency and efficiency of the method, which seems to be accurate. Comparative tests using a least squares method are presented too.

The Method

Let be a real valued function defined on a rectangular grid, and let

(1)

where, for a generic real , , with the integer such that . For simplicity of notation, in (1) as in the following, we adopt the convention that in any formula the indices run over the grid where the formula is defined. We refer to and as the unwrapped and the wrapped phase functions respectively. The inversion of (1), that is the reconstruction of from is the phase unwrapping process.

Let us define:

(2)

(3)

When and , we have that and . We assume that these relations hold in the majority of cases. In general we can restate the phase unwrapping problem of inverting (1) as the problem of finding the following residuals:

(4)

(5)

from which the neighboring pixel differences of the unwrapped phase can be calculated; then, through their integration, the unwrapped phase is reconstructed up to an additive constant which is an integer multiple of .

Let and be given non-negative real numbers weighting the a priori confidence that the residuals and must be small; when such an priori knowledge is not available, all the weights and are chosen equal to 1. We can estimate the residuals and as the solution of the following minimization problem:

(6)

subject to the constraints

(7)

 

(8)

(9)

The objective function to be minimized in (6) comes from the assumption that the residuals and are zero almost always; the absolute value function is chosen as the error criterion because it allows an efficient solution of the minimization problem, through the transformation shown below. The constraints in (7) express the property that and represent the neighboring pixel differences of a function (the unknown function ), as it results from (4) and (5); these constraints ensure that the unwrapped phase reconstructed from the solution of the minimization problem does not depend on the integration path. The constraints in (8) and (9) are a consequence of the fact that and take integer values, as it can be seen from (1), (2), (3), (4) and (5); as a consequence, the unwrapped phase achievable from the solution of the problem above is identical to the original wrapped phase when re-wrapped.

The problem given in (6), (7), (8) and (9) is a non-linear minimization problem with integer variables. Consider the following change of variables:

(10)

(11)

It can be seen that, through (10) and (11), the problem stated in (6), (7), (8) and (9) can be transformed so that it defines a minimum cost flow problem on a network, with the new variables representing the flow along the arcs of the network. In Figure 1 it is shown the network associated with the phase unwrapping problem.

Figure 1: The network associated with the phase unwrapping problem: the circles and the arrows represent the nodes and the arcs of the network respectively, while the boundary arcs are connected to the "earth" (by analogy with electrical networks) node

In the transformed problem the objective function in (6) becomes the total cost of the flow; the constraints corresponding to those given in (7) express the conservation of flow at the nodes; finally, the constraints in (8) and (9) are replaced by constraints defining the capacities of the arcs.

The transformation of the problem defined in (6), (7), (8) and (9) into a minimum cost flow problem on a network makes for an efficient solution, both as regards the memory and the computation required. An exhaustive review of algorithms for the solution of minimum cost network flow problems can be found for example in (Ahuja et al., 1993).

Experimental Results

The examples chosen to test our algorithm come from SAR interferometry, which provides one of the most difficult and interesting application of phase unwrapping. The first data set is a simulation of the wrapped phase obtained from SAR interferometric images relative to a scene consisting of a cone on a flat surface. The parameters of this experiment are chosen to simulate a real experiment with the ERS-1 SAR; for clarity, we have not introduced noise in the data. The simulated wrapped interferometric phase relative to the cone scene is shown in Figure 2. The neighboring pixel differences of the unwrapped phase estimated from the wrapped phase are somewhere inconsistent because the unwrapped phase is discontinuous (due to the simulation of the SAR phenomenon of layover). In Figure 3 it is shown where these inconsistencies are located. The unwrapped phase reconstructed with our algorithm (without weighting) is shown in Figure 4. It is perfectly identical to the "true" unwrapped phase resulting from the simulation. For a comparison, in Figure 5 we show the results of the unwrapping obtained using an unweighted least squares algorithm. Note the propagation of errors from the regions where the estimated neighboring pixel differences of the unwrapped phase are inconsistent.

Figure 2: The wrapped phase obtained from a simulated interferometric experiment (gray levels representation)

Figure 3: The regions where the neighboring pixel differences of the unwrapped phase estimated from the data of Figure 2 are inconsistent (in white)

Figure 4: The unwrapped phase reconstructed from the data of Figure 2 with our algorithm (perspective view); as the reconstruction is perfect, this is the "true" unwrapped phase obtained from the simulation

Figure 5: The unwrapped phase reconstructed from the data in Figure 2 with a least squares algorithm (perspective view)

The second data set is the wrapped phase obtained from the interferometry of two real images taken during the tandem campaign by the ERS-1 SAR and the ERS-2 SAR over the region of Etna volcano, Sicily, Italy. The wrapped phase is shown in Figure 6. Although the quality of the data is quite good, there are many inconsistencies in the neighboring pixel differences of the unwrapped phase estimated from the wrapped phase. The inconsistencies are due to phenomena well known in SAR interferometry: in particular can be recognized low coherence corresponding to the sea region at the right bottom of the image, layover in correspondence of the peak of the volcano, and maybe both of these phenomena corresponding to the mountainous region at the top of the figure. The location of the inconsistencies is shown in Figure 7. In Figure 8 we show the unwrapped phase reconstructed with our algorithm (without weighting). The result is visibly good: in particular, very sharp details can be appreciated. Finally, in Figure 9 we show the results of re-wrapping the unwrapped phase reconstructed using an unweighted least squares algorithm. By comparison with Figure 6, it can be noted that this wrapped phase is significantly different from the original one, for example in correspondence of the peak of the mountain. On the contrary, remember that the unwrapped phase obtained with the method proposed in this paper is identical to the original wrapped phase when re-wrapped.

Figure 6: The wrapped phase obtained from the interferometry of a pair of SAR images of the Etna volcano region, Sicily, Italy (gray levels representation)

Figure 7: The regions where the neighboring pixel differences of the unwrapped phase estimated from the data of Figure 6 are inconsistent (in white)

Figure 8: The unwrapped phase reconstructed from the data of Figure 6 with our algorithm (perspective view; the gray level is given by the ERS-1 SAR image intensity)

Figure 9: The wrapped phase obtained re-wrapping the phase unwrapped from the data in Figure 6 with a least squares algorithm (gray levels representation)

The computation efficiency of our phase unwrapping technique depends on the minimum cost flow algorithm used. In fact, many different strategies are possible to solve minimum cost problems on a network. The algorithm employed in the above examples requires about 10 seconds for the pixels simulated data and about 2 hours for the pixels Etna data on a Silicon Graphics Power Onyx RE2 (using 1 cpu R8000 at 75 Mhz). Probably, the efficiency of our method can be improved using a technique optimized for the particular network resulting from the phase unwrapping problem.

Conclusions

We propose a new method for phase unwrapping which appears to be accurate and efficient. The key points are: to formulate the phase unwrapping problem exploiting globally its integer qualities, which ensures accurate results; and to recognize the network structure underlying our formulation of the phase unwrapping problem, which makes for an efficient solution. The tests performed demonstrate the validity of this approach.

Details on the technique proposed will be reported in a more complete paper in preparation. Further validation tests are in progress and possible efficiency improvements are under investigation.

Acknowledgments

We would like to acknowledge Karl Heinz Theil and Xiaoqing Wu, of the University of Stuttgart (Germany) Institute of Navigation, for providing the wrapped phase of Figure 6, and Fabrizio Rossi, of the UniversitÓ dell'Aquila (Italy) Dipartimento of Matematica Pura e Applicata, for useful discussions on network programming algorithms. We also would like to thank Alfonso Farina, of the Alenia Area Sistemi, Roma, and Francesco Zirilli, of the UniversitÓ di Roma Dipartimento di Matematica, for their advice and encouragement. Finally, we wish to thank Steve Coulson and the staff of ESA ESRIN, Frascati (Italy), for support during the research reported on this paper.

References

Ahuja, R. K.,Magnanti, T. L., Orlin, J. B., 1993:
Network Flows: Theory, Algorithms, and Applications. Prentice-Hall, Englewood Cliffs, New Jersey, USA.
Fried, D. L., 1977:
Least-Squares Fitting a Wave-Front Distortion Estimate to an Array of Phase-Difference Measurements. Journal of the Optical Society of America, 67, pp.370-375.
Ghiglia, D. C., and Romero, L. A., 1994:
Robust Two-Dimensional Weighted and Unweighted Phase Unwrapping Using Fast Transforms and Iterative Methods. Journal of the Optical Society of America A, 11, pp. 107-117.
Goldstein, R. M., Zebker, H. A., and Werner, C. L., 1988:
Satellite Radar Interferometry: Two-Dimensional Phase Unwrapping. Radio Science, 23, pp. 713-720.
Hudgin, R. H., 1977:
Wave-Front Reconstruction for Compensated Imaging. Journal of the Optical Society of America, 67, pp. 375-378.
Hunt, B. R., 1979:
Matrix Formulation of the Reconstruction of Phase Values from Phase Differences. Journal of the Optical Society of America, 69., pp. 393-399.
Prati, C., Giani, M., Leuratti, N., 1990:
SAR Interferometry: a 2D Phase Unwrapping Technique Based on Phase and Absolute Value Informations. IGARSS '90, pp. 2043-2046.
Pritt, M. D., 1996:
Phase Unwrapping by Means of Multigrid Techniques for Interferometric SAR. IEEE Transactions on Geoscience and Remote Sensing, 34, pp. 728-738.
Takajo, H. and Takahashi, T, 1988:
Noniterative Method for Obtaining the Exact Solution for the Normal Equation in Least-Squares Phase Estimation from the Phase Difference. Journal of the Optical Society of America A, 5., pp. 1818-1827.
Zebker, H. A., Goldstein, R. M., 1986:
Topographic Mapping from Interferometric Synthetic Aperture Radar Observations. Journal of Geophysical Research, 91, pp. 4993-4999.

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