A Phase Unwrapping Method Based on Network Programming
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.
Let be a real valued function defined on a rectangular grid, and let
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:
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:
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:
subject to the constraints
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:
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.
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).
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.
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.
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.
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.
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.
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
|Copyright 2000 - European Space Agency. All rights reserved.|