on the use of anisotropic covariance models in ... - Semantic Scholar
ON THE USE OF ANISOTROPIC COVARIANCE MODELS IN ESTIMATING ATMOSPHERIC DINSAR CONTRIBUTIONS A. Refice, A. Belmonte, F. Bovenga, and G. Pasquariello ISSIA-CNR, Via Amendola 122/d, 70126 Bari, Italy; E-mail: [email protected]
ABSTRACT When studying geophysical processes through Differential SAR Interferometry (DInSAR), it is often necessary to estimate and subtract signals due to atmospheric inhomogeneities. To this end, stochastic models are often used to describe atmospheric phase delays in DInSAR data. As a first approximation, these can be modelled as isotropic, though this is a simplification, because SAR interferograms often exhibit anisotropic atmospheric signals. In view of this, it is increasingly advocated the use of anisotropic models for atmospheric phase estimation. However, anisotropic models lead to increased computational complexity in estimating the correlation function parameters with respect to the isotropic case. Moreover, performances can degrade when use is made of interferograms with only a few sparse points usable for computations, such as in the case of persistent scatterers interferometry (PSI) applications, especially when this estimation has to be done in an automated way on several tens of interferograms. In the present work we critically analyse the main aspects connected with the use of anisotropic models for DInSAR atmospheric delays, and we evaluate the advantage given by anisotropic modeling of atmospheric phase in the case of sparse grids of points. Through analysis of APS simulated data, we observe that a slight increase in the performances of reconstruction approaches can be obtained when sufficient sampling densities are available; based on these results, some recommendations for operational atmospheric phase estimation procedures are proposed. Key words: InSAR; Atmospheric Phase Screen; variogram; anisotropy.
1.
INTRODUCTION
Anisotropy is often considered important in modeling 2D stochastic fields such as atmospheric phase contributions in SAR interferograms. In fact, InSAR atmospheric phase screen (APS) fields published in the literature often show evident anisotropy (see e.g. [1], or works such as [2, 3, 4, 5, 6]). _____________________________________________________ Proc. ‘Fringe 2009 Workshop’, Frascati, Italy, 30 November – 4 December 2009 (ESA SP-677, March 2010)
The concept of anisotropy is usually quantified by meaning that estimates of spatial variance, such as those obtained through the variogram function, exhibit significant differences when estimated along different directions. Anisotropy can be easily accounted for in theoretical models of spatial variance. The most common and simple form of anisotropic behavior is that referred to as geometric anisotropy, which can be described as a variability of some of the model parameters with the estimation direction. Extending a given model of spatial variance to a form which includes geometric anisotropy consists basically in adding a pair of parameters to the isotropic model, controlling the amount and direction of the anisotropy. Modeling of APS in SAR interferograms is important in quantitative applications, such as monitoring of surface displacements, in which the APS represents a nuisance. The most advanced of such applications are those allowing the monitoring of temporal trends of displacements through analysis of stable targets in series of SAR interferograms. In such applications, known as persistent scatterers interferometry (PSI), APS contributions are estimated by filtering over selected pixels containing stable radar targets, and then a continuous APS field is interpolated over the entire investigated area, to remove it from the data and so to allow more PS targets to be recognized. This interpolation is usually performed through advanced methods such as kriging, which benefit from the definition of analytical models for the spatial variance of the field to be interpolated (a concept known as model-based geostatistics [7]). Kriging of APS fields in PSI applications is usually performed by assuming a certain analytical stochastic model for the fields, then estimating the model parameters from the APS data through model fitting, and finally by performing the interpolation of the APS field over the entire grid using the analytical model with the fitted parameters. Using mode accurate models to describe spatial variation of APS fields seems thus desirable to improve the performances of PSI systems. In this paper we review some of the concepts involved in the anisotropic modeling of APS fields, focusing specifically on their use in PSI applications. We support our considerations by showing examples of real and simulated interferograms.
2.
kriging approach constitutes a classic solution, then longrange behavior of covariance model functions is less important than its asymptotic values near the origin [11].
VARIOGRAM ANISOTROPY
Given a regionalized variable z (r), and defined the (isotropic) variogram model V (h), as a function of the (scalar) lag h through V (h) = E (z(r) − z(r + h))2 , the corresponding anisotropic version can be constructed by letting h become a vector value h, and allowing for a generic elliptical deformation and rotation of the axes along which h is measured, according to the following formulae: V ′ (h) = V S · R · h (1)
For these reasons, in practical applications APS data is often modeled through a number of “well-behaved” model functions - such as the Mat´ern [12], Gauss (a special case of the Mat´ern model) or Bessel functions [13]. The Mat´ern model, with the variogram and the corresponding correlation function, can be written as follows: √ α 21−α 2 αh V (h) = N + S · 1 − · Γ(α) L √ 2 αh ; · Kα L C(h) = C(0) − V (h);
where
√ δ S= 0
cos θ 0√ ;R = − sin θ 1/ δ
sin θ . cos θ
(2)
The parameters δ and θ quantify the anisotropy of the field, i.e., as appears from (2), model the affine transformation –an axis deformation followed by a rotation– necessary to bring the V ′ variogram to the isotropic form V. Other forms of anisotropic variogram models can be conceived with more general forms of parameter variability with estimation direction, such as, for instance, multiplicative or additive models, i.e. products or sums of onedimensional models valid along orthogonal directions. This behaviour is often indicated as zonal anisotropy, to distinguish it from the geometric anisotropy paradigm stated above. However, using such models in practical applications requires particular care, since they may violate basic assumptions such as (second-order) stationarity [8] or pose problems in the kriging process itself [9]. In this work, we restrict to models extended through geometric anisotropy as in (1) and (2).
3.
MODELING ATMOSPHERIC SCREEN INSAR SIGNALS
PHASE
Hanssen [1] describes thoroughly the various components of the atmospheric phase contributions which can be expected in a SAR interferogram. Among these, the most important for the study of ground movements is the component due to spatial inhomogeneities in the water vapour content of the air. Hanssen’s review models the APS water vapour statistics through multi-fractal functions, relying on the Kolmogorov turbulence theory. The theory predicts a power-law behaviour for the spatial variation of atmospheric components, with different exponents for different scales. The multi-fractal paradigm has theoretical appeal, and is receiving attention as a viable tool for APS modeling [10], also in view of the fact that it solves some problems connected with pure power-law functions especially concerning finite asymptotic values at large lags [8]. On the other hand, if we focus on the reconstruction (i.e. interpolation from scattered data) problem, for which the
(3a)
(3b)
where Kα is the modified Bessel function of the second kind of order α. N , S, L and α are the model parameters. The Gaussian model is obtained for α = 2, and is thus based on only three parameters, N , S, and L. Another useful model is the Bessel: −α h V (h) = N + S · 1 − 2α Γ (α + 1) L (4) h · Jα , L where Jα is the Bessel function of the first kind of order α. The above models can be extended to the anisotropic form using (1) and (2), As mentioned, in the anisotropic versions two additional parameters, δ and θ, are needed to specify completely the statistics of the field.
4.
RECONSTRUCTION OF 2D APS FIELDS FROM SPARSE POINT SAMPLES
As mentioned in the introduction, real APS fields often exhibit evidently anisotropic characteristics, observation which has led to consider the use of anisotropic structure functions as more suitable descriptive tools than isotropic models [13]. As an example, we show in fig. 1(a) and (e) two ERS-1/2 tandem differential interferograms. Topographic contributions have been removed using DEM data from the SRTM mission. The SAR data are acquired over a flat region in southern Italy, thus residual contributions from different stratification of atmospheric layers, which could lead to components proportional to topography, are likely to be negligible. The short revisit time (1 day) also ensures that displacements between the two acquisitions can be ruled out. It can therefore be assumed that the whole interferometric signal is that relative to tropospheric inhomogeneities. Data have been multilooked with a factor of 2×10 (range×azimuth), and the differential phase has been unwrapped through a state-of-theart minimum cost flow algorithm (as implemented in the
[−90° − −45°) [−45° − 0°) [0° − 45°) [45° − 90°) Omindirectional
14 20
12
700
6
40 400 20 200
0 −400
4
5
800
8
variance (rad2)
10
600
10
2
Azimuth
500
abs. phase (rad)
15
400
Differences (z units)
300
variance (rad2)
200
bidimensional variogram
bidimensional variogram
Variogram 100
40 400 20 200
0 −400
0 0
0 200
400
600
800
50
1000
100
(a)
150 200 250 300 350 Spatial distance (samples)
400
450 x lag
(b)
0
−200 200
Range
0 −200
−200
2
900 1000
400
−400 x lag
(c)
y lag
400
−400
(d)
bidimensional variogram
Variogram
−200 200
y lag
bidimensional variogram
100
5
5
−5 600 700
−10
4 3 [−90° − −45°) [−45° − 0°) [0° − 45°) [45° − 90°) Omnidirectional
2
800
1 900
−15 1000
200
400
600
800
Range
(e)
1000
100
20 400 10 200
0 −400
0 −200
400 10 200
0 −400
0 0
−200 200 x lag
(f)
20
−200 0
200 300 400 Spatial distance (samples)
variance (rad2)
Azimuth
500
variance (rad2)
0 400
abs. phase (rad)
300
Differences (z2 units)
200
y lag
400
−200 200
−400
(g)
x lag
y lag
400
−400
(h)
Figure 1. ERS-1/2 unwrapped tandem differential interferograms and corresponding variograms estimated on different numbers of points. (a);(e) interferograms acquired on July 23-24, 1995, and March 24-25, 1996, respectively; color is proportional to absolute phase; multilook factor is 2×10 (range×azimuth). (b);(f) directional variograms of the interferogram in (a);(e), respectively; (c)-(d) 2-dimensional variogram of the interferogram in (a), estimated over pixels with coherence above 0.85 and 0.87, respectively. (g)-(h) same as (c)-(d), but for the interferogram in (e), with thresholds of 0.9 and 0.93, respectively. SNAPHU software [14]). In fig. 1, directional variograms, estimated by restricting the classical scalar estimator to discrete direction intervals, as well as omnidirectional variograms, are shown in panels (b) and (f), referring to the interferograms in (a) and (e), respectively. Two-dimensional variograms are also shown: panels (c)-(d) refer to the interferogram in (a), panels (g)-(h) to that in (e). The 2D variograms have been obtained from the two-dimensional version of the classical estimator: h i 2V (h) =E (z(r + h) − z(r))2 P 2 (5) ′ (r′ −r)∈Nh [z(r ) − z(r)] , = |Nh | where Nh = N(hx ,hy ) = {((x, y), (x′ , y ′ )) |(x′ − x) ≃ hx , (y ′ − y) ≃ hy }, the ≃ sign meaning equalilty within a certain tolerance interval, i.e. in practice distance vectors are binned in discrete intervals. |Nh | is the number of elements in Nh . The variograms in fig. 1 have been estimated from a subset of the pixels in the raster interferogram, selected according to their interferometric coherence (γ) value. In detail, panel (c) has been derived using points with coherence above 0.85, corresponding to an average density of about 6.5%, while in panel (d) points with γ > 0.87 have been used, corresponding to about 0.5% of the pixels. Analogously, the variogram in panel (g) has been estimated on points with γ > 0.9, corresponding to 11.8% of the pixels, while that in panel (h) from pixels with γ > 0.93 (∼ 4.5%). These examples show how decreasing point density can decrease variogram estimation performance. For the interferogram of fig. 1(a), which has a pronounced anisotropy, decreasing the point sampling by more than
an order of magnitude does not decrease significantly the estimation accuracy, since the general shape of the variogram is preserved, while in the case of the interferogram in fig. 1(d) the estimate with less points appears much harder to interpret, with spurious peaks which may hinder model fitting. It can obviously be expected that estimation procedures will perform worse when fitting a given variogram model obtained from less sample points. Estimating parameters for anisotropic versions of a certain model exhibits higher complexity with respect to the corresponding procedure using the isotropic model version (essentially given by the two additional model parameters). This may balance the fact that an anisotropic model could describe better the APS field than an isotropic model does, so that it is not a trivial point to compare, given the same point density, estimation procedures to fit experimental variograms to anisotropic and isotropic model versions. This issue will eventually reflect on the quality of APS fields interpolated over the entire interferogram grid, so that better parameters estimation should in principle lead to better reconstruction quality. In the following, we investigate on these points, proposing some observations about the performance of estimation and reconstruction algorithms, based on simulations.
5.
SIMULATIONS
Using simulated data allows to keep control of all the variables in complex problems. We generated series of random surfaces with assigned statistical characteristics; we performed model fitting from samples of these surfaces taken at random with decreasing density, and then
surface reconstruction based on the fitted models through kriging.
5.1. Surface generation Series of random surface were generated with assigned covariance functions. Conditional simulation of surfaces with given covariance models can be achieved through algorithmic tools such as Cholesky or eigendecomposition [15], assuming a Gaussian distribution. Other simulation methods for fractal surfaces rely on spectral filtering, based on the Wiener-Khintchine theorem. In the present work, series of surfaces of size 100×100 pixels were generated through spectral methods, as suggested in [1], with varying degrees of anisotropy; in detail, the results described in the following refer to series of isotropic surfaces, anisotropic surfaces with δ = 2 and random orientations θ uniformly distributed in the interval [−π, π), and another series of surfaces with δ = 10 and again random orientation angles. In view of the considerations in sect. 3 about model selection, a reference set of parameters has been estimated using variograms obtained from all samples of the simulated fields, and then through LMS fit. To investigate the robustness of fit methods to the number of samples available, variograms parameters have then been fitted to experimental variograms obtained from reduced sets of surface samples selected at random positions, with sampling spatial densities ranging from 10% to 0.5% of the surface pixels. Note that, for ERS sampling conditions at full resolution, a density of points of 1% corresponds to about 100 pixels per km2 . In the case of PSI applications, such a density of PS per km2 corresponds to conditions typical of urban environments [16], whereas on scarcely urbanised areas, much lower values are typically observed [17]. Higher values of PS densities, up to several percents, can be expected on urban areas from higher resolution SAR data such as in the case of TerraSAR-X [18] or COSMO/SkyMed. This is a huge increase in spatial density, considering that for a full-resolution TSX or CSK stripmap image with ground spacing at 3×3 m2 , 1% sampling corresponds to about 1000 pixels per km2 . Variogram model fitting in general can be a hard problem. Experienced geostatisticians use a series of practical procedures and tools to infer the best statistical model which fits experimental data. When faced with data exhibiting anisotropy, such empirical methods are an aid to infer suitable values for the anisotropy parameters and thus finally aid the model fitting procedure. However, within the specific field of PSI applications, APS raster field interpolation has to be performed tipically for each acquisition of a time series of several tens of images. Therefore, the emphasis is posed over reconstruction methods which can work in an automated way, with as little operator interventions as possible, to avoid excessive processing times. In this respect, automated parameter estimation methods should be viewed as the most useful. The simplest and best known of such methods is of course
least mean square (LMS) curve fitting. The method relies on the assumption that errors are Gaussian distributed around the “exact” values and independent. Weights can be used to guide the fit by increasing the importance of some points with respect to others. For the specific case of variogram fitting, typical weighting schemes proposed in the literature are e.g. those inversely proportional to the lag distance or its squared value, or proportional to the bin population [19].
5.2. Reconstruction of simulated surfaces In fig. 2, we show a pair of examples of model fit and reconstruction of simulated surfaces. It can be seen that in the first example (panels a and b, original surface simulated with δ = 10), the reconstruction at 5% sampling (panel a) is more faithful using the 2D model than the omnidirectional one, since in this case the high sampling density allows a good fit of the 2D variogram: the reconstructed surface better reproduces the anisotropic structure of the original. In lower sampling conditions (2%, panel b) the surface reconstructed through the 2D model does not differ too much from the one obtained by using the 1D model. In this case, then, using the more complicated 2D model seems unjustified. In the second example (panels c and d, δ = 2), which refers to a slightly anisotropic original surface, with a 10% sampling rate (panel c) the reconstruction is good with both 1D and 2D models, which give closely comparable results. If sampling is reduced to 3%, instead (panel d), the low number of available points causes the 2D model to give unreliable results, especially for the anisotropy δ and θ parameters, which results in a reconstruction which is much worse than that obtained through the 1D model. In this case, then, robustness is more important than model accuracy for the reconstruction purpose. The degradation of model fit parameters as sampling density gets lower can be appreciated by looking at the values in tab. 1. Here, the parameters for the surface analysed in fig. 2(a) and (b) are shown, obtained through WLMS fit starting from isotropic conditions, as sampling rate decreases from 100% to 0.5%. It can be noticed the decrease in accuracy, especially for the anisotropy parameters, and the corresponding increase in RMSE, as sampling rate decreases. This behaviour can be understood at least partially by making the following observations. 1. “Curse of dimensionality”: in practice, estimating the experimental variogram in a 2D domain requires more sampling points than using a 1D model, since the number of required bins increases exponentially when adopting higher-dimensional models. Consequently, it is expected that, for a certain sampling density, estimates of model parameters from lowerdimensional models are more robust than those obtained from higher-dimensional estimates, in spite of the fact that the surface a priori statistics could be
Experimental 2D variogram − all pts.
Kriged surface − 1D model
Original surface
5
8
100
Kriged surface − 2D model − true
20
40
60
80
100
Kriged surface − 2D model − iso
6
20
6
40
4
40
4
0
80 100
−2
20
40
60
80
100
40
2
60
20
40
60
80
y lags (samples)
1
3.5
3.5
2
3 0
2.5 2
0 x lags (samples)
4 3.5 3 0
2.5 2
1.5
1.5
1
1
100
50
Fitted 2D variogram − true
4
2.5
0 x lags (samples)
−50
4
3
50 −50
3
Fitted 2D variogram − iso
0
5 4
50 −50
50
1.5
−2
100
0 x lags (samples)
−50
0
80
6 0
2
50 −50
50
−50
8
20
2
20 30 lags (samples)
Fitted 2D variogram − all pts.
8
60
10
2
y lags (samples)
80
0
3
0
1
y lags (samples)
60
0
100
y lags (samples)
40
Exp. 1D Var. all pts. Exp. 1D Var. Fitted 1D Var. all pts. Fitted 1D Var.
0
0
20
2
1
80
100
3
2
2
80
y lags (samples)
4
60
7
4
variance
40 4
60
9 8
6
6
40
−50
4
20
20
Experimental 2D variogram
−50 5
8
50 −50
50
0 x lags (samples)
1 50 −50
50
0 x lags (samples)
50
(a) Original surface
Kriged surface − 1D model
Experimental 2D variogram − all pts.
6
Experimental 2D variogram
−50
−50
8 5
60
2
2
2
80
80
100 20
40
60
80
100
0
20
40
60
80
Exp. 1D Var. all pts. Exp. 1D Var. Fitted 1D Var. all pts. Fitted 1D Var.
1
0
0
100
3
100
0
Kriged surface − 2D model − true
10
20 30 lags (samples)
40
2
0 x lags (samples)
2
−50
4.5
4.5 4
3.5
4
40 60
2
3 0
2.5 2
3 0
2.5 2
3.5 3 0
2.5 2
1.5
80
0
100
20
40
60
80
1.5
1.5
80 100
100
0
3.5
y lags (samples)
60
50
4
y lags (samples)
4
0 x lags (samples)
Fitted 2D variogram − true
−50
4
40
50 −50
50
Fitted 2D variogram − iso
−50
6
20
10
5
50 −50
50
6
20
15
0
1
Fitted 2D variogram − all pts.
Kriged surface − 2D model − iso
3
0
y lags (samples)
60
y lags (samples)
variance
4
40 4
20
4
4
y lags (samples)
6
40
25
5
6
20
20
0
1 1
50 −50
20
40
60
80
0 x lags (samples)
100
50 −50
50
0 x lags (samples)
1 50 −50
50
0 x lags (samples)
50
(b) Original surface
Kriged surface − 1D model
Experimental 2D variogram − all pts. 8
Experimental 2D variogram −50
−50
10
6
60
4
6
6
8
8
5
6
60
4 3
4
2
20
40
60
80
100
Kriged surface − 2D model − true
0
100
20
40
60
80
8
4 2
80 100
0
20
40
60
80
100
0
10
20 30 lags (samples)
40
7 6
0
5 4 3
2
2
1 50 −50
50
0 x lags (samples)
1
50 −50
50
Fitted 2D variogram − iso
−50
0 x lags (samples)
50
Fitted 2D variogram − true
−50
−50
6
12
6
6
5.5
20
10 8
40
6
60
4 3
Fitted 2D variogram − all pts.
Kriged surface − 2D model − iso
10
40
0
100
12
20
1
6
60
4
5
y lags (samples)
100
0
2
80
0
Exp. 1D Var. all pts. Exp. 1D Var. Fitted 1D Var. all pts. Fitted 1D Var.
2
80
5
y lags (samples)
40
9
10
4.5 4
0
3.5 3
20
40
60
80
4
3
5
0
4
3
1.5
50 −50
100
0 x lags (samples)
2
2
2
0
100
0
2.5
2
80
5
y lags (samples)
8
40
20
7
y lags (samples)
10
7
y lags (samples)
20
12
variance
12
50 −50
50
0 x lags (samples)
50 −50
50
0 x lags (samples)
50
(c) Kriged surface − 1D model
Experimental 2D variogram − all pts. 8
10
40
8
8
20 40
6 60
4 2
80 100
0 20
40
60
80
100
Kriged surface − 2D model − true
4
60 80
18
6
6
16
4 3
2
Exp. 1D Var. all pts. Exp. 1D Var. Fitted 1D Var. all pts. Fitted 1D Var.
2 1
5 0
4 3
0
20
40
60
80
100
12 10 8 6 4 2
1 0
10
20 30 lags (samples)
40
50 −50
50
Fitted 2D variogram − all pts.
Kriged surface − 2D model − iso
14
0
2
0 100
20
7
5
6
variance
20
Experimental 2D variogram −50
−50
7
y lags (samples)
12
y lags (samples)
Original surface
0 x lags (samples)
50 −50
50
Fitted 2D variogram − iso
−50
0 x lags (samples)
50
Fitted 2D variogram − true
−50
−50
6
4
60
2
80
5
6
40
4
60
20
40
60
80
100
3.5 3
5 4
0
3
6 5 4
0
3
2
2.5
2
2
0 100
4
0
2
80
0 100
4.5
y lags (samples)
40
8
20
y lags (samples)
6
y lags (samples)
8
20
6
5.5
20
40
60
80
100
50 −50
1.5
0 x lags (samples)
50
1
50 −50
0 x lags (samples)
50
1
50 −50
0 x lags (samples)
50
(d) Figure 2. Examples of model fit and reconstruction of simulated surfaces. Each panel shows original and reconstructed surface through kriging using 1D and 2D model fit, 1D and 2d experimental variograms and theoretical Mat´ern models, according to figure titles. (a) and (b): strongly anisotropic random surface (δ = 10) sampled at 5% and 2%, respectively; (c) and (d): slightly anisotropic random surface (δ = 2) sampled at 10% and 3%, respectively
better described by the latter.
Figs. 3, 4, and 5 reports results of some of these experiments by using the Mat´ern, the Bessel, and the Gauss model, respectively. From the figures it can be seen that, in spite of the fact that the 2D variogram models in principle represent better the simulated data, the RMSE values of the reconstructed surfaces stay slightly lower when using 1D models than by adopting anisotropic models, as long as low anisotropy surfaces are considered (δ = 1 or δ = 2). When considering stronger anisotropies, a slight advantage can be seen in using the 2D models in some cases, although the error bars in all the plots overlap. Note also the relatively low sensitivity of the reconstruction performances to the particular model adopted (all the three figures show similar treds). In practice, then, stability of parameters estimation, which is higher for simpler 1D models, seems to overcome the advantage given by more flexible 2D models, which could in principle describe the data more accurately.
3 ’True’ anisotropic fit Anisotropic fit Isotropic fit
2 1 0 0
2
4
6
8
10
8
10
8
10
Surface RMSE (rad)
Anisotropic surface (1) 3 2 1 0 0
2
4
6
Surface RMSE (rad)
Anisotropic surface (2) 3 2 1 0 0
2
4 6 Sampling density (%)
Figure 3. RMSE error in reconstruction of simulated surfaces vs. density of sampled points. Plots refer to isotropic (top), anisotropic (middle), and strongly anisotropic (bottom) surfaces. Reconstruction through kriging via WLMS-fitted anisotropic Mat´ern model covariance function obtained starting from the “true” guess values (blue curves), isotropic guess values (red curves), and completely isotropic (1D) Mat´ern model (green curves).
Isotropic surface 3
Surface RMSE (rad)
The preceding observations explain what is observed by comparing original and reconstructed surfaces taking means of the root-mean-square reconstruction error (RMSE).
Surface RMSE (rad)
Isotropic surface
2. By construction, the anisotropy parameters δ and θ are shape parameters, which are better estimated “away from the origin”, i.e. from variogram samples with large distance values; the rest of the model parameters, instead, which are the same involved in the 1D form, have more influence close to the origin [11]. Consequently, weighting schemes suited for estimating anisotropy parameters, assigning significant weights to samples at large distances, are different from those suited to estimate the rest of the parameters, which should give more emphasis to the estimates close to the origin. This dichotomy can lead to errors in estimates of either type of parameters, regardless of the analytical model adopted.
’True’ anisotropic fit Anisotropic fit Isotropic fit
2 1 0 0
2
4
6
8
10
8
10
8
10
Anisotropic surface (1)
2 1 0 0
2
4
6
Anisotropic surface (2) 3
Surface RMSE (rad)
The above-mentioned results stem mainly from characteristics inherent to the (weighted) LMS fit procedure, which is affected by bin population, number of bin samples, and model complexity. To overcome at least some of the problems connected with LMS fit, more sophisticated tools could be used, although at the cost of increased computational complexity and time. For example, methods such as the restricted maximum likelihood (REML) [20] allow to obtain model parameter estimates without explicitly recurring to binning of the data distances.
Surface RMSE (rad)
3
2 1 0 0
2
4
6
Sampling density (%)
6.
CONCLUSIONS AND RECOMMENDATIONS
In this work, we have reported some observations related to the use of anisotropic models for atmospheric phase
Figure 4. RMSE error in reconstruction of simulated surfaces with Bessel model covariance function. Colors and plots as in fig. 3.
screen modeling and reconstruction, with reference to their use in persistent scatterers interferometry applications.
Surface RMSE (rad)
Isotropic surface 3 ’True’ anisotropic fit Anisotropic fit Isotropic fit
2 1 0 0
2
4
6
8
10
Surface RMSE (rad)
Anisotropic surface (1) 3 2 1 0 0
2
4
6
8
10
8
10
Surface RMSE (rad)
Anisotropic surface (2) 3 2 1 0 0
2
4 6 Sampling density (%)
Figure 5. RMSE error in reconstruction of simulated surfaces with Gauss model covariance function. Colors and plots as in fig. 3.
We have shown that two-dimensional variogram models can better represent anisotropic atmospheric phase screen fields, provided sufficient densities of sampling points are available. This is hardly the case with medium resolution sensors such as ERS or ENVISAT, while it may become a viable solution with high-resolution sensors such as TerraSAR-X or COSMO/SkyMed which promise much higher stable target densities. However, when dealing with sampling densities up to a few percents of the available pixels, LMS fit procedures can perform poorly when used with complex 2D models, due both to the increased size of the parameters state space, which would require higher numbers of experimental samples, and the increased computational complexity connected with estimating parameters of a 2D model, especially in view of the different importance of samples close and away to the origin in the estimation of parameters connected to the model 1D shape and its anisotropy degree. Consequently, when using the LMS-fitted models to reconstruct the entire raster APS field from sparse samples, the potential advantage given by the higher flexibility of 2D models can be reduced by the increased complexity just described, and cause results to be substantially similar to those obtained by using simpler one-dimensional variogram models. More sophisticated estimation methods than conventional LMS fit, e.g. REML, could perhaps be used to ensure that the use of anisotropic models give a tangible advantage with respect to 1D models. Future work will concentrate on implementing some of the above-mentioned solutions within PSI processing chains, in order to improve the APS estimation and removal step, also by better constraining expectations about phase noise reduction.
Table 1. Estimated Mat´ern parameters and surface RMSE at various sampling densities, for one simulated surface. % 100% 10% 8% 5% 3% 2% 1% 0.5%
N 0.75 0.49 0.53 0.63 0.22 0.39 0.078 10−6
S 3.54 3.87 3.93 3.69 4.21 4.28 3.54 3.90
L 8.34 8.05 8.63 8.25 7.68 8.73 7.76 4.21
α 2.59 2.64 2.27 10.00 2.97 3.58 10.00 10.00
δ 2.97 3.38 2.68 3.55 3.35 1.66 1.47 2.67
θ 2.67 −0.49 −0.50 2.75 −0.41 −0.31 −0.64 −3.14
RMSE 10−15 0.85 0.96 1.08 1.32 1.53 1.92 1.99
ACKNOWLEDGEMENTS ERS SAR data were provided by ESA through CAT-1 Project n. 5367, entitled “Subsidence monitoring in Daunia and Capitanata (Puglia Region, Italy) through multitemporal point-target DInSAR techniques”.
REFERENCES 1. R. F. Hanssen, Radar Interferometry: Data Interpretation and Error Analysis. Kluwer Academic Publishers, Dordrecht, The Netherlands, 2001. 2. A. Ferretti, F. Novali, E. Passera, and F. Rocca, “Statistical analysis of atmospheric components in ERS
SAR data,” in Proceedings of FRINGE 2005, ESAESRIN, Frascati, Italy, 28 Nov. - 2 Dec. 2005, (Online): http://earth.esa.int/workshops/fringe2005/. 3. F. Onn and H. A. Zebker, “Correction for interferometric synthetic aperture radar atmospheric phase artifacts using time series of zenith wet delay observations from a gps network,” Journal of Geophysical Research, vol. 111, 2006. 4. B. Puyss´egur, R. Michel, and J.-P. Avouac, “Tropospheric phase delay in interferometric synthetic aperture radar estimated from meteorological model and multispectral imagery,” Journal of Geophysical Research, vol. 112, 2007. 5. Z. W. Li, X. L. Ding, C. Huang, Z. R. Zou, and Y. L. Chen, “Atmospheric effects on repeat-pass InSAR measurements over Shanghai region,” Journal of Atmospheric and Solar-Terrestrial Physics, vol. 69, pp. 1344–56, 2007. 6. D. Perissin and C. Prati, “Comparison between SAR atmospheric phase screens at 30’ by means of ERS and ENVISAT data,” in Proceedings of IGARSS 2008, Boston MA, USA, 6–11 July 2008. 7. P. J. Diggle and P. J. Ribeiro Jr., “Model-based geostatistics,” Associac¸a˜ o Brasileira de Estat´ıstica, Course Notes presented at 14th SINAPE - Simp´osio Nacional de Probabilidade e Estat´ıstica, 2000, (Online: www.maths.lancs.ac.uk/∼diggle). 8. A. G. Journel and C. Huijbregts, Mining Geostatistics. London, UK: Academic Press, 1978. 9. D. E. Myers and A. Journel, “Variograms with zonal anisotropies and noninvertible kriging systems,” Mathematical Geology, vol. 22, no. 7, pp. 779–785, 1990. 10. J. P. Merryman Boncori and J. J. Mohr, “A tunable closed-form model for the structure function of tropospheric delay,” IEEE Geoscience and Remote Sensing Letters, vol. 5, no. 2, pp. 222–226, Apr. 2008. 11. M. L. Stein, Interpolation of Spatial Data: Some Theory for Kriging. New York: Springer, 1999. 12. R. Hanssen, A. Ferretti, M. Bianchi, R. Grebenitcharsky, F. Kleijer, and A. Elawar, “Atmospheric phase screen (APS) estimation and modeling for radar interferometry,” in Proceedings of FRINGE 2005, ESA-ESRIN, Frascati, Italy, 28 Nov. - 2 Dec. 2005, (Online): http://earth.esa.int/workshops/fringe2005/. 13. S. Knospe and S. J´onsson, “Covariance estimation, geostatistical prediction and simulation for insar observations in presence of strong atmospheric anisotropy,” in Proceedings of Fringe 2007 Workshop. 26-30 November, Frascati, Italy: ESA SP-649, ESA Publications Division, ESTEC, Noordwijk, The Netherland, 2007.
14. C. W. Chen and H. A. Zebker, “Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization,” Journal of the Optical Society of America A, vol. 18, no. 2, pp. 338–351, Feb. 2001. 15. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in C. Cambridge: Cambridge University Press, 1988. 16. A. Ferretti, C. Prati, and F. Rocca, “Permanent scatterers in SAR interferometry,” vol. 39, no. 1, pp. 8– 20, Jan. 2001. 17. F. Bovenga, A. Refice, R. Nutricato, L. Guerriero, and M. T. Chiaradia, “SPINUA: a flexible processing chain for ERS/ENVISAT long term interferometry,” in Proceedings of the 2004 Envisat & ERS Symposium, Salzburg, Austria, 2004. 18. D. O. Nitti, R. Nutricato, F. Bovenga, A. Refice, M. T. Chiaradia, and L. Guerriero, “Terrasar-x insar multi-pass analysis on venice (italy),” in Proceedins of SPIE Europe Remote Sensing, Berlin, Germany, Aug. 31–Sept. 3, 2009. 19. N. A. C. Cressie, Statistics for Spatial Data (revised edition). New York: John Wiley and Sons, Inc., 1993. 20. P. K. Kitanidis, “Parametric estimation of covariances of regionalized variables,” Water Resources Bulletin, vol. 23, pp. 557–567, 1987.
ABSTRACT When studying geophysical processes through Differential SAR Interferometry (DInSAR), it is often necessary to estimate and subtract signals due to atmospheric inhomogeneities. To this end, stochastic models are often used to describe atmospheric phase delays in DInSAR data. As a first approximation, these can be modelled as isotropic, though this is a simplification, because SAR interferograms often exhibit anisotropic atmospheric signals. In view of this, it is increasingly advocated the use of anisotropic models for atmospheric phase estimation. However, anisotropic models lead to increased computational complexity in estimating the correlation function parameters with respect to the isotropic case. Moreover, performances can degrade when use is made of interferograms with only a few sparse points usable for computations, such as in the case of persistent scatterers interferometry (PSI) applications, especially when this estimation has to be done in an automated way on several tens of interferograms. In the present work we critically analyse the main aspects connected with the use of anisotropic models for DInSAR atmospheric delays, and we evaluate the advantage given by anisotropic modeling of atmospheric phase in the case of sparse grids of points. Through analysis of APS simulated data, we observe that a slight increase in the performances of reconstruction approaches can be obtained when sufficient sampling densities are available; based on these results, some recommendations for operational atmospheric phase estimation procedures are proposed. Key words: InSAR; Atmospheric Phase Screen; variogram; anisotropy.
1.
INTRODUCTION
Anisotropy is often considered important in modeling 2D stochastic fields such as atmospheric phase contributions in SAR interferograms. In fact, InSAR atmospheric phase screen (APS) fields published in the literature often show evident anisotropy (see e.g. [1], or works such as [2, 3, 4, 5, 6]). _____________________________________________________ Proc. ‘Fringe 2009 Workshop’, Frascati, Italy, 30 November – 4 December 2009 (ESA SP-677, March 2010)
The concept of anisotropy is usually quantified by meaning that estimates of spatial variance, such as those obtained through the variogram function, exhibit significant differences when estimated along different directions. Anisotropy can be easily accounted for in theoretical models of spatial variance. The most common and simple form of anisotropic behavior is that referred to as geometric anisotropy, which can be described as a variability of some of the model parameters with the estimation direction. Extending a given model of spatial variance to a form which includes geometric anisotropy consists basically in adding a pair of parameters to the isotropic model, controlling the amount and direction of the anisotropy. Modeling of APS in SAR interferograms is important in quantitative applications, such as monitoring of surface displacements, in which the APS represents a nuisance. The most advanced of such applications are those allowing the monitoring of temporal trends of displacements through analysis of stable targets in series of SAR interferograms. In such applications, known as persistent scatterers interferometry (PSI), APS contributions are estimated by filtering over selected pixels containing stable radar targets, and then a continuous APS field is interpolated over the entire investigated area, to remove it from the data and so to allow more PS targets to be recognized. This interpolation is usually performed through advanced methods such as kriging, which benefit from the definition of analytical models for the spatial variance of the field to be interpolated (a concept known as model-based geostatistics [7]). Kriging of APS fields in PSI applications is usually performed by assuming a certain analytical stochastic model for the fields, then estimating the model parameters from the APS data through model fitting, and finally by performing the interpolation of the APS field over the entire grid using the analytical model with the fitted parameters. Using mode accurate models to describe spatial variation of APS fields seems thus desirable to improve the performances of PSI systems. In this paper we review some of the concepts involved in the anisotropic modeling of APS fields, focusing specifically on their use in PSI applications. We support our considerations by showing examples of real and simulated interferograms.
2.
kriging approach constitutes a classic solution, then longrange behavior of covariance model functions is less important than its asymptotic values near the origin [11].
VARIOGRAM ANISOTROPY
Given a regionalized variable z (r), and defined the (isotropic) variogram model V (h), as a function of the (scalar) lag h through V (h) = E (z(r) − z(r + h))2 , the corresponding anisotropic version can be constructed by letting h become a vector value h, and allowing for a generic elliptical deformation and rotation of the axes along which h is measured, according to the following formulae: V ′ (h) = V S · R · h (1)
For these reasons, in practical applications APS data is often modeled through a number of “well-behaved” model functions - such as the Mat´ern [12], Gauss (a special case of the Mat´ern model) or Bessel functions [13]. The Mat´ern model, with the variogram and the corresponding correlation function, can be written as follows: √ α 21−α 2 αh V (h) = N + S · 1 − · Γ(α) L √ 2 αh ; · Kα L C(h) = C(0) − V (h);
where
√ δ S= 0
cos θ 0√ ;R = − sin θ 1/ δ
sin θ . cos θ
(2)
The parameters δ and θ quantify the anisotropy of the field, i.e., as appears from (2), model the affine transformation –an axis deformation followed by a rotation– necessary to bring the V ′ variogram to the isotropic form V. Other forms of anisotropic variogram models can be conceived with more general forms of parameter variability with estimation direction, such as, for instance, multiplicative or additive models, i.e. products or sums of onedimensional models valid along orthogonal directions. This behaviour is often indicated as zonal anisotropy, to distinguish it from the geometric anisotropy paradigm stated above. However, using such models in practical applications requires particular care, since they may violate basic assumptions such as (second-order) stationarity [8] or pose problems in the kriging process itself [9]. In this work, we restrict to models extended through geometric anisotropy as in (1) and (2).
3.
MODELING ATMOSPHERIC SCREEN INSAR SIGNALS
PHASE
Hanssen [1] describes thoroughly the various components of the atmospheric phase contributions which can be expected in a SAR interferogram. Among these, the most important for the study of ground movements is the component due to spatial inhomogeneities in the water vapour content of the air. Hanssen’s review models the APS water vapour statistics through multi-fractal functions, relying on the Kolmogorov turbulence theory. The theory predicts a power-law behaviour for the spatial variation of atmospheric components, with different exponents for different scales. The multi-fractal paradigm has theoretical appeal, and is receiving attention as a viable tool for APS modeling [10], also in view of the fact that it solves some problems connected with pure power-law functions especially concerning finite asymptotic values at large lags [8]. On the other hand, if we focus on the reconstruction (i.e. interpolation from scattered data) problem, for which the
(3a)
(3b)
where Kα is the modified Bessel function of the second kind of order α. N , S, L and α are the model parameters. The Gaussian model is obtained for α = 2, and is thus based on only three parameters, N , S, and L. Another useful model is the Bessel: −α h V (h) = N + S · 1 − 2α Γ (α + 1) L (4) h · Jα , L where Jα is the Bessel function of the first kind of order α. The above models can be extended to the anisotropic form using (1) and (2), As mentioned, in the anisotropic versions two additional parameters, δ and θ, are needed to specify completely the statistics of the field.
4.
RECONSTRUCTION OF 2D APS FIELDS FROM SPARSE POINT SAMPLES
As mentioned in the introduction, real APS fields often exhibit evidently anisotropic characteristics, observation which has led to consider the use of anisotropic structure functions as more suitable descriptive tools than isotropic models [13]. As an example, we show in fig. 1(a) and (e) two ERS-1/2 tandem differential interferograms. Topographic contributions have been removed using DEM data from the SRTM mission. The SAR data are acquired over a flat region in southern Italy, thus residual contributions from different stratification of atmospheric layers, which could lead to components proportional to topography, are likely to be negligible. The short revisit time (1 day) also ensures that displacements between the two acquisitions can be ruled out. It can therefore be assumed that the whole interferometric signal is that relative to tropospheric inhomogeneities. Data have been multilooked with a factor of 2×10 (range×azimuth), and the differential phase has been unwrapped through a state-of-theart minimum cost flow algorithm (as implemented in the
[−90° − −45°) [−45° − 0°) [0° − 45°) [45° − 90°) Omindirectional
14 20
12
700
6
40 400 20 200
0 −400
4
5
800
8
variance (rad2)
10
600
10
2
Azimuth
500
abs. phase (rad)
15
400
Differences (z units)
300
variance (rad2)
200
bidimensional variogram
bidimensional variogram
Variogram 100
40 400 20 200
0 −400
0 0
0 200
400
600
800
50
1000
100
(a)
150 200 250 300 350 Spatial distance (samples)
400
450 x lag
(b)
0
−200 200
Range
0 −200
−200
2
900 1000
400
−400 x lag
(c)
y lag
400
−400
(d)
bidimensional variogram
Variogram
−200 200
y lag
bidimensional variogram
100
5
5
−5 600 700
−10
4 3 [−90° − −45°) [−45° − 0°) [0° − 45°) [45° − 90°) Omnidirectional
2
800
1 900
−15 1000
200
400
600
800
Range
(e)
1000
100
20 400 10 200
0 −400
0 −200
400 10 200
0 −400
0 0
−200 200 x lag
(f)
20
−200 0
200 300 400 Spatial distance (samples)
variance (rad2)
Azimuth
500
variance (rad2)
0 400
abs. phase (rad)
300
Differences (z2 units)
200
y lag
400
−200 200
−400
(g)
x lag
y lag
400
−400
(h)
Figure 1. ERS-1/2 unwrapped tandem differential interferograms and corresponding variograms estimated on different numbers of points. (a);(e) interferograms acquired on July 23-24, 1995, and March 24-25, 1996, respectively; color is proportional to absolute phase; multilook factor is 2×10 (range×azimuth). (b);(f) directional variograms of the interferogram in (a);(e), respectively; (c)-(d) 2-dimensional variogram of the interferogram in (a), estimated over pixels with coherence above 0.85 and 0.87, respectively. (g)-(h) same as (c)-(d), but for the interferogram in (e), with thresholds of 0.9 and 0.93, respectively. SNAPHU software [14]). In fig. 1, directional variograms, estimated by restricting the classical scalar estimator to discrete direction intervals, as well as omnidirectional variograms, are shown in panels (b) and (f), referring to the interferograms in (a) and (e), respectively. Two-dimensional variograms are also shown: panels (c)-(d) refer to the interferogram in (a), panels (g)-(h) to that in (e). The 2D variograms have been obtained from the two-dimensional version of the classical estimator: h i 2V (h) =E (z(r + h) − z(r))2 P 2 (5) ′ (r′ −r)∈Nh [z(r ) − z(r)] , = |Nh | where Nh = N(hx ,hy ) = {((x, y), (x′ , y ′ )) |(x′ − x) ≃ hx , (y ′ − y) ≃ hy }, the ≃ sign meaning equalilty within a certain tolerance interval, i.e. in practice distance vectors are binned in discrete intervals. |Nh | is the number of elements in Nh . The variograms in fig. 1 have been estimated from a subset of the pixels in the raster interferogram, selected according to their interferometric coherence (γ) value. In detail, panel (c) has been derived using points with coherence above 0.85, corresponding to an average density of about 6.5%, while in panel (d) points with γ > 0.87 have been used, corresponding to about 0.5% of the pixels. Analogously, the variogram in panel (g) has been estimated on points with γ > 0.9, corresponding to 11.8% of the pixels, while that in panel (h) from pixels with γ > 0.93 (∼ 4.5%). These examples show how decreasing point density can decrease variogram estimation performance. For the interferogram of fig. 1(a), which has a pronounced anisotropy, decreasing the point sampling by more than
an order of magnitude does not decrease significantly the estimation accuracy, since the general shape of the variogram is preserved, while in the case of the interferogram in fig. 1(d) the estimate with less points appears much harder to interpret, with spurious peaks which may hinder model fitting. It can obviously be expected that estimation procedures will perform worse when fitting a given variogram model obtained from less sample points. Estimating parameters for anisotropic versions of a certain model exhibits higher complexity with respect to the corresponding procedure using the isotropic model version (essentially given by the two additional model parameters). This may balance the fact that an anisotropic model could describe better the APS field than an isotropic model does, so that it is not a trivial point to compare, given the same point density, estimation procedures to fit experimental variograms to anisotropic and isotropic model versions. This issue will eventually reflect on the quality of APS fields interpolated over the entire interferogram grid, so that better parameters estimation should in principle lead to better reconstruction quality. In the following, we investigate on these points, proposing some observations about the performance of estimation and reconstruction algorithms, based on simulations.
5.
SIMULATIONS
Using simulated data allows to keep control of all the variables in complex problems. We generated series of random surfaces with assigned statistical characteristics; we performed model fitting from samples of these surfaces taken at random with decreasing density, and then
surface reconstruction based on the fitted models through kriging.
5.1. Surface generation Series of random surface were generated with assigned covariance functions. Conditional simulation of surfaces with given covariance models can be achieved through algorithmic tools such as Cholesky or eigendecomposition [15], assuming a Gaussian distribution. Other simulation methods for fractal surfaces rely on spectral filtering, based on the Wiener-Khintchine theorem. In the present work, series of surfaces of size 100×100 pixels were generated through spectral methods, as suggested in [1], with varying degrees of anisotropy; in detail, the results described in the following refer to series of isotropic surfaces, anisotropic surfaces with δ = 2 and random orientations θ uniformly distributed in the interval [−π, π), and another series of surfaces with δ = 10 and again random orientation angles. In view of the considerations in sect. 3 about model selection, a reference set of parameters has been estimated using variograms obtained from all samples of the simulated fields, and then through LMS fit. To investigate the robustness of fit methods to the number of samples available, variograms parameters have then been fitted to experimental variograms obtained from reduced sets of surface samples selected at random positions, with sampling spatial densities ranging from 10% to 0.5% of the surface pixels. Note that, for ERS sampling conditions at full resolution, a density of points of 1% corresponds to about 100 pixels per km2 . In the case of PSI applications, such a density of PS per km2 corresponds to conditions typical of urban environments [16], whereas on scarcely urbanised areas, much lower values are typically observed [17]. Higher values of PS densities, up to several percents, can be expected on urban areas from higher resolution SAR data such as in the case of TerraSAR-X [18] or COSMO/SkyMed. This is a huge increase in spatial density, considering that for a full-resolution TSX or CSK stripmap image with ground spacing at 3×3 m2 , 1% sampling corresponds to about 1000 pixels per km2 . Variogram model fitting in general can be a hard problem. Experienced geostatisticians use a series of practical procedures and tools to infer the best statistical model which fits experimental data. When faced with data exhibiting anisotropy, such empirical methods are an aid to infer suitable values for the anisotropy parameters and thus finally aid the model fitting procedure. However, within the specific field of PSI applications, APS raster field interpolation has to be performed tipically for each acquisition of a time series of several tens of images. Therefore, the emphasis is posed over reconstruction methods which can work in an automated way, with as little operator interventions as possible, to avoid excessive processing times. In this respect, automated parameter estimation methods should be viewed as the most useful. The simplest and best known of such methods is of course
least mean square (LMS) curve fitting. The method relies on the assumption that errors are Gaussian distributed around the “exact” values and independent. Weights can be used to guide the fit by increasing the importance of some points with respect to others. For the specific case of variogram fitting, typical weighting schemes proposed in the literature are e.g. those inversely proportional to the lag distance or its squared value, or proportional to the bin population [19].
5.2. Reconstruction of simulated surfaces In fig. 2, we show a pair of examples of model fit and reconstruction of simulated surfaces. It can be seen that in the first example (panels a and b, original surface simulated with δ = 10), the reconstruction at 5% sampling (panel a) is more faithful using the 2D model than the omnidirectional one, since in this case the high sampling density allows a good fit of the 2D variogram: the reconstructed surface better reproduces the anisotropic structure of the original. In lower sampling conditions (2%, panel b) the surface reconstructed through the 2D model does not differ too much from the one obtained by using the 1D model. In this case, then, using the more complicated 2D model seems unjustified. In the second example (panels c and d, δ = 2), which refers to a slightly anisotropic original surface, with a 10% sampling rate (panel c) the reconstruction is good with both 1D and 2D models, which give closely comparable results. If sampling is reduced to 3%, instead (panel d), the low number of available points causes the 2D model to give unreliable results, especially for the anisotropy δ and θ parameters, which results in a reconstruction which is much worse than that obtained through the 1D model. In this case, then, robustness is more important than model accuracy for the reconstruction purpose. The degradation of model fit parameters as sampling density gets lower can be appreciated by looking at the values in tab. 1. Here, the parameters for the surface analysed in fig. 2(a) and (b) are shown, obtained through WLMS fit starting from isotropic conditions, as sampling rate decreases from 100% to 0.5%. It can be noticed the decrease in accuracy, especially for the anisotropy parameters, and the corresponding increase in RMSE, as sampling rate decreases. This behaviour can be understood at least partially by making the following observations. 1. “Curse of dimensionality”: in practice, estimating the experimental variogram in a 2D domain requires more sampling points than using a 1D model, since the number of required bins increases exponentially when adopting higher-dimensional models. Consequently, it is expected that, for a certain sampling density, estimates of model parameters from lowerdimensional models are more robust than those obtained from higher-dimensional estimates, in spite of the fact that the surface a priori statistics could be
Experimental 2D variogram − all pts.
Kriged surface − 1D model
Original surface
5
8
100
Kriged surface − 2D model − true
20
40
60
80
100
Kriged surface − 2D model − iso
6
20
6
40
4
40
4
0
80 100
−2
20
40
60
80
100
40
2
60
20
40
60
80
y lags (samples)
1
3.5
3.5
2
3 0
2.5 2
0 x lags (samples)
4 3.5 3 0
2.5 2
1.5
1.5
1
1
100
50
Fitted 2D variogram − true
4
2.5
0 x lags (samples)
−50
4
3
50 −50
3
Fitted 2D variogram − iso
0
5 4
50 −50
50
1.5
−2
100
0 x lags (samples)
−50
0
80
6 0
2
50 −50
50
−50
8
20
2
20 30 lags (samples)
Fitted 2D variogram − all pts.
8
60
10
2
y lags (samples)
80
0
3
0
1
y lags (samples)
60
0
100
y lags (samples)
40
Exp. 1D Var. all pts. Exp. 1D Var. Fitted 1D Var. all pts. Fitted 1D Var.
0
0
20
2
1
80
100
3
2
2
80
y lags (samples)
4
60
7
4
variance
40 4
60
9 8
6
6
40
−50
4
20
20
Experimental 2D variogram
−50 5
8
50 −50
50
0 x lags (samples)
1 50 −50
50
0 x lags (samples)
50
(a) Original surface
Kriged surface − 1D model
Experimental 2D variogram − all pts.
6
Experimental 2D variogram
−50
−50
8 5
60
2
2
2
80
80
100 20
40
60
80
100
0
20
40
60
80
Exp. 1D Var. all pts. Exp. 1D Var. Fitted 1D Var. all pts. Fitted 1D Var.
1
0
0
100
3
100
0
Kriged surface − 2D model − true
10
20 30 lags (samples)
40
2
0 x lags (samples)
2
−50
4.5
4.5 4
3.5
4
40 60
2
3 0
2.5 2
3 0
2.5 2
3.5 3 0
2.5 2
1.5
80
0
100
20
40
60
80
1.5
1.5
80 100
100
0
3.5
y lags (samples)
60
50
4
y lags (samples)
4
0 x lags (samples)
Fitted 2D variogram − true
−50
4
40
50 −50
50
Fitted 2D variogram − iso
−50
6
20
10
5
50 −50
50
6
20
15
0
1
Fitted 2D variogram − all pts.
Kriged surface − 2D model − iso
3
0
y lags (samples)
60
y lags (samples)
variance
4
40 4
20
4
4
y lags (samples)
6
40
25
5
6
20
20
0
1 1
50 −50
20
40
60
80
0 x lags (samples)
100
50 −50
50
0 x lags (samples)
1 50 −50
50
0 x lags (samples)
50
(b) Original surface
Kriged surface − 1D model
Experimental 2D variogram − all pts. 8
Experimental 2D variogram −50
−50
10
6
60
4
6
6
8
8
5
6
60
4 3
4
2
20
40
60
80
100
Kriged surface − 2D model − true
0
100
20
40
60
80
8
4 2
80 100
0
20
40
60
80
100
0
10
20 30 lags (samples)
40
7 6
0
5 4 3
2
2
1 50 −50
50
0 x lags (samples)
1
50 −50
50
Fitted 2D variogram − iso
−50
0 x lags (samples)
50
Fitted 2D variogram − true
−50
−50
6
12
6
6
5.5
20
10 8
40
6
60
4 3
Fitted 2D variogram − all pts.
Kriged surface − 2D model − iso
10
40
0
100
12
20
1
6
60
4
5
y lags (samples)
100
0
2
80
0
Exp. 1D Var. all pts. Exp. 1D Var. Fitted 1D Var. all pts. Fitted 1D Var.
2
80
5
y lags (samples)
40
9
10
4.5 4
0
3.5 3
20
40
60
80
4
3
5
0
4
3
1.5
50 −50
100
0 x lags (samples)
2
2
2
0
100
0
2.5
2
80
5
y lags (samples)
8
40
20
7
y lags (samples)
10
7
y lags (samples)
20
12
variance
12
50 −50
50
0 x lags (samples)
50 −50
50
0 x lags (samples)
50
(c) Kriged surface − 1D model
Experimental 2D variogram − all pts. 8
10
40
8
8
20 40
6 60
4 2
80 100
0 20
40
60
80
100
Kriged surface − 2D model − true
4
60 80
18
6
6
16
4 3
2
Exp. 1D Var. all pts. Exp. 1D Var. Fitted 1D Var. all pts. Fitted 1D Var.
2 1
5 0
4 3
0
20
40
60
80
100
12 10 8 6 4 2
1 0
10
20 30 lags (samples)
40
50 −50
50
Fitted 2D variogram − all pts.
Kriged surface − 2D model − iso
14
0
2
0 100
20
7
5
6
variance
20
Experimental 2D variogram −50
−50
7
y lags (samples)
12
y lags (samples)
Original surface
0 x lags (samples)
50 −50
50
Fitted 2D variogram − iso
−50
0 x lags (samples)
50
Fitted 2D variogram − true
−50
−50
6
4
60
2
80
5
6
40
4
60
20
40
60
80
100
3.5 3
5 4
0
3
6 5 4
0
3
2
2.5
2
2
0 100
4
0
2
80
0 100
4.5
y lags (samples)
40
8
20
y lags (samples)
6
y lags (samples)
8
20
6
5.5
20
40
60
80
100
50 −50
1.5
0 x lags (samples)
50
1
50 −50
0 x lags (samples)
50
1
50 −50
0 x lags (samples)
50
(d) Figure 2. Examples of model fit and reconstruction of simulated surfaces. Each panel shows original and reconstructed surface through kriging using 1D and 2D model fit, 1D and 2d experimental variograms and theoretical Mat´ern models, according to figure titles. (a) and (b): strongly anisotropic random surface (δ = 10) sampled at 5% and 2%, respectively; (c) and (d): slightly anisotropic random surface (δ = 2) sampled at 10% and 3%, respectively
better described by the latter.
Figs. 3, 4, and 5 reports results of some of these experiments by using the Mat´ern, the Bessel, and the Gauss model, respectively. From the figures it can be seen that, in spite of the fact that the 2D variogram models in principle represent better the simulated data, the RMSE values of the reconstructed surfaces stay slightly lower when using 1D models than by adopting anisotropic models, as long as low anisotropy surfaces are considered (δ = 1 or δ = 2). When considering stronger anisotropies, a slight advantage can be seen in using the 2D models in some cases, although the error bars in all the plots overlap. Note also the relatively low sensitivity of the reconstruction performances to the particular model adopted (all the three figures show similar treds). In practice, then, stability of parameters estimation, which is higher for simpler 1D models, seems to overcome the advantage given by more flexible 2D models, which could in principle describe the data more accurately.
3 ’True’ anisotropic fit Anisotropic fit Isotropic fit
2 1 0 0
2
4
6
8
10
8
10
8
10
Surface RMSE (rad)
Anisotropic surface (1) 3 2 1 0 0
2
4
6
Surface RMSE (rad)
Anisotropic surface (2) 3 2 1 0 0
2
4 6 Sampling density (%)
Figure 3. RMSE error in reconstruction of simulated surfaces vs. density of sampled points. Plots refer to isotropic (top), anisotropic (middle), and strongly anisotropic (bottom) surfaces. Reconstruction through kriging via WLMS-fitted anisotropic Mat´ern model covariance function obtained starting from the “true” guess values (blue curves), isotropic guess values (red curves), and completely isotropic (1D) Mat´ern model (green curves).
Isotropic surface 3
Surface RMSE (rad)
The preceding observations explain what is observed by comparing original and reconstructed surfaces taking means of the root-mean-square reconstruction error (RMSE).
Surface RMSE (rad)
Isotropic surface
2. By construction, the anisotropy parameters δ and θ are shape parameters, which are better estimated “away from the origin”, i.e. from variogram samples with large distance values; the rest of the model parameters, instead, which are the same involved in the 1D form, have more influence close to the origin [11]. Consequently, weighting schemes suited for estimating anisotropy parameters, assigning significant weights to samples at large distances, are different from those suited to estimate the rest of the parameters, which should give more emphasis to the estimates close to the origin. This dichotomy can lead to errors in estimates of either type of parameters, regardless of the analytical model adopted.
’True’ anisotropic fit Anisotropic fit Isotropic fit
2 1 0 0
2
4
6
8
10
8
10
8
10
Anisotropic surface (1)
2 1 0 0
2
4
6
Anisotropic surface (2) 3
Surface RMSE (rad)
The above-mentioned results stem mainly from characteristics inherent to the (weighted) LMS fit procedure, which is affected by bin population, number of bin samples, and model complexity. To overcome at least some of the problems connected with LMS fit, more sophisticated tools could be used, although at the cost of increased computational complexity and time. For example, methods such as the restricted maximum likelihood (REML) [20] allow to obtain model parameter estimates without explicitly recurring to binning of the data distances.
Surface RMSE (rad)
3
2 1 0 0
2
4
6
Sampling density (%)
6.
CONCLUSIONS AND RECOMMENDATIONS
In this work, we have reported some observations related to the use of anisotropic models for atmospheric phase
Figure 4. RMSE error in reconstruction of simulated surfaces with Bessel model covariance function. Colors and plots as in fig. 3.
screen modeling and reconstruction, with reference to their use in persistent scatterers interferometry applications.
Surface RMSE (rad)
Isotropic surface 3 ’True’ anisotropic fit Anisotropic fit Isotropic fit
2 1 0 0
2
4
6
8
10
Surface RMSE (rad)
Anisotropic surface (1) 3 2 1 0 0
2
4
6
8
10
8
10
Surface RMSE (rad)
Anisotropic surface (2) 3 2 1 0 0
2
4 6 Sampling density (%)
Figure 5. RMSE error in reconstruction of simulated surfaces with Gauss model covariance function. Colors and plots as in fig. 3.
We have shown that two-dimensional variogram models can better represent anisotropic atmospheric phase screen fields, provided sufficient densities of sampling points are available. This is hardly the case with medium resolution sensors such as ERS or ENVISAT, while it may become a viable solution with high-resolution sensors such as TerraSAR-X or COSMO/SkyMed which promise much higher stable target densities. However, when dealing with sampling densities up to a few percents of the available pixels, LMS fit procedures can perform poorly when used with complex 2D models, due both to the increased size of the parameters state space, which would require higher numbers of experimental samples, and the increased computational complexity connected with estimating parameters of a 2D model, especially in view of the different importance of samples close and away to the origin in the estimation of parameters connected to the model 1D shape and its anisotropy degree. Consequently, when using the LMS-fitted models to reconstruct the entire raster APS field from sparse samples, the potential advantage given by the higher flexibility of 2D models can be reduced by the increased complexity just described, and cause results to be substantially similar to those obtained by using simpler one-dimensional variogram models. More sophisticated estimation methods than conventional LMS fit, e.g. REML, could perhaps be used to ensure that the use of anisotropic models give a tangible advantage with respect to 1D models. Future work will concentrate on implementing some of the above-mentioned solutions within PSI processing chains, in order to improve the APS estimation and removal step, also by better constraining expectations about phase noise reduction.
Table 1. Estimated Mat´ern parameters and surface RMSE at various sampling densities, for one simulated surface. % 100% 10% 8% 5% 3% 2% 1% 0.5%
N 0.75 0.49 0.53 0.63 0.22 0.39 0.078 10−6
S 3.54 3.87 3.93 3.69 4.21 4.28 3.54 3.90
L 8.34 8.05 8.63 8.25 7.68 8.73 7.76 4.21
α 2.59 2.64 2.27 10.00 2.97 3.58 10.00 10.00
δ 2.97 3.38 2.68 3.55 3.35 1.66 1.47 2.67
θ 2.67 −0.49 −0.50 2.75 −0.41 −0.31 −0.64 −3.14
RMSE 10−15 0.85 0.96 1.08 1.32 1.53 1.92 1.99
ACKNOWLEDGEMENTS ERS SAR data were provided by ESA through CAT-1 Project n. 5367, entitled “Subsidence monitoring in Daunia and Capitanata (Puglia Region, Italy) through multitemporal point-target DInSAR techniques”.
REFERENCES 1. R. F. Hanssen, Radar Interferometry: Data Interpretation and Error Analysis. Kluwer Academic Publishers, Dordrecht, The Netherlands, 2001. 2. A. Ferretti, F. Novali, E. Passera, and F. Rocca, “Statistical analysis of atmospheric components in ERS
SAR data,” in Proceedings of FRINGE 2005, ESAESRIN, Frascati, Italy, 28 Nov. - 2 Dec. 2005, (Online): http://earth.esa.int/workshops/fringe2005/. 3. F. Onn and H. A. Zebker, “Correction for interferometric synthetic aperture radar atmospheric phase artifacts using time series of zenith wet delay observations from a gps network,” Journal of Geophysical Research, vol. 111, 2006. 4. B. Puyss´egur, R. Michel, and J.-P. Avouac, “Tropospheric phase delay in interferometric synthetic aperture radar estimated from meteorological model and multispectral imagery,” Journal of Geophysical Research, vol. 112, 2007. 5. Z. W. Li, X. L. Ding, C. Huang, Z. R. Zou, and Y. L. Chen, “Atmospheric effects on repeat-pass InSAR measurements over Shanghai region,” Journal of Atmospheric and Solar-Terrestrial Physics, vol. 69, pp. 1344–56, 2007. 6. D. Perissin and C. Prati, “Comparison between SAR atmospheric phase screens at 30’ by means of ERS and ENVISAT data,” in Proceedings of IGARSS 2008, Boston MA, USA, 6–11 July 2008. 7. P. J. Diggle and P. J. Ribeiro Jr., “Model-based geostatistics,” Associac¸a˜ o Brasileira de Estat´ıstica, Course Notes presented at 14th SINAPE - Simp´osio Nacional de Probabilidade e Estat´ıstica, 2000, (Online: www.maths.lancs.ac.uk/∼diggle). 8. A. G. Journel and C. Huijbregts, Mining Geostatistics. London, UK: Academic Press, 1978. 9. D. E. Myers and A. Journel, “Variograms with zonal anisotropies and noninvertible kriging systems,” Mathematical Geology, vol. 22, no. 7, pp. 779–785, 1990. 10. J. P. Merryman Boncori and J. J. Mohr, “A tunable closed-form model for the structure function of tropospheric delay,” IEEE Geoscience and Remote Sensing Letters, vol. 5, no. 2, pp. 222–226, Apr. 2008. 11. M. L. Stein, Interpolation of Spatial Data: Some Theory for Kriging. New York: Springer, 1999. 12. R. Hanssen, A. Ferretti, M. Bianchi, R. Grebenitcharsky, F. Kleijer, and A. Elawar, “Atmospheric phase screen (APS) estimation and modeling for radar interferometry,” in Proceedings of FRINGE 2005, ESA-ESRIN, Frascati, Italy, 28 Nov. - 2 Dec. 2005, (Online): http://earth.esa.int/workshops/fringe2005/. 13. S. Knospe and S. J´onsson, “Covariance estimation, geostatistical prediction and simulation for insar observations in presence of strong atmospheric anisotropy,” in Proceedings of Fringe 2007 Workshop. 26-30 November, Frascati, Italy: ESA SP-649, ESA Publications Division, ESTEC, Noordwijk, The Netherland, 2007.
14. C. W. Chen and H. A. Zebker, “Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization,” Journal of the Optical Society of America A, vol. 18, no. 2, pp. 338–351, Feb. 2001. 15. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in C. Cambridge: Cambridge University Press, 1988. 16. A. Ferretti, C. Prati, and F. Rocca, “Permanent scatterers in SAR interferometry,” vol. 39, no. 1, pp. 8– 20, Jan. 2001. 17. F. Bovenga, A. Refice, R. Nutricato, L. Guerriero, and M. T. Chiaradia, “SPINUA: a flexible processing chain for ERS/ENVISAT long term interferometry,” in Proceedings of the 2004 Envisat & ERS Symposium, Salzburg, Austria, 2004. 18. D. O. Nitti, R. Nutricato, F. Bovenga, A. Refice, M. T. Chiaradia, and L. Guerriero, “Terrasar-x insar multi-pass analysis on venice (italy),” in Proceedins of SPIE Europe Remote Sensing, Berlin, Germany, Aug. 31–Sept. 3, 2009. 19. N. A. C. Cressie, Statistics for Spatial Data (revised edition). New York: John Wiley and Sons, Inc., 1993. 20. P. K. Kitanidis, “Parametric estimation of covariances of regionalized variables,” Water Resources Bulletin, vol. 23, pp. 557–567, 1987.
