COMPARISON OF TWO PHASE UNWRAPPING METHODS FOR PSP TIME SERIES ANALYSIS

Phase unwrapping is a key issue for PSP (Permanent Scatters Pairs) time series analysis. In this paper, the LLL method and grid search method are applied for determining elevation error and velocity difference in LOS for each PSP, and the RMS is used to assess the correctness of the solution. The Delaunay triangulation is used to formulate PSPs, which are used to unwrap elevations and velocities of all PSs with respect to a given reference point. The simulation examples with real SAR satellite orbits show that if the atmosphere noise is less than 2 radians and the random noise is less than 10 degrees, LLL method can obtain correct solution while grid search can obtain solution with a small error probability.


INTRODUCTION
Differential synthetic aperture radar interferometry (DInSAR) can detect large-scale surface deformation with millimeter accuracy (Goldstein et al.1989).DInSAR is widely used in the field of earth sciences such as monitoring deformation of seismic fault, ground subsidence, volcanic activity, land slide and so on (Colesanti et al. 2006;Zebker et al. 1994;Motagh et al. 2007;Pepe et al. 2008;Prati et al. 2010).
Time series analysis of interferometric phases provides a new idea to avoid the temporal decorrelation and atmospheric inhomogeneities during the DInSAR procedure (Ferretti et al. 2001;Kampes 2006).The time series method concentrates on those points maintaining good coherence during a long time called coherent targets or persistent scatters (PS).Based on the credible interferometric phase information of PSs, subsidence information can be obtained accurately (Ferretti et al. 2001;Kampes 2006).
Phase unwrapping is a key issue for PSP (Permanent Scatters Pairs) time series analysis.Phase obtained from the interferogram is wrapped, varying in [-π, π).The procedure of finding the lost 2kπ, getting the continuously changing real phase, is called unwrapping.Once the wrapped phase is unwrapped, the change of range in Line of Sight (LOS) can be obtained, then the height and subsidence velocities of PS points can be obtained.
Phase unwrapping appeared in the late 1960s, which is in one-dimensional.Since the late 1970s, two-dimensional phase unwrapping gradually flourished.Especially in 1988, Goldstein et al proposed branch-cut method (Branch-Cut), since then, various unwrapping algorithm emerged.These algorithms can be roughly divided into two categories: the first category is based on the path tracking algorithms, mainly including branch-cut method, quality map guidance algorithm, mask cutting (Mask-cut) algorithm and Flynn minimum discontinuity method, etc; the second category is based on the minimum norm algorithm, including non-weighted least-squares algorithm, preconditioned conjugate gradient (PCG) algorithm, weighted multi-grid algorithm, minimum LP norm algorithm, etc. Path tracking method can get the unique integer solution with high precision, but sometimes its continuity is so poor that there may be some wrapped islands; minimum norm unwrapping method has a strong continuity, but it may cause gross error be distributed to every point.
In this paper, we focus on the phase unwrapping of PSP time series by using the LLL method and grid search method.At first, the methods of the LLL and grid search are introduced, the simulation examples are then used to demonstrate the effectiveness of the methods, and some conclusions are given at last.

Basic principle of LLL algorithm
The LLL algorithm, proposed by Lenstra, Lenstra and Lovász (Lenstra et al. 1982) as a lattice space compute algorithm, can effectively reduce the search radius and narrow the search area.It has been applied to a variety of aspects such as integer programming (Finck et al. 1984), cryptography (Lagarias, 1984;Coppersmith, 1997) and number theory (Kaltofen et al. 1991 ).In GNSS and InSAR, a mixed integer least squares problem has to be solved in order to fix integer ambiguities and calculate positions.The basic idea of LLL algorithm is to separate the unknowns into integer part and real part, so as to solve integer first and float second.As a result, it is proved that LLL algorithm can provide fast and numerically reliable routines to a mixed integer least squares problem.
First, let's state our problem, If z is fixed, there must be a appropriate that enable the first term in (4) is zero, therefore the problem can be decomposed into the following two small problems: A. An ordinary integer least squares problem to calculate ̂ min (5) Specifically, a reduction algorithm and a search algorithm are presented to get the integer z which satisfies (5).With z known, (2) becomes a least squares problem.
B. With ̂ back into (4) and let the first term be zero, we get from For simplicity, we note the above problem A as: Obviously, y is a known vector, z is the least squares solution required, Bz is a vector in the grid, thus seeking the solution in (7) is of searching the nearest grid vector to y in fact.This is a CVP(Closest Vector Problem) problem and has been proven to be an NP-hard problem.To make the search process simple and efficient, a lot of reduction methods are proposed.The LLL method is an outstanding one of them.Reduction: First, with the minimum main-element method, do QRZ decomposition to matrix B to turn it into an upper triangular matrix R and an orthogonal matrix Q.Second, reduce the non-diagonal elements in R using integer Gaussian transform, so as to remove correlation and enable efficient search.Third, rearrange columns using minimum main elements principles to meet the LLL conditions.Search: After reduction, we need to search the optimal integer solution to satisfy min .Given a threshold β, assuming optimal integer solution z satisfies: (8) This corresponds to search the optimal solution within an ellipsoid.Then decompose R[r ij ] into first n-1 order sub-matrix and the last line, y into the (n-1) dimensional sub-vector and the last element.So, To satisfy ( 8), there should be: : , (11) is then an n-1 dimensional integer least squares problem, and the corresponding search radius is .The integer solution to (10) falls in / , / .Recursively using this algorithm, we can solve the upper triangular integer least squares problem.
Once the unknown integer solution ̂ is solved, we can use the following upper triangular matrix to solve the k corresponding real parameters: ̂ , where 1, ,1

