INSAR UNWRAPPING ERROR CORRECTION BASED ON QUASI-ACCURATE DETECTION OF GROSS ERRORS ( QUAD )

Unwrapping error is a common error in the InSAR processing, which will seriously degrade the accuracy of the monitoring results. Based on a gross error correction method, Quasi-accurate detection (QUAD), the method for unwrapping errors automatic correction is established in this paper. This method identifies and corrects the unwrapping errors by establishing a functional model between the true errors and interferograms. The basic principle and processing steps are presented. Then this method is compared with the L1-norm method with simulated data. Results show that both methods can effectively suppress the unwrapping error when the ratio of the unwrapping errors is low, and the two methods can complement each other when the ratio of the unwrapping errors is relatively high. At last the real SAR data is tested for the phase unwrapping error correction. Results show that this new method can correct the phase unwrapping errors successfully in the practical application. * Corresponding author


INTRODUCTION
In order to acquire high accuracy monitoring results with Small Baseline Subsets (SBAS)-InSAR technology, various errors should be estimated and eliminated (Hanssen, 2001), among which the phase unwrapping error is the phase jump of the integer times of 2π due to the low coherence of the interferogram.The unwrapping error can seriously decrease the accuracy of InSAR results (Zebker et al. 1992).Pepe and Lanari (2006) presented an extension of the minimum cost flow algorithm for unwrapping of multitemporal differential SAR interferograms by assuming a temporal displacement model.Yu et al. (2013a) presented a fast phase unwrapping method for large-scale interferograms, which could obtain the unwrapped phase quickly and accurately.Liu et al. (2015) presented a cluster-analysis-based noise-robust phase-unwrapping algorithm.Compared with the conventional cluster-analysisbased method, this method improves noise robustness significantly.However, the unwrapping errors can hardly be solved completely (Yu. 2013).In the SBAS InSAR processing, a closed loop detection and correction of phase unwrapping error is often applied (Biggs et al., 2007, Wang et al., 2012).However, it still needs manual operation in some cases, and it can hardly operate automatically (Li et al., 2014).Lauknes et al. (2011) presented an L1-norm based method to accurately estimate the time series deformation by an iteratively reweighted least squares algorithm.The L1-norm based method suppresses the unwrapping error by decreasing its weight rather than correcting the error.A novel method, Quasi-accurate detection (QUAD), to correct the phase unwrapping errors in the interferograms will be proposed in this study, which was firstly proposed by Ou (1999) for gross error correction in geodetic measurements.

Principle of QUAD
The phase unwrapping errors occurred on interferograms results in closed residual errors, which can be regarded as gross error.So the QUAD method is designed to detect and correct the unwrapping error simultaneously and automatically under the frame of SBAS-InSAR method.For one generic pixel, the following observation equation can be formulated, where A is the design matrix with m n × dimensions and rank m ; 0 X , the cumulative deformation, is the n-dimensional unknown parameter vector; L is the unwrapped phase with m dimension, Δ is the true error with m dimensions.The formula (1) is rewritten as, where 0

X
∧ is the parameters to be estimated.V is the residual of the observations, which satisfies: where R is the adjustment factor, and , I is a unit matrix.However, due to the rank deficiency of R , the true error Δ cannot be estimated uniquely.According to the election group fitting principle, it is assumed that r quasi-accurate observations are selected, the corresponding coefficient matrix is r A .Then the additional condition equation is as follows, where According to formula (3) and ( 4), it can be derived that, ( ) Therefore, when the quasi-accurate observations are determined reasonably, the true error will appear as clusters as the unwrapping error is usually larger than other errors.Then the position of the unwrapping error can be easily determined.Assuming that z unwrapping errors are found, then the matrix R and L can be divided into four and two blocks as follows. , At last, based on the formula of gross error estimation in QUAD (Ou. 1999), the unwrapping error could be estimated according to the following formula.
where * z Δ is an estimated unwrapping error.

