The proposed method for phase unwrapping is based on a global analysis of the interferometrical phase. The underlying principle is that the interferogram is partitioned such that the unwrapped-phase function on each element can be locally modelled by the mean values of the phase difference between neighbouring pixels in azimuth and range directions. Using this local information and a 

least-squares algorithm (Gauss-Seidel relaxation), an approximate model of the unwrapped phase is then generated and tested by calculating the “residue image” defined as the difference between the original interferogram and the model itself. If there are residual fringes, then the result must be iteratively refined applying the method to the residue image. The accuracy of the proposed estimation depends on the dimensions of the elements and the dynamic content of the phase, i.e. on the “roughness” of the ground surface

Bo, G.

Dellepiane, S.

English

Scientific Open-access Literature Archive and Repository

IL NUOVO CIMENTO Vol. 24 C, N. 1 Gennaio-Febbraio 2001Weighted multi-resolution phase-unwrapping method (∗)G. Bo and S. DellepianeDipartimento di Ingegneria Biofisica ed Elettronica - Universita` degli Studi di Genova, Italy(ricevuto il 25 Novembre 1999; revisionato il 2 Gennaio 2001; approvato il 23 Gennaio 2001)Summary. — The proposed method for phase unwrapping is based on a globalanalysis of the interferometrical phase. The underlying principle is that the inter-ferogram is partitioned such that the unwrapped-phase function on each elementcan be locally modelled by the mean values of the phase difference between neigh-bouring pixels in azimuth and range directions. Using this local information anda least-squares algorithm (Gauss-Seidel relaxation), an approximate model of theunwrapped phase is then generated and tested by calculating the “residue image”defined as the difference between the original interferogram and the model itself. Ifthere are residual fringes, then the result must be iteratively refined applying themethod to the residue image. The accuracy of the proposed estimation depends onthe dimensions of the elements and the dynamic content of the phase, i.e. on the“roughness” of the ground surface.PACS 93.85.+ – Instrumentation and techniques for geophysical research.PACS 93.65.+ – Data acquisition and storage.PACS 91.10.Jf – Topography; geometric observations.1. – IntroductionIn the last years several efforts have been made to solve the problem of phase un-wrapping in SAR interferometry [1-3] . The case of recovering a phase function fromhigh-quality data can be regarded as worked out, but when noise or aliasing corrupts theinterferogram many difficulties arise.In this article, an algorithm is proposed that is more robust to noise than classicalapproaches. This algorithm is based on the method described in [4] which generatesa global model for the unwrapped phase, while limiting the propagation of errors dueto the presence of noisy areas in the interferometrical image. The key idea is that theinterferogram can be subdivided into polygonal elements, on each of which the unwrappedphase is modelled by the local mean slopes. Starting from this information, it is possible(∗) Paper presented at the Workshop on Synthetic Aperture Radar (SAR), Florence, 25-26February, 1998.c© Societa` Italiana di Fisica 123124 G. BO and S. DELLEPIANEto locally approximate the true phase with linear functions, which in two dimensionsare planes, and then a global model is reconstructed by linking the nodes of the planesfound.Suppose a 2-D continuous space and let φ(x, y) and φp(x, y) represent, respectively,the unwrapped and the interferometrical phase functions, where x and y indicate theazimuth and range directions. Since the interferogram contains the principal values ofthe phase, i.e. modulo 2π, we may writeφ(x, y) = φp(x, y) + 2πm1(x, y) ,(1)where m1(x, y) is an integer-valued function of azimuth and range.If we now consider a function µ(x, y) approximating the true phase, it is possible tocalculate the “residue” associated to this model as the principal value of the differencebetween the interferogram and the model itself:R(x, y) = [φp(x, y)− µ(x, y)]p(2)Combining eqs. (1) and (2) yields[φp(x, y)− µ(x, y)]p = [φ(x, y)− 2πm1(x, y)− µ(x, y)]p ,(3)leading to[φp(x, y)− µ(x, y)]p = [φ(x, y)− µ(x, y)− 2πm1(x, y)] + 2πm2(x, y) ,(4)and thenφ(x, y) = [φp(x, y)− µ(x, y)]p + µ(x, y) + 2π(m1(x, y)−m2(x, y)) .(5)Assuming the model and the residue to be continuous, then 2π(m1(x, y)−m2(x, y))is a continuous function of integers and, therefore, must be a constant:φ(x, y) = µ(x, y) + [φp(x, y)− µ(x, y)]p + 2πC .(6)This means that, if one is able to generate a model with a continuous residue, thenthe unwrapped phase is given by the sum of the model with the residue image and aconstant. This constant can be calculated if the elevation of some point in the scene isknown.In the real case the interferogram is not continuous, but discrete. This implies thepresence of discontinuities between near points both in the model and in the residue,but if those discontinuities are limited and have a modulo value smaller than π, it is stillpossible to unwrap the phase.WEIGHTED MULTI-RESOLUTION PHASE-UNWRAPPING METHOD 1252. – The proposed approachIn order to determine a surface approximating the true phase function, the best so-lution would be the use of a finite-element based modelisation technique, dividing thedomain with triangles of various shapes and dimensions that best fit the local charac-teristic of the interferogram. Instead of this choice, in our approach we preferred to usesquared elements, mainly in the aim of improve the computational efficiency. The coarse-ness of the grid used for the partition of the interferogram is determined as a compromisebetween computational burden and precision of the result. Big squares result in a ratherpoor accuracy, in contrast with small squares. Once the size of the squares has beenfixed, the phase model generation goes through the following steps:– estimation of the mean slopes in azimuth and range for every square;– estimation of the phase values at the nodes of the grid through the solution of aleast mean squares problem;– computation of the model in each point using bilinear interpolation;– computation of the residue image.If the obtained solution is not accurate enough, an improvement is possible iterativelyapplying the algorithm to the residue image with a partition of different size.For the evaluation of the accuracy of the model, a criterion based on the number offringes in the residue image can be used. The interferometric fringes in the interferogramare those points that are characterised by a modulo of the gradient greater than π andhave therefore a discontinuity in phase. A point that belongs to a fringe therefore satisfiesthe following property [4]:[([φp(i+ 1, j)− φp(i, j)]p)2 + ([φp(i, j + 1)− φp(i, j)]p)2]1/2  π .(7)The objective of the process of phase unwrapping now becomes to find the model thateliminates the fringes in the residue image, i.e. that eliminates the points satisfying (7).Since in real world image one has to deal with noise and other phenomena that makesthis impossible, this constraint is relaxed to minimising the number of points that verify(7).2.1. The computation of the phase model . – With the approach proposed in [4] themean slopes, in azimuth and range, are determined adjusting a bilinear model (i.e. aplane) to each elementary square. The algorithm searches for the plane that best fitsthe behaviour of the phase in every square of the grid. To this end, a general plane isconstructed:µloc(i, j) = slopeazimuth × (i− i0) + sloperange × (j − j0) , i, j ∈ ES ,(8)where i0 and j0 indicate the location of an arbitrary initial point and ES indicates theelementary square. The parameters slopeazimuth and sloperange should be chosen suchthat they minimise the number of fringes in the residue which is associated to the square:R(i, j) = [φp(i, j)− slopeazimuth × (i− i0)− sloperange × (j − j0)]p , i, j ∈ ES .(9)126 G. BO and S. DELLEPIANEA successive minimisation of the variance on the elementary residue allows refiningthe obtained slopes. This technique for estimating the local information, though robustto noise, results in a noticeable computational burden and causes the process of phaseunwrapping to be slow. In order to solve this problem, we propose an alternative methodfor computing local slopes that, though faster and simpler, maintains the ability oflimiting the propagation of errors due to the presence of noise.For the purpose of understanding well the proposed modification, we consider for amoment the 1-D case and assume that the phase is represented by a continuous function.If the derivative of the phase function is known, the phase of a generic location can becomputed by integration:φ(x) = φ(x0) + ∫x0φ′(x)dx .(10)The derivative, φ′(x), describes the local variations and therefore the local slopes ofthe phase function. If we assume that the behaviour of the phase is characterised by acomponent with a constant slope to which a dynamic component is added that accountsfor small variations, we can rewrite the first derivative asφ′(x) = φ′cost + φ′dynam(x) ,(11)leading toφ(x) = φ(x0) + φ′cost × (x− x0) + ∫x0φ′dynam(x)dx .(12)Considering the discrete case, if we assign the mean slope of the phase function tothe constant component of the first derivative, we obtainφ(xn) = φ(x0) + meanslope× (n∆x) + w(x0, xn) ,(13)where w(x0, xn) represent the dynamic content of the phase and ∆x is the spatial reso-lution.Consider now the discrete form of (10):φ(xn) = φ(x0) +∑n−1∆xi ,(14)where ∆xi is the principal value of the phase difference between neighbouring pixels:∆xi =W [φ(xi+1)− φ(xi)] , i = 1, . . . , n− 1 ,(15)W being the wrapping operator.The mean value of the local phase differences for all the points is given byE{∆x} = (∑n−1∆xi)/(n− 1) = (∑n−1W [φ(xi+1)− φ(xi)])/(n− 1) .(16)WEIGHTED MULTI-RESOLUTION PHASE-UNWRAPPING METHOD 127At this point, adding and subtracting E{∆x} to and from each term in the summation(14) yieldsφ(xn) = φ(x0) +∑n−1(∆xi − E{∆x}+ E{∆x})(17)andφ(xn) = φ(x0) + n× E{∆x}+∑n−1(∆xi − E{∆x}) ,(18)where the spatial resolution is supposed to be unitary.Comparing (18) and (13) leads to the observation that the mean of the local phasedifferences can be used as an estimation of the mean slope. If in (18) the dynamic termis neglected, we may writeφ(xn) ∼= φ(x0) + n× E{∆x} .(19)Obviously, the accuracy of this approximation depends both on the dynamic contentof the phase and the length of the analysed sequence.In the 2-D case, we can now estimate the mean slope both in azimuth and in range,using the previous observations:meanslopeazimuth = E{∆x} =∑i∑j∆xij/N(M − 1) =(20a)=∑i∑jW [φi+1j − φij ]/N(M − 1) , i=1, . . . ,M−1, j=1, . . . , N ,meansloperange = E{∆y} =∑i∑j∆yij/M(N − 1) =(20b)=∑i∑jW [φij+1 − φij ]/M(N − 1) , i=1, . . . ,M, j=1, . . . , N−1 ,where M and N are the dimensions of the rectangle on which the phase function is tobe estimated. If the amount of points affected by noise or aliasing is not great comparedto the total number of pixels in the image, big squares are less sensitive to noise. On thecontrary, the use of small elements seems to be a better solution for capturing the detailsof the phase function, but results in more important effects of errors in the evaluationof local slopes. A compromise between the two possibilities could be reached using aniterative algorithm that involves a partition of different size for each iteration.Once the mean slopes in azimuth and range are known, the phase at the vertices ofthe model can be computed. This is done by minimising, for every set of vertices, thedifference between the phase variation on the model and the variation imposed by theslopes computed in the previous point. The solution is obtained by solving a least-mean-squares problem. From the physical point of view, this can be translated in the necessityof obtaining a continuous topographic model. Starting from the values of the phase inthe vertices, one can estimate the values of the model in all the points of the domainusing bi-linear interpolation.128 G. BO and S. DELLEPIANEFig. 1. – Result showing the capability of the method in limiting error propagation.a) Rewrapped model obtained by applying a classic least-mean-squares algorithm to a syntheticfringe pattern with a noisy area. b) The same as in a) but exploiting the good properties of theproposed method. c), d) The correspondent residual images.2.2. Implementation. – In summary, the following are the various steps necessary toperform phase unwrapping with the discussed approach. Let DIM be the size (in pixels)of the squares during a certain iteration:1) DIM = 64;2) divide the interferogram into squares of DIM×DIM pixel in size;3) compute, for each square, the mean of the difference in phase between neighbouringpixels, for both the azimuth and range direction (this determines the mean slope);4) compute the value of the phase in the vertices of the grid solving a least meansquares problem;5) compute the phase for all the points of the model using bilinear interpolation;WEIGHTED MULTI-RESOLUTION PHASE-UNWRAPPING METHOD 129Fig. 2. – Result showing the behaviour of the discussed approach in the presence of aliasingin the interferogram. A “pyramidal” phase function with an horizontal cut was used. a) Theoriginal synthetic interferogram. b) The rewrapped phase model. Also in this case the algorithmis able to limit the influence of defects in the interferometrical phase.6) compute the residue image;7) DIM = DIM/2;8) unwrapping of the residue image repeating the steps 2 through 7 until DIM = 2;9) successive unwrapping with squares growing in size until DIM=64.Let us notice the iterative character of the algorithm. At the end of each iteration thealgorithm generates an updated model and a residue image: if fringes are still presentin the residue image, this will be considered as a new interferogram to be unwrapped bythe same algorithm. The reconstructed phase function results from the sum of all thepartial models obtained with all iterations.In order to evaluate the accuracy of the model obtained with the steps describedabove, visual inspection of the residue image could be sufficient because the presence ofresidue fringes implicates the necessity of further processing.3. – Experimental resultsIn order to illustrate the usefulness of the method, both synthetic and real data setshave been used to test the method. Figure 1 shows the results obtained by applying botha least-square algorithm and the proposed method to a synthetic fringe pattern with anoisy area simulating the presence of low coherence in real data. In order to evaluatethese results, we compared the original interferogram with the rewrapped phase models.The distortions clearly visible in the image obtained using the least-square approach arestrongly reduced by the described algorithm. This means that the influence of noise isstrictly limited to the neighbourhood of noisy points.The fringe patterns in fig. 2 represents what has been obtained by using an interfer-ogram derived from a synthetic “pyramidal” phase function with an horizontal cut on a130 G. BO and S. DELLEPIANEFig. 3. – Sequence of images illustrating the process of development of the phase model.Rewrapped phase model (left), residue image (center) and a 3-D reconstruction of the un-wrapped phase (right), respectively, after the first (a), the third (b) and the fifth (c) iterationof the described approach. In d) the original interferogram is shown. It is evident how eachiteration introduces more details.WEIGHTED MULTI-RESOLUTION PHASE-UNWRAPPING METHOD 131face, in order to simulate an aliasing case. Also in this experiment the algorithm showsa good capability in preserving correct phase values from being corrupted by local errorspropagation.The objective of fig. 3 is to explain how the algorithm comes to the final result.In this case a real dataset has been used. A grey-levels representation of the originalwrapped phase is shown in fig. 3d. It has been obtained by applying the interferometricalprocess to a pair of complex images acquired by the ERS-1 and ERS-2 satellites, duringthe tandem campaign, over the Etna volcano in Sicily (Italy). The quality of the data isoverall rather high: in particular the concentric fringes corresponding to the mountain areclearly visible and appear almost free of noise. Notwithstanding these good properties,several regions in the interferogram are characterised by low coherence and inconsistenciesin the phase values. For instance, in the upper part of the image fringes are damaged bynoise and, close to the Etna’s peak, a layover area can be observed. Moreover, in the rightbottom portion of the interferogram a low-coherence area is visible, which correspond tothe sea. In such a situation it is fundamental to unwrap phases by means of an algorithmable to deal with the problem of error propagation. The obtained results confirm thatour approach has the capability of limiting the reduction of reliability in the final phasemodel, due to the presence of inconsistencies in the interferometric data. In fig. 3 therewrapped phase model (left), the residue image (centre) and a 3-D reconstruction ofthe unwrapped phase (right), are shown, respectively after the first (fig. 3a), the third(fig. 3b) and the fifth (fig. 3c) iteration of the described approach. The first iterationshows clearly the effect of the applied grid: the phase model is approximated by largeportions of flat surfaces (fig. 3a, on the right). As the model gets refined, more detailsare introduced in the model. As a consequence, fringes in the rewrapped phase imagebecome more and more complicated. By comparing the original interferogram (fig. 3d)and the rewrapped true phase after the fifth iteration (fig. 3c, on the left), it is possible tonotice the high similarity of fringes in the noise-free regions. On the contrary the resultis less reliable where inconsistencies are present but, in any case, errors are confinedin the noisy areas. In the residue image, after each iteration, the number of residualfringes decreases and they appear wider. This means that the significant information leftin the residue image is progressively reduced. Since after the fifth iteration some largefringes are still present (fig. 3c, centre), the algorithm should be iteratively applied tothe residue image in order to extract the remaining information and improve the finalphase model.4. – Discussion and conclusionsIn this paper a phase unwrapping method has been described which shows a goodaccuracy and a noticeable robustness to noise in reconstructing the true phase modelstarting from the interferometric wrapped values. The algorithm is based on the recon-struction of a global unwrapped phase function exploiting local information about theshape of the function itself. The approach is iterative in the sense that an improvementof the phase model can be obtained successively applying the algorithm to the residueimage: in this way it is possible to capture the details of the scene. Another importantcharacteristic of the algorithm is the possibility for the user to keep control over theprocessing during the iterations and to stop the algorithm when the residue image doesnot contain any fringes.Summarising, the proposed algorithm allows the user to unwrap the interferogramin an intuitive and to the point manner. The method limits error propagation, and132 G. BO and S. DELLEPIANEhas the advantage that no processing parameters have to be defined. Moreover, thestraightforward computation of the local phase guarantees reduction of computationtime in respect to other method present in literature.REFERENCES[1] Ghiglia D. C. and Romero L. A., J. Opt. Soc. Am., 11-1 (1994) 107-117.[2] Goldstein R. M., Zebker H. A. and Werner C. L., Radio Sci., 23-4 (1988) 713-720.[3] Pritt M. D., IEEE Trans. Geosci. Remote Sensing, 34-3 (1996) 728-738.[4] Tarayre H., Massonet D. and Sirat J. A., Proc. SPIE, Vol. 2548 (1995) 310-138.

Weighted multi-resolution phase-unwrapping method

Abstract

Similar works

Full text

Available Versions

Scientific Open-access Literature Archive and Repository

Archivio istituzionale della ricerca - Università di Genova