Validation of a Novel Phase Unwrap Algorithm Using True and Simulated ERS Tandem SAR Interferometric Data

Abstract

Phase unwrapping is a key problem in all quantitative applications of SAR (Synthetic Aperture Radar) interferometry. We have developed a phase unwrapping method that exploits globally the integer qualities of the problem. Recognizing the network structure underlying our formulation of the phase unwrapping problem makes for an efficient solution. We present the results of some tests performed on true and simulated ERS tandem interferometric SAR data, which confirm the validity of our approach.
Keywords: phase unwrapping, network programming, SAR interferometry.

Introduction

Phase unwrapping is a key problem in all quantitative applications of SAR (Synthetic Aperture Radar) interferometry (Zebker and Godstein, 1986), and in several other fields (Oppenheim and Lim, 1981). 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.

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.

Recently has been pointed out how more sophisticated estimations of the neighboring pixel differences of the unwrapped phase can reduce errors (Davidson and Bamler, 1996). On the other hand, a different unwrapping approach has been proposed for the case when several related interferograms are available (Ferretti et al., 1996).

We propose (Costantini, 1996) 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 formulating the phase unwrapping model as the problem of minimizing the weighted 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 . With this constraint, the unwrapping results do not depend critically on the weighting mask used, and errors are prevented to spread.

Minimization problems with integer variables are usually very complex computationally. However, recognizing the network structure underlying the phase unwrapping problem makes it possible to employ very efficient strategies for its solution. In fact, it 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.

We present tests of our method performed on simulated and true ERS tandem SAR interferometric data. The results are very encouraging and demonstrate the robustness, accuracy and efficiency of the method.

The Phase Unwrap Method

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

,

(1)

where is an 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 pixel 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 preliminary estimates of the neighboring pixel differences of the unwrapped phase according to:

,

(2)

,

(3)

where and are integers selected a priori so that and hold in the majority of cases. In the experiments described here, and have been chosen to have and .

In general and cannot be consistently interpreted as neighboring pixel differences of the unwrapped phase and "integrated" to recover the unwrapped phase because their integral depends on the integration path. It is convenient to restate the phase unwrapping problem of inverting (1) as the problem of finding the following neighboring pixel difference residuals:

,

(4)

,

(5)

from which the neighboring pixel differences of the unwrapped phase are 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. In the examples shown in this paper we use weights and based on the consistency of and as neighboring pixel differences of a function (that is on if their "integral" is independent of the integration path), and on the coherence of the interferometric phase, which is a measure of reliability available for data coming from interferometric experiments. 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 to be estimated 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 are 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)

,

(12)

.

(13)

It can be seen that, through (10), (11), (12) and (13), the problem stated in (6), (7), (8) and (9) can be transformed so that it defines a minimum cost flow problem on a network (figure 1), with the new variables representing the flow along the arcs of the network.

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

Finally, it is important to note that, at the cost of finding a suboptimal solution instead of the optimal one, a block subdivision of the minimization problem can be implemented elegantly; that is enforcing to each block to be processed further constraints depending on the solution found in the blocks already processed. A block subdivision allows to unwrap phase data virtually size unlimited with a fixed memory requirement and a computational time linearly increasing with the size. On the other hand, the suboptimal solution found is practically the same quality as the optimal one when the block size is chosen sufficiently large. For example, the unwrap results shown in the following section have been obtained by splitting the pixel data in four blocks of pixels each.

With the implementation chosen, unwrapping the pixel data has taken about two hours on a Silicon Graphics Power Onyx RE2 (using one cpu R8000 at 75 Mhz). Note, however, that many different strategies can be used to solve minimum cost flow problems on networks, and the efficiency of our method can be improved using a technique optimized for the particular network resulting from the phase unwrapping problem.

The Validation Results

SAR interferometry provides one of the most difficult and interesting application of phase unwrapping. We have tested our algorithm on simulated and true SAR interferometric phase: the simulated data allow a quantitative validation of our unwrapping method, while the true data are useful to verify the robustness of the algorithm in a less controlled situation.

Figure 2. (a) The wrapped phase simulating the interferogram phase obtained from the 5 and 6 September 1995 ERS-1 and ERS-2 SAR images of the Etna volcano in Sicily, Italy (the image size is pixels; the gray scale represents the interval ). (b) The true data corresponding to those simulated in (a).

Both the simulated and true interferometric phases to be unwrapped (figure 2) have been generated by means of the interferometric processor DIAPASON (CNES, 1996), starting from the SAR images taken on September 5 and 6 of 1995 by the ERS-1 and ERS-2 satellites over the Etna volcano in Sicily (Italy), and using a digital elevation model of the same region produced from SPOT images. The simulated phase has been made more realistic by introducing noise based on the coherence of the corresponding true interferometric phase, which is a measure of the correlation of the two SAR images (Lee et al., 1994).

Figure 3. (a) The inconsistencies in the preliminary estimates (2) and (3) of the neighboring pixel differences of the unwrapped phase of figure 2(a) (the white points correspond to the non-zero "integrals" along elementary closed paths of those preliminary estimates; the dashed rectangle identifies the area enlarged in figure 4). (b) The same as (a) but referring to figure 2(b).

