LEAST SQUARES SOLUTION OF SMALL SAMPLE MULTIPLE-MASTER PSINSAR SYSTEM Lei Zhang(1) , Xiao Li Ding(1), Zhong Lu(2) (1)Department of Land Surveying and Geo-Informatics, The Hong Kong Polytechnic University, Hung Hom, KLN, Hong Kong, Email: [email protected]; [email protected] (2)U.S. Geological Survey, Vancouver, Washington, USA, Email: [email protected] ABSTRACT In this paper we propose a least squares based approach for multi-temporal SAR interferometry that allows to estimate the deformation rate with no need of phase unwrapping. The approach utilizes a series of multi-master wrapped differential interferograms with short baselines and only focuses on the arcs constructed by two nearby points at which there are no phase ambiguities. During the estimation an outlier detector is used to identify and remove the arcs with phase ambiguities, and pseudoinverse of priori variance component matrix is taken as the weight of correlated observations in the model. The parameters at points can be obtained by an indirect adjustment model with constraints when several reference points are available. The proposed approach is verified by a set of simulated data. 1. INTRODUCTION How to reduce or avoid the unwrapping errors is a challenge that all multi-temporal InSAR techniques [1-3] must overcome. In the application of such techniques we have found that for a set of multi-master interferograms with short baselines, there are usually a sufficient number of arcs at which the double-differential phase components contributed by topography errors and atmospheric artifacts appear to be very low and the relative deformation rate between two connected points is small. These arcs can be safely considered to be immune of phase ambiguities. If only these arcs are taken as observations for estimating the DEM errors and deformation rate, the complexity of parameter solution can be reduced significantly since there is no need of phase unwrapping. In this paper we propose a least squares based method which can select the arcs without phase ambiguities efficiently and solve the linear and non-linear parameters at temporarily coherent points (TCP) [4] reliably. We first use an efficient network construction strategy which performs _____________________________________________________ Proc. ‘Fringe 2009 Workshop’, Frascati, Italy, 30 November – 4 December 2009 (ESA SP-677, March 2010) the Delaunay triangulation locally to ensure that coherent points can be connected extensively without an increase of computational complexity. Then a least squares based model is proposed for parameter estimation. During the parameter estimation, an outlier detector is used to identify and remove the arcs with ambiguities according to the residuals raised by the first estimation. Considering the stochastic components of SAR acquisitions, we introduce the weights to interferometric phases at arcs by applying the law of propagation of variances [5-7]. Since it is likely that some of the interferograms are correlated each other, pseudoinverse of priori variance components is taken as the weight of the observations in our model. The parameters at TCP are obtained by spatial integration in the frame of least squares.. The proposed approach has been tested with a set of simulated SAR data to ensure the proposed algorithm functions as expected under controlled circumstances [5]. A validation has been carried out by comparing the estimated parameters with their true values. 2. MOLDING SAR INTERFEROGRAMS Considering J 1 SAR images which were acquired at the ordered time (t0 t1 t J ) , we generate I interferograms with short baselines (say, less than 150 m). In each interferogram i , the line-of-sight (LOS) displacement of TCP (l , m) can be given by a linear combination of the mean deformation rate (v j ) between time-adjacent acquisitions and the corresponding time span. Given two acquisitions, one is master ( M ) image and the other is slave ( S ) image and M was acquired later than S , i.e., t M i tSi . The LOS deformation during this period ( rli,m ) can be expressed as Mi rli,m r (t M i , l , m) r (tSi , l , m) (t j t j 1 )v j j Si 1 (1) Because Eq. (1) is a combination of LOS deformation li,m ,l ',m ' li,m hl ,m,l ',m ' i V wli ,m,l ',m ' estimated at the highest time resolution, it has a risk of over-parameterization, which should be carefully dealt with in real application. If a linear deformation rate ( v ) during the whole time span is assumed, then v1 v j v . Eq. (1) can also be tailored as any combination of deformation rates and time intervals in order to compare with field measurements which, for example, are preformed annually. The corresponding phase item is 4 i 4 M i rl ,m (t t ) v j Si 1 j j 1 j iV i defo ,l ,m i i i wli,m,l ',m ' atmo ,l ,m ,l ',m ' orbit ,l ,m ,l ',m ' noise ,l ,m ,l ' m ' (6) T V v1l ,m ,l ',m ' v l2,m ,l ',m ' v lI,m ,l ',m ' . Since atmospheric artifacts are strongly correlated in space, the differential atmospheric contributions relative to a close TCP will be extremely low (for points less than 1km apart is usually lower than 0.1 rad2 [7]). The differential where atmo (2) orbital component generally has a similar characteristic. Moreover, if both the connected TCP are not heavily i affected by decorrelation, noise ,l ,m ,l ' m ' will also show a low where is the wavelength, V is a vector of deformation rates, 4 i Ti and Ti is a vector of time combination whose variance as well. Therefore wli,m,l ',m ' can be safely taken as a elements correspond to the deformation rates within the respective time intervals. It has a form as follows Ti 0 t2 t1 t j t j 1 0 1J random variable with an expectation E ( wli,m,l ',m ' ) 0 . For arcs without phase ambiguities, the system of observations can be written as h A l ,m,l ,m W V (3) l1,m ,l,m The wrapped phase of a TCP with a pixel coordinate of (l , m) can be written as i i i i i li,m topo ,l ,m defo ,l ,m atmo ,l ,m orbit ,l ,m noise ,l ,m A (4) 1 l ,m (7) 2 l ,m T I l ,m T i i atmo ,l ,m is the atmospheric delay , orbit ,l ,m is the phase caused i by orbit error, and the noise ,l ,m is the additive noise term. The i topo ,l ,m term has a direct relationship with the height error hl ,m i 4 B ,l ,m hl ,m rli,m sin li,m 1 2 I i where topo ,l ,m is the phase related to the topographic residual, i topo ,l ,m l2,m,l ,m lI,m,l ,m W wl1,m,l ,m where is a vector containing phase differences at an arc in a total of I interferograms. A is design matrix including height-to-phase conversion factors and time combination matrix. W is a stochastic vector. 3. (5) hl ,m i l ,m where Bi ,l ,m is the local perpendicular baseline, rli,m is the range from master sensor to the target, li,m is the local incident angle. wl2,m,l ,m wlI,m ,l ,m LEAST SQUARES SOLUTION Introducing operators for the expectation ( E ) and the ), the functional and the stochastic model of dispersion ( D the system of observations can be expressed as follows: hˆ E = A l ,m,l ,m = - W ˆ (8) V D Q dd The phase difference between two TCPs located at (l , m) and (l ' , m' ) is given by where Q dd is an I I covariance matrix of observations. The function model reflects the linear or linearized relationship between measurements and the unknowns and the stochastic model describes the accuracy of measurements and their correlation between each other. In this section the weighted least squares estimator is used to solve the parameters in Eq. (7). 3.1. Priori variance components For the J 1 images considered, the variance matrix of i ) on the i-th TCP can be expressed as random noises ( Q noise i Q noise 2 noise i 0 2 noise i J 1 3.2. Least squares solution The least squares solution of observation equations is [7] hˆl ,m,l ,m 1 T dd T dd A P A A P ˆ V ˆ A AT P dd A1 AT P dd (12) v A AT P dd A AT P dd 1 where the circumflex ˆ denotes estimated quantities. The corresponding covariance matrices are (9) It is assumed that the noise is uncorrelated among SAR images. Since the I interferograms were combined from the J 1 images, the variance matrix of interferograms is given according to the law of propagation of variances as hˆ 1 D l ,m,l ,m Qxxˆ ˆ AT P dd A ˆ V ˆ Q ˆ ˆ A AT P dd A1 AT D (13) D v Qvv Qdd A AT P dd A AT 1 3.3. Outlier detection in Q noise DQnoise D T (10) where D is combination matrix indicating which pair of SLC images are used to generate the interferograms. As the priori variance, we assume that the TCP noise level at 2 H each SLC image is the same, i.e., Q1noise Qnoise Q noise , where H is the number of TCP. The variance-covariance (VC) matrix for the double-difference phases becomes in Q dd 2Qnoise 2 DQnoise D T (11) Considering the impact of atmospheric artifact on the double-difference phases is rather limited, we do not model it separately. Instead, we treat it as a stochastic component behaving like noise but with less magnitude. Least squares estimation has a basic assumption that all the gross errors and systematic effects have been eliminated before the adjustment is performed and only random errors affect the data. However during the initial least squares estimation all arcs including those at which there are phase ambiguities in some interferograms are taken as observations for the model. Herein we need to remove these unwanted arcs, i.e., outliers. To detect outliers in the observations, a statistical test of the estimated residuals is commonly used [5]. However the hypothesis testing is usually performed in a low efficiency. Considering our case that outliers mean double differences on arcs with phase ambiguities ( N 2 , N 1,2, , n ), which would cause the magnitude of residuals to grow dramatically, here we use a simplified detector for outliers [11] Max ( vi ) c Max ((Q dd )ii ) 2 Max((Q ˆ ˆ )ii ) The weights ( P dd ) for the double-difference observations can be obtained by taking the inverse of the VC matrix. However, it is possible that interferograms are correlated. Therefore, the VC matrix may not be a full rank. In this case we use the pseudoinverse of VC matrix instead, i.e., P dd (Q dd ) , which can be obtained by the singular value where Max () stands for getting the maximum value from the vector or matrix, the constant c may be chosen as 3 or 4. If the Eq. (14) holds true, the i-th observation is an outlier at the 95% confidence interval. decomposition (SVD) [8]. 3.4. Final solution (14) After removing the arcs possessing phase ambiguities, the least squares adjustment is performed again to estimate the parameters at arcs. Compared with the integer least squares estimator and the maximization of the ensemble coherence, our proposed method can determine the DEM error and difference deformation rates at arcs more efficiently and reliably since there is no need to perform a search of phase ambiguities in the solution space. Once the parameters at arcs are determined, the parameters at points can be obtained by spatial integration, which can also be performed under the least squares framework. The relationship between arcs and points can be expressed by a design matrix U . Let the i-th point be the reference point with known parameters ( Ri ), maximum magnitude of 72mm/year (Fig.2 (b)) to test the performance of the model for areas undergoing relatively large deformation. In order to get the phase differences at point pairs, we then adopted the local triangulation strategy to generate the network. The TCP were connected into 20091 arcs (Fig. 3(a)). The longest arc has a distance of 1454 m. By applying the outlier detector, 9711 arcs were removed out from the network (Fig. 3(b)). The number of arcs with phase ambiguities is 9535, all of which were successfully identified and it also means that 176 arcs were misidentified. After the outlier detection, the least squares estimator was performed again. and its corresponding column in U is Si , the observation equations are L UX Si Ri (15) where L is the parameters at arcs, and X is the parameters at points without the reference one, i.e. X x1 x2 xi 1 xi 1 xH (16) By introducing LL with LL L Si Ri , we obtain LL UUX (17) where UU is an updated design matrix in which the i-th column has been removed. The least squares solution can be obtained as X UU T PUU UU T P 1 (18) where P I . 4. VALIDATION WITH SIMULATED DATA A dataset consisting of 21 simulated C-band images have been used to validate the proposed approach. The advantage of using simulated data is that the estimated parameters can be compared with their true values which are often not known in case of real data [4]. During the simulation, the o standard deviation of the random noise is set to 15 for all images and the atmospheric phase is simulated using fractal surfaces with a dimension of 2.67. The details of simulation can be found in [4] and [9]. Examples of the simulated noise and atmospheric signal are shown in Fig.1. From the 21 images, we produced 44 interferograms with perpendicular baselines of less than 150m and temporal baselines of less than two years. Within an area of 5×5 km2 , 1,500 TCP are selected. The simulated DEM errors at TCP are shown in Fig. 2(a). We simulated the linear deformation signal with the Figure 1. Examples of simulated signals. (a) is the simulated noise in a certain SLC image. (b) is the simulated atmospheric artifact. The unit is rad. Using the priori VC matrix for the double-difference phase observations, we can get the VC matrix of the estimated parameters (the DEM error and deformation rate (differences)), see Eq. (13). 2.7062 0.0914 Qh ,v 0.0914 0.0188 (19) Here, the DEM error is first estimated in meters and deformation rate is the second parameter in mm/year. From the Eq. (19) we can find that the standard deviations of the estimated difference parameters are thus 1.64 m and 0.137 mm/year. Once the double difference parameters at arcs are determined, the parameters at points can be obtained by spatial integration (one reference point assumed). Fig. 3 shows the adjusted residuals. Figure 2. Parameters to be estimated. (a) DEM errors. (b) Linear deformation rate. The statistics of the errors on the estimated parameters is reported in Tab. 1. The standard deviation from the propagated VC matrix is given in parentheses. It can be found that the estimated DEM error is not as precise as that estimated by PSInSAR. It is mainly due to the fact that only interferograms with short baselines are selected for the model leading to the coefficient corresponding to DEM error in design matrix not sensitive to the errors. Table 1 Statistics of errors on the estimated parameters compared with their true values. The standard deviation from Eq. (13) is given in parentheses. Figure 3.(a) Coherent points are connected by local triangulation;(b) TCP network after removing arcs with phase ambiguities detected by the outlier detector from (a) min max mean std DEM error(m) -8.1 2.4 -2.6 1.72(1.64) Linear defo. Rate(mm/y) -0.45 0.41 -0.008 0.164(0.137) (a) thanks also go to Dr. B. Kampes for releasing the demo codes of STUN algorithm with [1]. REFERENCES (b) Fig.3 (a) the residuals of estimated DEM errors at TCP compared with their true values; (b) the residuals of estimated deformation rate at TCPs . 5. CONCLUSIONS A least squares model for the parameter estimation from multi-temporal SAR data was developed and a validation test was carried out using simulated data. Results shown in terms of statistical deviations of the estimated parameters, i.e., the DEM error and the linear deformation rate with respect to their true values indicate that the proposed model works well in areas either with moderate ground deformation rates or with complex deformation pattern. In the simulation the point density is 60 TCP/km2. It should be noted that the proposed model can also have good performance with sparse TCP density under the condition that enough arcs can be constructed. ACKNOWLEDGMENT The authors would like to thank Prof. R. Hanssen and Dr. B. Kampes for sharing the InSAR simulation package. Many [1] Ferretti, A., Prati, C., Rocca, F., & Elettronica, D. di (2000). Nonlinear subsidence rate estimation using permanent scatterers indifferential SAR interferometry, IEEE Trans. Geosci. Remote Sensing. 38(5), 2202-2212. [2] Lanari, R., Mora, O., Manunta, M., Mallorqui, J., Berardino,P., & Sansosti, E. (2004). A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms. IEEE Trans. Geosci. Remote Sensing. 42(7), 1377-1386. [3] Hooper, A., Zebker, H., Segall, P., & Kampes, B. (2004). A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers. Geophys. Res. Lett, 31, 611-615. [4] Zhang, L., Ding, X.L., & Lu, Z. Ground settlement monitoring from temporarily coherent points between two SAR acquisitions. ISPRS-J. Photogramm. Remote Sens., submitted. [5] Kampes, B. (2006). Radar Interferometry: Persistent Scatterer Technique, Springer. Dordrecht, The Netherlands [6] Teunissen, P. J. G. (2000). Adjustment Theory: An Introduction, 1st ed., Delft Univ. Press. Delft, The Netherlands. [7] Koch, K.(1999). Parameter Estimation and Hypothesis Testing in Linear Models, 2nd ed., Springer. Berlin, Germany. [8] Williams, S., Bock, Y., & Fang, P. (1998). Integrated satellite interferometry: Tropospheric noise, GPS estimates and implications for interferometric synthetic aperture radar products. J. Geophys. Res.-Solid Earth, 103(B11), 27051-27067. [9] Rao, C., & Mitra, S. (1971). Generalized Inverse of a Matrix and its Applications, J. Wiley. New York. [10] Hanssen, R. (2001). Radar interferometry: Data interpretation and error analysis, Kluwer Academic Press. Dordrecht, Netherlands. [11] Jia, P. Z., & Tao, Z. Z. (1984). Optimal Estimation and its Applications, Science Press. Beijing, China.
© Copyright 2024