Processing Steps
The general processing step is presented with the simulated data, which are generated through the actual temporal and perpendicular baseline of 129 interferograms by setting thresholds as 200 days and 80 meters for 42 sentinel SAR data (shown in figure 1).The unwrapping error is identified and corrected pixel by pixel as follows: (1) The establishment of observation equations.
With the same pre-processing as SBAS, SAR data registration, interferograms generation and phase unwrapping are processed firstly.Then 129 observation equations will be established according to equation ( 2) for 129 unwrapped interferograms, and the residual V of each observation L is calculated under the Least Square norm (as shown in figure 2).(2) The preliminary selection of quasi-accurate observations.Based on formula (2), the adjustment factor matrix can be calculated as ( ) . Then some evaluation factors are given as follows Where ii r and ij r are elements of R .The observations can be divided into four categories, that is, ; 4) Type 2, the rest observations.The type 0 observations may contain unwrapping errors, the structure of the type 1 observations is not good, and the type 3 observations are considered to be good observations without unwrapping errors.In order to avoid rank-deficiency, more than m quasi-accurate observations should be selected.If the number of type 3 observations is less than m , type 2 observations with relative small residual errors should be selected as the quasi-accurate observations.Finally 2 m + quasiaccurate observations are selected in this study.
(3) The selection of final quasi-accurate observations.Based on the quasi-accurate observations determined in the second step, the true errors Δ can be obtained through formula (4) and ( 5).L and Δ are sorted according to the value of Δ from small to large, then the difference d Δ between the adjacent elements of Δ are calculated.In this step, gradual accumulation of quasi-accurate observations is adopted.The number of quasi-accurate observations is gradually added.Once the quasi-accurate observations changed, the Δ will be calculated and sorted again.If the number of quasi-accurate observations is equal to the position of maximum d Δ , the accumulation stops.It is assumed that k quasi-accurate observations are acquired finally.
(4) The estimation of unwrapping error.Based on the quasi-accurate observations determined in the step (3), the true error k Δ can be calculated by the formula (4) and (5).We set ( ) W is larger than 3, the corresponding observations are considered to contain the unwrapping error.We assume that z observations contain unwrapping errors are detected, and then the unwrapping error * z Δ can be estimated by the formula (7).Final, * z Δ are rounded to an integer multiple of 2π (as shown in the figure 2).We simulated the unwrapping error with an 18% uniform distribution, the 10m DEM error, the temporal decorrelation noise with a critical temporal baseline of 600 days, and an atmospheric delay with a standard deviation of 2mm.A is the residual error V obtained by LS; B is the true error Δ calculated by the preliminary quasi-accurate observations; C is the true error Δ calculated by the final quasi-accurate observations (step 3); D is the unwrapping errors estimated in the step 4; E is the simulated unwrapping errors.

