A method for the phase unwrapping in interferometric synthetic aperture radars (InSAR) and an algorithm for its implementation are proposed, which localizes an unwrapping error in a neighborhood of a point of discontinuity. The method includes the iterative phase discontinuity correction by implementation of the phase pseudo discontinuities of opposite directions. Application of the algorithm allows one to obtain a continuous phase function, which can then be converted into an absolute by the simple unwrapping by a linear path. The method was tested on models of typical discontinuities (isolated phase gap, phase dipole, phase aliasing) and on the interferogram obtained by the ALOS PALSAR. The method demonstrated the best results in comparison with traditional unwrapping methods, i.e. SHAPHU (MCF) and Region Growing ones. It is also shown that the method can be easily implemented in parallel processing systems.
Digital elevation models (DEM) and displacement maps are widely used in various scienti c and technical areas, i.e. cartography, geodesy, geology, environmental monitoring of mining areas, monitoring the transport communications, etc. [1{5]. Radar remote sensing is performed by interferometric SAR techniques (InSAR and DInSAR) that allows one to obtain both types of elevation data with semi-automatic data processing, which makes it very attractive for use in these tasks. However, the phase unwrapping stage of interferometric processing, which converts a relative phase de ned on the interval [ ; ], into an absolute phase, which is approximately proportionally related to the surface topography (for InSAR) or the relief displacement (DInSAR), is the obvious \bottleneck" of the whole radar interferometry. Existing algorithms for its solving are generally based on utilization of the optimization techniques (Minimal Cost Flow, Integer optimization), on search for the optimal integration path (Goldstein residue cut), on solution of large systems of equations (least square method), etc. Such techniques have low computational e ciency (typically quadratic complexity) and are di cult for parallel execution. Also, the most part of the existing methods is aimed to building the eld absolute phase, which is congruent to the relative phase eld. But it is actually useless in terms of side-looking radar geometry, which causes an irreversible damage to the local data parts.
Phase discontinuity is an element of the relative (wrapped) interferometric phase, which leads to dependence of the absolute (unwrapped) phase on the integration contour shape, and, so, the absolute phase cannot be restored uniquely. It is possible to identify the following elements of the discontinuity: two discontinuity points (so-called residues), where the cumulative sum on elementary path (4 adjacent elements) of the phase gradient is not equal to zero, and a discontinuity line, which virtually connects these points (Fig. 1a). It is possible to allocate all discontinuity points by calculating the residue function for the whole interferogram, but not the discontinuity line. This fact makes the solution of unwrapping problem ambiguous. Analysis of di erent interferogram types allows one to distinguish, at least, 3 types of phase discontinuities that may request di erent unwrapping techniques.
1. An elementary discontinuity is the phase discontinuity, which is caused by phase noise peaks. Such discontinuity covers two adjacent elements of the interferogram, and the their phase di erence exceeds . A simple unwrapping procedure (phase unwrapping by linear paths) passing through both these elements will result in the absolute phase error of the value 2 .
2. Phase discontinuity caused by layover. The echo signals layover occurs when the terrain slope angle exceeds the elevation angle of the SAR carrier, and the SAR echoes from the di erent surface elements return simultaneously and can not be resolved.
3. Phase discontinuity caused by aliasing. It occurs in high-slope terrain due to the discrete nature of the interferogram. In such scene elements, the adjacent interferometric fringes disappear that leads to a signi cant unwrapping error. 2
An ambiguity resolving method for phase
discontinuities
The layover discontinuity (type 2 in the classi cation above) is the most frequent
and complex discontinuity type, and, so, it's reasonable to explore it closely. Such
discontinuity may be easily simulated in a complex domain as a function like
I_(zm;n) = exp j arg
zm;n
zm;n
z01
zp1
;
(
Its structure may be better shown on a three-dimensional phase image. In
the case of single discontinuity point (Fig. 1b), the phase acquires a form of the
vortex evenly gluing the branches 3D-phase; and for two discontinuity points
the superposition of two distant oncoming vortices takes place, which forms a
jumper between unambiguous 3D-phase branches (Fig. 1c).
(
a b c
Let's use the circumstance that in the case of type 2 discontinuity the pair of opposite-signed residues in points z01 and zp1 it remains localized in some neighborhood of these points. Such discontinuity forms two arti cial phase vortices C_ 0(zm;n) and C_ p(zm;n) (pseudo-discontinuities) with centers at the points z01 and zp1 that have inverse directions a b
The SAR interferograms usually contains multiple chaotic located discontinuities, and, so, it is essential, at rst, to localize them by the residues function, and then form the product of the complex images of inverse vortices for each point
P_ (zm;n) = C_ 01(zm;n)C_ 02(zm;n) ::: C_ 0M
C_ p1(zm;n)C_ p2(zm;n) ::: C_ pN : (
Let us call P_ (zm;n) the inverse vortex phase eld. A dot product of complex interferogram and inverse vortex phase eld
I_c(zm;n) = I_(zm;n)P_ (zm;n)
(
On the basis of proposed phase ambiguity resolving technique, let us formulate an algorithm for the phase unwrapping, which would include the following steps:
1. generation of interferogram residues function | Rm;n.
If Rm;n = 0 8(m; n), then go to step 7;
2. detection of elementary discontinuities fWe(m; n)g according to the criterion of 8-neibourhood of two opposite-signed residues, and their correction by phase zeroing;
3. regeneration of interferogram residues function | Rmc;n.
If Rmc;n = 0 8(m; n), then go to step 7;
4. generation inverse vortex phase eld for remaining residues points fz01; z02; :::; z0M g and fzp1; zp2; :::; zpN g
P_ (zm;n) = exp j arg (zm;n (zm;n zp1) (zm;n z01) (zm;n zp2) ::: (zm;n z02) ::: (zm;n zpn) z0n) 5. interferogram correction with inverse vortex phase eld
I_(zm;n) ! I_(zm;n) P_ (zm;n); 6. regeneration of the interferogram residues function | Rmc00;n; If Rmc00;n 6= 0 8(m; n), then go to step 4, else go to step 7; 7. simple phase unwrapping by the linear path.
Steps 4{5 have the most computational complexity, but their performance
may be improved by the following measures.
;
(
2. The inversed vortex P 1(zm;n) for a single point can be computed once for the eld of 3M 3N size and saved into the processing device memory, and at, step 4, the vortex fragment of M N size with center at the point corresponding to the zero z0i or the pole zpi can be simply read from the memory.
3. Despite the algorithm is global, the inverse vortex phase eld is constructed by independent repeated complex multiplications (or additions in the relative phase domain) of the deterministic function, and, so, it can be easily implemented through a parallel execution of computing devices (within one iteration). Computations of residues function on step 1{3 and 6 are local and, so, can be distributed by di erent computational devices; computations for steps 4{5 must be made on the interferogram of the original size, but discontinuities may be passed in any order by di erent computational devices, and the partial results may be summarized.
The proposed unwrapping technique will further be called the inverse vortex phase eld method (IVPF). The proof of the algorithm convergence is not in the scope of this paper, but it should be noted that the algorithm does not diverge in any of the following cases. A similar approach to phase unwrapping in laser interferometry applications in a simpli ed form was previously proposed Aoki et al. [6] and further studied by Tomioka [8]. However, it does not obtain noticeable development due to higher complexity of shapes in the interferometry of the \small forms"; but it seems to be more applicable to the radar interferometry with its peculiarities [9]. 3
An e ciency analysis of phase unwrapping by the IVPF algorithm Let us use the following models of interferometric phase for comparative analysis of the IVPF algorithm e ciency.
1. A \lake" model, which simulates an uncorrelated SAR interferogram of M N size and includes an area with uniformly distributed phase (Fig. 3a). Such model simulates numerous elementary discontinuities.
2. A \Gauss hill" model (Fig. 3b), which simulates the phase discontinuity caused by layover.
3. A \steep slope" model (Fig. 3c), which simulates the phase discontinuity caused by aliasing.
E ciency is estimated by accuracy of the absolute phase restoration and by the required computer time. The major evaluation of the restoration accuracy is performed by the standard deviation of the simulated absolute phase and restored absolute phase = s
1 M N
X 1 m;n ^m;n m;n 2 :
However, the accuracy estimation based on the standard deviation criterion is insu cient, because in presence of propagating unwrapping error, the deviation becomes strongly dependent on the size of the interferogram. For larger interferogram sizes, the standard deviation may be small despite the fact that the restored absolute phase can have an obvious and signi cant damage. Therefore, let us introduce two additional accuracy criteria for the unwrapped phase: | stripe coe cient of the propagating error, which is equal to proportion of incorrectly unwrapped elements to the linear interferogram size 0 = lim
MN!!11
N (M; N ) pM N
; 1 = lim
MN!!11
N (M; N )
M N
:
where N is the number of incorrectly unwrapped elements. If the propagating
error is localized in a stripe of constant width, limit of (
The following phase unwrapping algorithms were researched:
(
Experiments were conducted on square form simulated phase models (M = N ). The results for 3 models are presented in Tables 1|5 and Fig. 4; the results for \lake" model were calculated for the scene part lying outside the damaged area. The following notes should be done according to the experiment results.
1. For model 1 (\lake"), the damaged area was unwrapped correctly only by two algorithms: the MCF and IVPF. The other two algorithms generate propagating error, and in the case of Region Growing algorithm, the error diverges with a 1=0.14. For the IVPF, a phase error around the damaged area occurs, which is rapidly decreasing with a distance from the center of the damaged area. The MCF algorithm in the majority of cases doesn't produce any error, but in one experiment a propagating error with 1=0.21 has appeared (Fig. 4b).
2. For model 2 (\Gauss hill"), none of the algorithms has recovered the absolute phase accurately. However the MCF and IVPF algorithms do not lead to the propagating errors. The IVPF algorithm distorts the neighborhood of the discontinuity such that the phase inside it tends to a zero mean value (Fig. 4c). The MCF algorithm connects the edges of the discontinuity region by the strips with approximately constant phase.
3. For model 3 (\steep slope"), none of the algorithm does not recover the absolute phase accurately. The MCF algorithm generates an unwrapping error; the IVPF distorts the discontinuity neighborhood, but demonstrates the lowest phase error.
4. The speed of the IVPF and MCF algorithms is determined not only by the interferogram size, but, also, by the number of discontinuities (\lake" model). For the size of 1500 1500 elements the IVPF wins the MCF in performance by 24%.
Algorithm
SU RG MCF IVPF
An experiment with ALOS PALSAR interferogram phase unwrapping was conducted with application of the MCF and IVPF algorithms (Simple unwrapping and Region Growing were useless here). The interferogram (Fig. 5a) has 13072 3600 elements and it was previously ltered by the Goldstein-Baran spectral lter. An accuracy estimation was performed by the inverse transformation of the reference DEM [10]. Both algorithms show comparable results in accuracy of the absolute phase restoration: 7.99 m for the MCF and 7.19 m for the IVPF, and the computational time is 923 sec for the IVPF and 3140 sec for the MCF. The Inverse vortex phase eld method (IVPF) and an algorithm for its implementation are proposed for the phase unwrapping in interferometric synthetic aperture radars (InSAR/DInSAR) data processing. Three test phase models for typical phase discontinuities were simulated, and the analysis of the algorithm e ciency was performed. Also, two additional criteria for accuracy estimation of the unwrapped phase were proposed taking into account the in uence of the propagating unwrapping error, which can't be estimated correctly by the standard deviation. It is shown that application of the IVPF algorithm does not lead to appearance of the propagating errors, and its accuracy is slight better than the same of the Minimum cost ow method (MCF/SNAPHU); but the computational speed is faster up to 3 times for the large scenes with numerous discontinuities. It is also shown that the algorithm allows parallel executing on multiple computing devices (within one iteration) for the further performance improvement.
Acknowledgments. The work was supported by Act 211 Government of the Russian Federation, contract 02.A03.21.0006.