A parallel implementation of the inverse vortex phase eld algorithm is proposed for unwrapping the large-size interferograms obtained by the space synthesized aperture radars. It is shown that the parallelization of calculations in solving the unwrapping problem by this method is possible due to the fact that the phase discontinuity compensation operations are linear. So,the computations can be carried out in any sequence and they can be arbitrarily distributed among the individual computing devices. The method of interferogram splitting into blocks of arbitrary size (not exceeding the memory size of computing devices) is also proposed, and it allows solving the problem of unwrapping for the large-size interferograms. The results of processing the ALOS PALSAR interferograms with the size of 14x3 thousands elements on a quad-core computing device are presented.
Digital elevation models (DEM) and displacement maps are widely used in various scienti c and technical areas, such as cartography, geodesy, geology, environmental monitoring of mining areas, monitoring the transport communications, etc. [1{5]. The radar remote sensing is performed by interferometric SAR techniques (InSAR and DInSAR) that allows one to obtain both types of elevation data with the semi-automatic data processing, which makes it very attractive for use in these tasks. However, the phase unwrapping stage of interferometric processing converts a relative phase de ned on the interval [ ; ] into an absolute phase. The latter is approximately proportionally related to the surface topography (for InSAR) or the relief displacement (DInSAR). So, the unwrapping stage is the obvious \bottleneck" of the whole radar interferometry. Existing algorithms for its solving are generally based on application of the optimization techniques (Minimal Cost Flow, Integer optimization), searching 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 the parallel execution. Also, the most part of the existing methods is aimed to building the absolute phase eld, which is congruent to the relative phase eld. But it is actually useless in terms of the 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. 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. As it was shown earlier [11], an inverse vortex phase eld algorithm for phase unwrapping allows one to obtain a stable continous solution of the problem, and it may be suitable for the parallel and block implementation for the large-size data processing. 2
Inverse vortex phase eld method (IVPF) uses an arti cial phase discotinuity (phase singularity) of the opposite sign (pseudo-discontinuity) to compensate the interferogram phase discontinuity. The most widespread type of the phase discontinuity is the layover, which can be described as an angle of rational function 1 in a complex domain
I_(zm;n) = exp j arg
zm;n
zm;n
z01
zp1
;
;
(
Such discontinuity may be compensated by two arti cial phase vortices C_ 0(zm;n) and C_ p(zm;n) with centers at the points z01 and zp1 that have inverse directions C_ 0(zm;n) = exp j arg
zm;n C_ p(zm;n) = exp fj arg [zm;n 1
z0 zp]g ; where C_ 0(zm;n) and C_ p(zm;n) are the inverse phase vortices.
A new interferogram I_c(zm;n) is formed after the compensation of all phase discontinuities
I_c(zm;n) = I_(zm;n)C_ 0(zm;n)C_ p(zm;n):
As it was determined earlier by experiments [11], there is no neccessity to identify pairs of discontinuity points, and all such points may be processed seperately in any order. The SAR interferograms usually contain the multiple chaotic located discontinuities. So, it is essential, at rst, to localize them by the residues function, and then form the product of the complex images of the 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 ; (
I_c(zm;n) = I_(zm;n)P_ (zm;n) forms the corrected interferogram I_c(zm;n), where the phase ambiguity is not obligatory fully resolved because the inverse vortices may lead to occurrence of new phase discontinuities or to moving it into a new position. The procedure should be iteratively repeated until all ambiguities to be fully resolved.
For elementary discontinuities, such correction is unwanted because the application of the inverse vortex will here lead to additional distortion of the unwrapped phase. On the other hand, a simple zeroing of such discontinuities instead of inverse vortex correction does not produce an unwrapping error and allows one to increase the computational speed.
So, the algorithm for the phase unwrapping will include the following steps: 1. generation of interferogram residues function | Rm;n; If Rm;n = 0 for all (m; n), then go to step 7;
2. detection of elementary discontinuities fWe(m; n)g according to the criterion of the 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 for all (m; n), then go to step 7;
4. generation the 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 by the 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= 0f orall(m; n), then go to step 4, else go to step 7;
7. simple phase unwrapping by the linear path.
(
Steps 4{5 of the algorithm have the most computational complexity, but their performance may be improved by the following measures.
1. Multiplication of the complex functions I_(zm;n) and P_ (zm;n) may be unambiguously replaced by summation of their arguments, and the arguments can be evaluated once before the rst iteration.
2. The inversed vortex P 1(zm;n) for a single point can be computed once for the eld of 2M 2N size and saved into the processing device memory, and on step 4, the vortex fragment of M N size with the centre 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. So, it can be easily implemented by the a parallel execution of computing devices (within one iteration).
Computations of residues function on step 1{3 and 6 are local. Therefore, it 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.
But the following problem remains. Usually, raw SAR interferograms (without data cropping and multilooking) have sizes of thousands elements by each side, or more than 100 million (sometimes, more than billion) pixels in the complex single-precision format (more 1{5 GB or more). During the algorithm execution, these data may be reproduced for several times, requiring big RAM volumes. So, to make the IVPF applicable to the large-size interferograms, the following implementation strategy (the block-IVPF algorithm) may be o ered.
1. Splitting the whole interferogram by the blocks Bi of relatively small size (for example, 1000 1000 elements), that are easily processed by a single processing unit. The right-sided and bottom blocks may have the distinctive size. 2. Localization of the discontinuity points in the Bj th block.
3. Processing the Bi block:
{ generation of the inversed vortex of the Bj block (double-sized respectively to
the block size) in the borders of Bi block by the suitable part of (
j pairs are processed. 5. Block merging.
6. Localisation of the discontinuity points. If any discontinuity points remain, then return to the step 1, else, stop the algorithm.
Execution of step 3 may be implemented on the parallel computing devices, because, as it mentioned earlier, all operations in the algorithm are commutative and may be executed independently on each other. For the parallel execution, it is required to tramsmit the block data (or group of blocks data) and all discontinuity points coordinates to the elementary device and to transmit the processed block back to the central device for merging. Experimental imagery data was received by the ALOS PALSAR radar and are represented by an interferometric pair of SAR images for the polygon contained surfaces of di erent types ( elds, forests, urban, and mining areas). Data samples were obtained in dual polarization mode, and the HH-polarization was used. The phase eld has an ambiguity height of 17.2 m (at the near stage); the spatial sampling interval was 15:0 3:1 m. Two data subsets were used for the unwrapping experiments 1. scene with the size of 3000 3000 samples; 2. scene with the size of 14000 3000 samples (Fig. 2a).
The block-IVPF algorithms was implemented in the MATLAB environment. It includes the main program and two sub-programs, which have 65 code strings in total, only. The codes were executed on a system with the Intel Core-i5 processor with 2.67 GHz clock frequency, the RAM size was 12 GB.
The rst scene was processed by the both initial and block versions of the IVPF, the interferogram was splitted into 16 blocks in the block version. The results turned out to be equal, the computational time was 420 sec for the initial version and 730 sec for the block version. The last was slower, because only single computing device was used. Both algorithms converged in 6 iterations. The accuracy of the unwrapping result was earlier proved (for the initial version) in [9, 10]. The second scene (Fig. 2b) was processed by the block-version alogorithm only, because the initial version, as well as the Minimum cost ow and Region growing algorithms (in SARscape implementation) fail causing \Out of memory" error; the computational time was about 6 hours, the algorithm also converged in 6 iterations. 5
The block implementation for the Inverse vortex phase eld algorithm for the InSAR data processing was proposed. Two sets of the ALOS PALSAR data were processed, and it was shown that the algorithm gives the same results, as the initial algorithm version. The algorithm allows one to process large data sets, which are unsupportable by the initial algorithm and some other processing algorithms. The algorithm has very compact code and is suitable for execution on parallel computing systems. It is because the algorithm uses the commutative computational operations, and the data blocks may be randomly distributed between the computational devices.
Acknowledgments. The work was supported by the Act 211 Government of the Russian Federation, contract A 02.A03.21.0006.
a b