Although the quality of the data is quite good, there are a lot of inconsistencies in the preliminary estimates (2) and (3) of the neighboring pixel differences of the unwrapped phase (figure 3). The failures of those preliminary estimates can be related to phenomena well known in SAR interferometry: in particular can be seen low coherence regions at the left corner of the image, corresponding to the sea, and at the right side, corresponding to mountains, while layover structures are recognizable at the right side and at the center, corresponding to peaks of mountains and the volcano.

Figure 4. (a) The dashed area of figure 3(a) enlarged. (b) The dashed area of figure 3(b) enlarged. (c) The neighboring pixel difference residuals (4) and (5) found for the data of figure 2(a) by solving the minimization problem (6), (7), (8) and (9) (the white points correspond to the non-zero difference residuals). (d) The same as (c) but referring to figure 2(b).

The inconsistencies in the preliminary estimates (2) and (3) are overcome by solving the minimization problem (6), (7), (8) and (9) to determine where these estimates must be wrong; that is where the difference residuals (4) and (5) are found to be different from zero because the neighboring pixel differences of the unwrapped phase must be greater than p in absolute value (figures 4).

Figure 5. (a) The phase of figure 2(a) unwrapped by "integrating" its neighboring pixel differences determined by substituting in (4) and (5) the difference residuals which solve the minimization problem (6), (7), (8) and (9) (the color scale represent the interval ). (b) The same as (a) but referring to figure 2(b).

Finally (figures 5), the unwrapped phase is reconstructed by "integrating" its neighboring pixel differences obtained by substituting in (4) and (5) the difference residuals which solve the minimization problem (6), (7), (8) and (9). Remember that the reconstructed unwrapped phase is identical to the original wrapped phase when re-wrapped.

Figure 6. The difference between the reconstructed unwrapped phase of figure 5(a) and the simulated unwrapped phase (the gray scale represent the interval ).

To quantify the accuracy of the reconstructed unwrapped phase we have calculated the difference between this and the known simulated one. When unwrapping the simulated phase, the error committed is concentrated in layover and very low coherence areas, where interferometry does not furnish actual information, while is zero for the 99% of pixels (figure 6).

Figure 7. The difference between the reconstructed unwrapped phase of figure 5(b) and the simulated unwrapped phase (the gray scale represent the interval ; the contour lines are spaced of ).

Figure 8: The wrapped difference between the true wrapped phase of figure 2(a) and the simulated one shown in figure 2(b) (the gray scale represent the interval ).

In the case of the true data, the difference between the unwrapped phase and the simulated unwrapped phase is greater than half a cycle for a large portion of the data (figure 7). However, this difference is not due to phase unwrapping errors, but to the discrepancy (due to DEM errors and atmospheric artifacts) between the simulated and true wrapped phases (figure 8): in fact, it is evident (compare figures 7 and 8) that they matches very well. We recall that wrapping the difference between the unwrapped true SAR interferometric phase and the simulated unwrapped phase gives by construction the wrapped difference between the simulated and true wrapped phases.

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.

Acknowledgements

We would like to acknowledge Didier Massonet, of the CNES, Toulouse (France), for providing us with the Etna DEM used for the simulation of figure 2, and Fabrizio Rossi, of the University of L’Aquila (Italy) pure and applied mathematics department, for useful discussions on network programming algorithms. We also would like to thank Alfonso Farina, of the Alenia systems area, Rome, and Francesco Zirilli, of the University of Rome mathematics department, for their advice. Finally, we wish to thank Betlem Rosich Tell, Steve Coulson, and the staff of ESA ESRIN, Frascati (Italy), for support during the research reported here.

References

Ahuja, R. K., Magnanti, T. L., Orlin, J. B., 1993:
Network Flows: Theory, Algorithms, and Applications, Prentice-Hall, Englewood Cliffs, New Jersey, USA.

CNES, 1996:
Logiciel de traitement interferometrique automatise DIAPASON, CNES/833/2/95/0145, Toulouse, France.

Costantini, M., 1996:
A Phase Unwrapping Method Based on Network Programming, Proceedings of the "Fringe ‘96" Workshop, Zurich, Switzerland, 1996, ESA SP-406 (http://www.geo.unizh.ch/rsl/fringe96/papers/costantini/).

Davidson, G. W., and Bamler, R., 1996:
Robust 2-D Phase Unwrapping Based on Multiresolution, SPIE Proceedings Series, 2958, pp. 226-237.

Ferretti, A., Monti Guarnieri, A., Prati, C., and Rocca, F., 1996:
Multi Baseline Interferometric Techniques and Applications, Proceedings of the "Fringe ‘96" Workshop, Zurich, Switzerland, 1996, ESA SP-406 (http://www.geo.unizh.ch/rsl/fringe96/papers/ferretti-et-al/).

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.

Lee, J. S., Hoppel, K. W., Mango, S. A., and Miller, A. R., 1994:
Intensity and Phase Statistics of Multilook Polarimetric and Interferometric SAR Imagery, IEEE Transactions on Geoscience and Remote Sensing, 32, pp. 1017-1028.

Oppenheim, A. V., and Lim, J. S., 1981:
The Importance of Phase in Signals, Proceedings of the IEEE, 69, pp. 529-541.

Prati, C., Giani, M., Leuratti, N., 1990:
SAR Interferometry: a 2D Phase Unwrapping Technique Based on Phase and Absolute Value Informations, Proceedings of the 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