LLL method used for InSAR
Assume that there are N SAR images of the same area, obtained at time t 1 ,⋯, tn.One image is chosen as the master and the other N-1 images are chosen as the slaves, so as to take N-1 interferograms.The unwrapped phase between pixel i and pixel j in the interferogram k is represented as where v and h are the relative displacement rate and relative error of elevation respectively.β changes with vertical baseline and α changes with time baseline.a0 is initial phase difference due to hardware defects.z is the number of unknown circles and δ is noise corresponding to decorrelation error, unmodeled nonlinear settlement, hardware defects and so on.Note that the atmosphere phase has been International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-7/W1, 3rd ISPRS IWIDF 2013, 20 -22 August 2013, Antu, Jilin Province, PR China considered to be spatial correlated which can be greatly reduced through near points differential We can assume relative error of elevation h , relative displacement rate v , an initial phase a0 and certain noise δ at the point i and j, generate a simulation phase with (12).Later wrap the phase and treat it as double differential phase observations.The LLL algorithm is then used to solve the integer least-squares problem for the parameters of v , h , a0 and ambuguity integer z.The reliability and efficiency of the solution is done by comparing the solution obtained with the known simulated values.
Suppose there are k +1 SAR images, generating k interferogram images.As one point in one interferogram has one ambiguity integer z, together with v , h , a0 , there are k+3 unknowns.As there are only k observations, hence the number of error equation is less than the number of unknowns and the coefficient matrix is rank-defect.For obtaining least-squares solution, we assume three initial values for the three unknown parameters as: to make the coefficient matrix full-rank.
Of course, these initial values are unreasonable, but we can update these initial unknowns by iterations.The new observation equation can be expressed as follows: where A1 has k rows, its three columns are the elevation-phase conversion factor, the sedimentation-phase conversion factor and 1. A2 is a 3*3 identity matrix; Real unknown x represents v , h and a0 .B1 has k rows and it's an identity matrix times 2π.B2 is a zero matrix.If signal to noise ratio of observations is high, the best solution can be found with a few iteration loops.

Grid search algorithm
The basic principle of multi-parameter grid search method is to divide possible variable range of each unknown parameter into a series of discrete values, ranging from small to large.By calculating the target function values of each possible parameter combinations, the parameter combinations with the minimum target function value is defined as the optimum solution.The grid search method can guarantee the solution is the global optimal, however, its accuracy is determined by limited search steps, and the computer load will increase dramatically when parameters are larger and the searching step becomes small.Grid search method can be used to solve the relative displacement rate and the error of elevation among PSPs.The target function of the grid search is defined as the root mean squares (RMS) of residuals of the warped double phase differences.So there is no need to do phase unwrapping for the grid search method, thus error during unwrapping process is avoided.What's more, grid search method is relatively simple and easy to be implemented.

Program reliability test
The land subsidence of Shanghai is currently stabilized, sedimentation rate in general is less than 10mm/y.Assume that there are only small amount of high buildings in place of interest, so the relative error of elevation ranges from -30m to 30m.
First, generate 61 * 21 * 31 =39711 samples.For each group (v , h , a0 ), a simulated interferometric phase time series is given with certain noise.Then the reliability of the program will be evaluated by whether the ambiguity integer z is solved correctly.It is proved that for all these 39711 combinations, all z obtained by LLL are correct and the solution is very close to the true value.Figure 1 shows corresponding phase RMS of each combination test.The histogram of difference between the simulation relative sediment velocity and LLL solution is shown in Figure 2 while the histogram of difference between the relative simulation elevation error and LLL solution is shown in Figure 3.As a result, the phase standard deviation of all combinations are less than 2δ (δ is the phase error.)and it takes few minutes to work out all the 39711 combination tests.All the differences between the simulation sediment velocity and LLL solution are less than 0.05mm/y while most of the differences between the simulation elevation error and LLL solution are less than 0.05m.Thus the program is reliable and efficient.

