1. science
2. mathematics
3. statistics

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: l.zhang@polyu.edu.hk; lsxlding@polyu.edu.hk
(2)U.S. Geological Survey, Vancouver, Washington, USA, Email: lu@usgs.gov
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 1J
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 Bi ,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 A1 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 A1 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 
Qh ,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.