Simulated Data
A set of simulated data are used to verify the effectiveness of this method by comparing with LS and the L1-norm methods.
The root-mean-square-error (RMSE) is used as a measure of goodness of fit, defined as follows.
where N the number of SAR images, i φ is the simulated time series deformation, * i φ is the real time series deformation.The new method can get the corrected unwrapping interferograms directly, so these interferograms will be used to obtain the time series deformation through LS then.
Three different deformation modes are designed to verify the appearance of three methods, that is, (1) Linear deformation with the deformation rate of 2 cm/year; (2) Linear deformation with the deformation rate of 2cm/year plus a sigmoidal drop of about 5 mm; (3) Linear deformation with the deformation rate of 2cm/year plus a periodic component of 5mm amplitude in one year.
A series of Monte Carlo simulation are implemented 500 iterations using different percentage of unwrapping errors with uniform distribution (0%, 3%, 6%, 9%, 12%, 15%, 18% and 21%).Two different atmospheric noise levels (2 and 4 mm) are set in this simulations.And the mean value of RMSEs of 500 Monte Carlo simulations is taken as the evaluation standard.Unwrapping errors are simulated by adding a phase of ±2π in some selected parts of the interferograms.Same simulation method as in Lauknes et al. (2011), the atmospheric delay is simulated with zero mean Gauss distribution on each SAR image and the temporal decorrelation noise with 300 days and 600 days critical baseline are simulated separately.
It can be found that both the QUAD and L1-norm method can reduce the influence of the unwrapping error effectively on the time series deformation when the ratio of unwrapping error is less than 21%.It should be noted that QUAD can get the correct unwrapped interferograms directly.However, the L1-norm gets the robust time series deformation directly, so extra processing is still needed to get the interferograms.In addition, with the increasing of the unwrapping error ratio, the RMSE of the LS method presents a rising trend with the fluctuation.This is because the effect of the unwrapping error on the final result is related to the number and the location of the unwrapping errors of the interferograms (figure 1).When the unwrapping error is located in the interferogram with less correlated observations, it will have a larger impact on the result.It can be deduced that baseline network is very important for the accuracy of monitoring results.Then we conducted a test with the unwrapping error ratio of 20%, 25%, and 30% respectively.The critical temporal baseline for 600 days, the atmospheric delay with 2mm noise level, and the linear with periodic deformation mode were set up for 500 times Monte Carlo simulations.As shown in Figure 4, it can be found that the RMSE of the two methods has different distributions, which indicates that two methods have different results of the unwrapping error correction.Based on the simulation above, we suggest that the unwrapping errors are completely identified and corrected (suppressed) when the RMSE of the corrected result is less than 3mm, and the unwrapping errors are partly corrected when the RMSE of the original result is 2mm larger than the RMSE of the corrected result.We calculated the number of times of complete and partial unwrapping errors correction by L1-norm, QUAD and either of two methods, as shown in the table 1.It can be found that these two methods can complement each other, and when the unwrapping error ratio is 25% and 30%, the number of times of unwrapping errors completely corrected by QUAD is larger than that completely corrected by the L1-norm method.When the unwrapping error ratio reaches 20%, the two methods can correct all the unwrapping errors for 285 times, which improves the number of complete correction for one method by 24%.It means that if one method fails, and another method may work.Table 1.The number of complete and partial unwrapping errors correction in the 500 simulations.

Real Data
The Maoxian landslide, Sichuan, China is selected as the real SAR SAR data test, which was occurred in June 24, 2017, and caused huge casualties and property losses.Forty-four archived Sentinel-1A and Sentinel-1B SAR data spanning from October 2014 to May 2017 are acquired over the study area.A multilook factor of 4 (4 pixels in range and 1 pixels in azimuth directions) is applied to generate SAR images.The baseline configuration of the interferograms is shown in figure 1. One-arc-second SRTM is used for the topographic phase removal.In the first step, 12 original interferograms with unwrapping errors in the landslide source area can be visually identified.
Then QUAD is to detect and correct the unwrapping errors pixel by pixel.Finally, all these 12 interferograms are completely corrected.We select 6 corrected interferograms and extract the phase before and after the correction along the profile AB, as shown in Figure 5.

CONCLUSION
A method for the unwrapping error automatic detection and correction based on QUAD is presented in this paper.Firstly, the effectiveness of the method has been validated with simulated data.Finally, the experiment of real SAR data shows that the new method can perform well in actual InSAR application.

Figure 1 .
Figure 1.The configuration of 129 interferograms.Every dot represents a SAR image, while each line corresponds to one interferogram.

Figure 2 .
Figure 2. The error distribution diagram of the interferograms.We simulated the unwrapping error with an 18% uniform distribution, the 10m DEM error, the temporal decorrelation noise with a critical temporal baseline of 600 days, and an atmospheric delay with a standard deviation of 2mm.A is the

Figure 3 .
Figure 3. Monte Carlo simulation results.The solid line corresponds to 2mm atmospheric noise while the dashed line corresponds to the 4mm atmospheric noise level.

Figure 4 .
Figure 4.The RMSE of the different methods with a ratio of 20%, 25% and 30% unwrapping errors, a critical temporal baseline of 600 days, and an atmospheric level of 2mm.

Figure 5 .
Figure 5. (a)-(f) are the original unwrapped interferograms, (g)-(i) the corrected unwrapped interferograms, (m)-(r) are crossections of the unwrapped phase before and after the correction along the profile AB.