Grid search program test
Assume the searching step of sediment rate is 1mm / y while the step of elevation error is 1m.On the assumption that the ground subsidence is linear, then for arbitrary v , h , a0 combinations ( v : ( -30mm/y, 20mm/y ) , h :(-30m,30m), a0 =0), the histogram of difference between the simulation relative sediment velocity and grid search solution is shown in Figure 4 while the histogram of difference between the relative simulation elevation error and grid search solution is shown in Figure 5.
From Figure 4 and Figure 5, we can see that most of the differences between the simulation sediment velocities and that that of grid search solutions are less than 1mm/y while most of the differences between the simulation relative elevation errors and that of grid search solutions are less than 2m.Thus the grid search program is reliable and efficient.

Delaunay triangle net test with two methods
In order to approximate the real data, the ENVISAT ASAR data from European Space Agency is used.There are 31 ASAR images of Shanghai area from February 2007 to May 2010.Image of August 2008 is selected as master image, resulting in 30 interferograms with time and space baselines.The incident angle θ and distance from antenna to scattering target are the average of 31 images.Figure 6 shows the distribution of baseline (m).Gaussian noise is considered and the standard deviation of noise is simulated less than 10 degrees.The max atmosphere phase is set as 2radian and maximum DEM error 20m, maximum sediment 20mm.Figure 7 shows the Delaunay triangular net.54 PSPs are formed.Table 1 shows the basic parameters in the simulation.Figure 8 shows the elevation difference and sediment rate difference of all PSPs calculated by the LLL algorithm and the grid search method.Figure 9 shows the absolute elevation and sediment rate of PS points calculated by the two methods with respect to the last PS point.(Assume that the absolute elevation and sediment rate of 20th point are zero.)Table 2 shows the comparison of true relative elevation errors and sediment rate differences of 54 PSPs and their corresponding solutions.8, among 54 PSPs, the LLL method can always obtain correct solution with decimeter elevation difference and 0.1mm/y sediment rate difference while grid search method can almost get correct solution except for one wrong solution.As a whole, the accuracy of the LLL method is better than that of grid search method.See from Figure 9, all the differences between true absolute elevations and that calculated by the LLL method is less than 1m while the differences between true absolute elevations and that calculated by the grid search method is mostly less than 5m.All the differences between true absolute deformation rate and that calculated by LLL method is less than 1mm/y while the differences between true absolute deformation rate and that calculated by grid search method is less than 3mm/y.Taking the accuracy of the existing InSAR technique into account, both results are acceptable.
From Table 2, we can see that both LLL method and grid search method can reach accurate solutions.However, the grid search method makes great mistake at dealing with the 14th PSP.In fact, the grid search method cannot give integer ambiguity because its assumption on search is that the phase difference between arbitrary point pairs is within [-π,π).For example, if the true phase difference value is 2.8rad, with certain noise (e.g.,0.4rad) added, the simulation phase difference observation will be 2.8+0.4-2*π=-3.04rad,unwrapping using -3.04rad will lead to great error, thus double phase this -3.04rad is mirrored to be 3.04rad which is closer to the true value.For the 14th PSP in Table 2, after the preprocess, grid search method gives correct solution (dH=-4.2m,dv=19mm/y) finally.

CONCLUSION
This paper presents two methods for PSP (Permanent Scatters Pairs) time series analysis, namely, LLL algorithm and grid search method.31 SAR images of the Shanghai area from ENVISAT are used for taking simulation test.The results show that both methods can obtain correct elevation and sediment rate of PSPs with acceptable accuracy if the noise including atmosphere noise and random noise is less than 20 degrees.The LLL method has greater dependence on initial unknowns but can obtain accurate solution if iteration process is convergent.The grid search method is relatively simple but is time-consuming if search step is small, and preprocess should be take for observations which values are close to -π or π.We plan to combine these two methods together for PSP unwrapping in the future.

Figure 1 .
Figure 1.Phase RMS of each combination

Figure 3 .
Figure 3. Histogram of difference between the simulation elevation error and LLL solution

Figure 4 Figure 5 Figure 6
Figure 4 Histogram of difference between the simulation sediment velocity and grid search solution

Figure 7
Figure 7 The Delaunay triangle net

Figure 8
Figure 8 Elevation difference and sediment rate difference of arcs calculated by two methods before adjustment

Table 1 .
The basic parameters in simulation Overview of parameters:

Table 2 .
Comparison of true relative elevation errors and sediment rate differences of PSPs with solutions Comparison of true v , h with solutions by the two methods Start point ID End point ID Red row in the table shows the only wrong solution in the grid search method.h represents the relative elevation error of PSP and v represents sediment rate difference of PSP.International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-7/W1, 3rd ISPRS IWIDF 2013, 20 -22 August 2013, Antu, Jilin Province, PR China