Abdourrahmane M Atto
email: abdourrahmane.atto@univ-savoie.fr

Emmanuel Trouvé
email: emmanuel.trouve@univ-savoie.fr

Yannick Berthoumieu
email: yannick.berthoumieu@ims-bordeaux.fr

Grégoire Mercier
email: mercier@telecom-bretagne.fr

Emmanuel Trouv

Gr Egoire Mercier

Multi-Date Divergence Matrices for the Analysis of SAR Image Time Series

Keywords: Image Time Series, Wavelets, Curvelets, Parametric modeling, Kullback-Leibler Divergence, Change Detection

come

I. Introduction

I MAGE analysis from remote sensing time series plays an essential role in many applications involving monitoring of regions that are dicult to access, all over the world.

The analysis of remote sensing time series has been addressed in the literature by exploiting some specic (class dependent) image features and the signal coherency in time. For instance, the Synthetic Aperture Radar (SAR) image features proposed in [START_REF] Bruzzone | advanced system for the automatic classication of multitemporal SAR images[END_REF] consist of long-term coherence and backscattering temporal variability learned from dierent land-cover classes, for a classication purpose. For ground deformation monitoring, [START_REF] Guoxiang | \Estimating spatiotemporal ground deformation with improved permanent-scatterer radar interferometry[END_REF] has proposed estimating some linear and nonlinear deformation features from neighborhood networking of temporally coherent radar targets (permanent-scatterer) in interferometric SAR time series. For the detection of buildings destroyed in an earthquake, [START_REF] Brunner | \Earthquake damage assessment of buildings using VHR optical and SAR imagery[END_REF] has constructed some building features consisting in non-overlapping rectangular footprints and has proposed a post-event building feature predictor as a benchmark relating building states (destroyed or nondestroyed) over time.

The analysis of remote sensing time series has also been addressed by transforming the observed data in some appropriate representation spaces. As a matter of example, the detection of spatial similarities has been addressed in the spectral domain in [START_REF] Lhermitte | \Hierarchical image segmentation based on similarity of NDVI time series[END_REF] and [START_REF] Bastarrika | \Automatic burned land mapping from MODIS time series images: Assessment in mediterranean ecosystems[END_REF]. This detection has been performed upon a compressed representation of the change information in [START_REF] Bovolo | framework for automatic and unsupervised detection of multiple changes in multitemporal images[END_REF]. This detection has also been addressed by focusing on a low frequency representation of the SAR signal from some iterative ltering of the observed data to remove speckle in [START_REF] Bazi | unsupervised approach based on the generalized gaussian model to automatic change detection in multitemporal SAR images[END_REF].

In addition, the selection of suitable similarity measures with respect to the remote sensing data distribution has deserved much interest in the literature on multi-temporal change detection. In particular, we have:

• the time warping distance associated with a diskbased sux tree indexing is a accurate in comparing sequences of dierent lengths and/or dierent sampling rates, see [START_REF] Park | \Similarity-based subsequence search in image sequence databases[END_REF],

• the Bhattacharyya distance between the probability density functions is ecient for measuring contrast similarity, see [START_REF] Morio | \Information theory-based approach for contrast analysis in polarimetric and/or interferometric sar images[END_REF],

• the mutual information between Wishart processes over time is suitable for characterizing temporal polarimetric and interferometric informations in SAR data, see [START_REF] Erten | new coherent similarity measure for temporal multichannel scene characterization[END_REF]. Exploiting new generation remote sensing images is actually facing three major challenges: high resolution, large database (long acquisition sequence) and speckle eect. The conjunction of these facts makes the pursuit of a temporal event intricate on a long time series, when no a priori is available concerning i) the spatial location of changing areas and ii) the change information. Indeed, local analysis from overlapped sliding windows all over the scene and all along the whole time series is not thinkable due to computational complexity. This analysis is however necessary for the concerns of early detection of some physical phenomenon such as ground displacement or melting cycle of glaciers that can be evaluated through long image time series.

The approach proposed in this work for simplifying the analysis of image time series relates to the description In this diagram, transform T designates either the SWT or the FDCT. Specically, the term KL(T X/T Y) designates KL based similarity measurements from the SWT or the FDCT decompositions of X and Y, see for instance [START_REF] Brunner | \Earthquake damage assessment of buildings using VHR optical and SAR imagery[END_REF] and [START_REF] Lhermitte | \Hierarchical image segmentation based on similarity of NDVI time series[END_REF]. Temporal analysis involves large size image (or subimage) comparison whereas spatial analysis applies only for specic images at the pixel level. The MDDM shown in the center of the gure has been computed upon a multi-date Gaussian random eld composed with 3 stationary subsequences.

I j-1 I j I j+1 I k-1 I k I k+1 I -1 I I +1 > t • • • • • • • • • • • • • • • • • • • • • Subsequence j Subsequence k Subsequence
of any image/subimage by a set of essential description parameters obtained from representations having suitable statistical and geometrical properties. These representations are the discrete Stationary Wavelet Transform, SWT, [START_REF] Coifman | Translation invariant denoising[END_REF] and the Fast Discrete Curvelet Transform, FDCT, [START_REF] Cand | \Curvelets -a surprisingly eective nonadaptive representation for objects with edges[END_REF] [13], [START_REF] Cand | \Fast discrete curvelet transforms[END_REF].

The SWT and FDCT coecients of the SAR images exhibit stochastic nature because of the presence of speckle in radar images. In this respect, we consider parameter extraction from stochastic modeling of the coecients of these representations. The modeling is achieved by using 1) Edgeworth expansions for the SWT and FDCT approximations and 2) Generalized Gaussian, Log-Normal and Weibull distributions for the details and orientations of these transforms. These distributions are used in a joint framework requiring similarity evaluation from distribution models pertained to dierent parametric families: the closed-form expressions for inter-family parameter

Input image Im SWT representation I W of Im Approximation c J,0 κ 1 [c J,0 ] κ 2 [c J,0 ] κ 3 [c J,0 ] κ 4 [c J,0 ] Details c j,n 1 j J n∈{1,2,3}
Model selection for c j,n

ρ 1 c j,n ρ 2 c j,n FDCT representation I C of Im Approximation I app C κ 1 [I app C ] κ 2 [I app C ] κ 3 [I app C ] κ 4 [I app C ] Orientations I ang C [q] q∈{0,1,...,Q-1}
Model selection for I ang

C [q] ρ 1 I ang C [q] ρ 2 I ang C [q]
Fig. comparison are then derived in the paper.

The proposed similarity measurements are based on the Kullback-Leibler (KL) divergences and applied in a two stage approach illustrated in Fig. 1. In the rst stage, a global temporal analysis is performed at the image/subimage level, i.e between any pair of image description parameters. These KL based similarity measurements are gathered to form a Multi-Date Divergence Matrix (MDDM). The analysis of the MDDM eases the selection of the acquisitions-of-interest: these acquisitions are either change-images, i.e images relating the transitions between the stationary subsequences composing the image time series or image-outliers, i.e acquisitions that exhibit a suspicious behavior with respect to their neighborhood in time.

In the second stage, a local spatial analysis is performed at the pixel level between acquisitions-of-interest derived from the MDDM analysis. Pixelwise similarity measurements are computed within a sliding window in order to build change maps and detect the changing areas.

The presentation of the paper is as follows. Section II addresses the extraction of image description parameters and the similarity measurements in the wavelet and curvelet domains. Section III presents the MDDM and provides an interpretation of this matrix in term of abrupt and progressive changes, as well as multiple change-images in time. Section IV is dedicated to experimental results on multi-temporal change detection in SAR time series from TerraSAR-X (TS-X) images. Section V addresses spatial change detection from MDDM critical acquisitions. Section VI concludes the work.

I m = I[t m ]
where t m is the acquisition time.

The method proposed below for analyzing I relates to a parsimonious description of every I m by a set of P essential description parameters extracted from the image features. This set is required to be parsimonious in the sense that it consists of an accurate modeling with P N elements, where N is the number of pixels of I m (images are assumed to have the same sample size). In this respect, the set of parameters used hereafter derives from parametric modeling of the statistical distributions of wavelet and curvelet contributions. Indeed, in presence of non-stationarities (the modeling concerns an image/subimage and not a stationary region from this image), wavelet frames are relevant for building parsimonious models by reducing several types of lower and higher order dependencies (some statistical properties of wavelet transforms for analyzing non-stationary random processes are given in [START_REF] Atto | \Wavelet packets of nonstationary random processes: Contributing factors for stationarity and decorrelation[END_REF]).

The main steps of this modeling are described in the block diagram given by Fig. 2 and are developed below.

B. SWT coecients -FDCT subimages

Among the variants of wavelet frames, we consider the SWT because of its suitable statistical properties [START_REF] Coifman | Translation invariant denoising[END_REF], [START_REF] Atto | \How to perform texture recognition from stochastic modeling in the wavelet domain[END_REF], [START_REF] Atto | \Structuring of large and heterogeneous texture databases[END_REF] and the FDCT for its sensitivity to image geometry: the transform is a quasi-optimal representation for 2 di-TABLE I Average retrieval rates of some distribution families over the SWT and FDCT contributions for two TS-X time series over the Argenti ere glacier. A retrieval rate per family indicates the percentage of occurrence of the distributions pertaining to this family as those with the smallest Kolmogorov stochasticity parameter, among the distribution families considered.

The sequences of TS-X images considered are 13 images in ascending orbit and 10 images in descending orbit. mensional regular signals admitting some singularities like edges and contours for natural images [START_REF] Cand | \Curvelets -a surprisingly eective nonadaptive representation for objects with edges[END_REF], [START_REF] Cand | \New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF], [START_REF] Cand | \Fast discrete curvelet transforms[END_REF].

The J-level SWT and FDCT representations of an image I have the following form:

I = k 1 ,k 2 I, Φ J,[k 1 ,k 2 ] Φ J,[k 1 ,k 2 ] + J j=1 n 1 ,n 2 k 1 ,k 2 I, Ψ j,[n 1 ,n 2 ],[k 1 ,k 2 ] Ψ j,[n 1 ,n 2 ],[k 1 ,k 2 ] . (1)
In this expansion, Φ J,[k 1 ,k 2 ] denotes translated versions of a scaling function Φ J and the wavelet/curvelet waveforms

Ψ j,[n 1 ,n 2 ],[k 1 ,k 2 ]
are generated at decomposition level j from the scaled, translated/oriented versions of a waveform Ψ.

When the SWT is concerned, I will be denoted by I W and the SWT approximation coecients are denoted by

c J,0 [k 1 , k 2 ] = I, Φ J,[k 1 ,k 2 ]
whereas the SWT detail coecients are obtained from

c j,n [k 1 , k 2 ] = I, Ψ j,n,[k 1 ,k 2 ] .
In the above notations, n ∈ {0, 1, 2, 3} is a capital parameter associated with indices n 1 , n 2 ∈ {0, 1}. We will use the standard terminology subband wavelet coecients to designate the subimage of coecients

c j,n = (c j,n [k 1 , k 2 ]) (k 1 ,k 2 )∈Z×Z .
Parametric modeling of this SWT representation concerns subband coecients c j,n . In the experimental results, we consider an SWT with decomposition level is J = 4 and the \Symlet" wavelet of order 8 (a wavelet with high order have more stringent statistical properties with respect to the parametric modeling concerned, see [START_REF] Atto | \Wavelet packets of nonstationary random processes: Contributing factors for stationarity and decorrelation[END_REF], [START_REF] Atto | \How to perform texture recognition from stochastic modeling in the wavelet domain[END_REF]).

When the FDCT is concerned, I will be denoted by I C . In FDCT implementation [START_REF] Cand | \Fast discrete curvelet transforms[END_REF], detail subbands located at a given decomposition level j are generated from scaled and oriented versions of the curvelet waveform Ψ, where orientations are associated with a rotation by an angle θ n = 2π2 -j/2 n, with n ∈ 0, 1, . . . , 2 j/2 -1 (with the correspondence n 1 , n 2 → n, as above).

FDCT coecients form a huge number of subbands, many of these subbands having a small number of coefcients when the decomposition level and the number of orientations are large. Thus, we need to pre-process FDCT in order to make suitable, the curvelet based parametric modeling. In this respect, we consider synthesizing from the curvelet coecients, an approximation and some orientation subimages with the same sizes as the input image.

The curvelet approximation subimage is computed by setting to zero all the curvelet coecients except the approximation ones and by applying the \inverse" FDCT to reconstruct from these coecients:

I app C = k 1 ,k 2 I, Φ J,[k 1 ,k 2 ] Φ J,[k 1 ,k 2 ] .
A curvelet orientation subimage is obtained by aggregating the dierent curvelet subbands lying in a given angular wedge and reconstructing the corresponding subimage by using the inverse FDCT. Assume that the total number of orientations at coarse scale is Q. Then the following steps are used for selecting the FDCT coecients involved in a subimage synthesis. Computation of curvelet orientation subimages: For every q ∈ {0, 1, . . . , Q -1}:

• Identify for every j = 1, 2, . . . , J, the set of polar wedges:

A[q] = {(j, n) : [θ n (j), θ n+1 (j)[ ⊂ [θ q , θ q+1 [} .
• Then, set to zero all the curvelet coecients, except those pertaining to the above wedges, and reconstruct the subimage I ang C [q] by using the \inverse" FDCT:

I ang C [q] = (j,n)∈A[q] k 1 ,k 2 I, Ψ j,n,[k 1 ,k 2 ] Ψ j,n,[k 1 ,k 2 ]
It follows from the above decomposition and by taking into account (1) that

I C = I app C + q∈{0,1,...,Q-1} I ang C [q] = I.
In the following experimental results, we consider a number Q = 8 orientations.

C. Parametric modeling of the SWT and FDCT contributions The parametric modeling is addressed by using 1) Edgeworth expansions for describing SWT and FDCT

[I app C ] 1 p 4 {κp[c J,0 ]} 1 p 4 GG / LOGN / WBL ρ 1 I ang C [q] , ρ 2 I ang C [q] q=0,1,...,Q-1 ρ 1 c j,n , ρ 2 c j,n 1 j J n∈{1,2,3} Number of parameters 2 × Q + 4 6 × J + 4
approximations and 2) a dictionary D composed with the Generalized Gaussian (GG), Log-Normal (LOGN) and Weibull (WBL) distributions for describing the amplitudes of SWT and FDCT details and orientations. The analytic forms of these distribution functions are given in Appendix A.

The SWT approximation subband c J,0 and FDCT approximation subimage I app C have specic behavior [START_REF] Atto | \Wavelet packets of nonstationary random processes: Contributing factors for stationarity and decorrelation[END_REF] because these coecients result from inner products with scaling functions (a scaling function has no vanishing moments). Finding a class of parametric models for the approximation contributions is very intricate because of this specicity. In this respect, we consider the Edgeworth expansion of order p = 4 in order to capture the information conveyed by the rst 4 cumulants of SWT c J,0 and FDCT I app C approximations. For high frequency SWT and FDCT contributions (details and orientations), we focus on GG, LOGN and WBL distributions from model validation results that are given in Table 1 I. This model validation has been performed upon distribution families with exponential decays so as to comply with the statistical properties of wavelet-based transforms. From this model validation step, it follows that: when compared with Gamma (GAM), Generalized Extreme Values (GEV) and Exponential (EXP) distributions, the GG, LOGN and WBL distributions yield more than 99% of retrieval results when the query concerns the \distribution with the smallest Kolmogorov stochasticity parameter" (see average retrieval results given in Table I). These 3 distribution families are thus aggregated to form the dictionary D of parametric models used in the sequel.

We then use the following model selection procedure to retrieve the parameters of the best distribution from this dictionary.

D. Model selection

Consider an SWT detail subband c j,n with some n = 0 or an FDCT orientation subimage

I ang C [q] for some q ∈ {0, 1, . . . , Q -1}.
The model selection step involved in the diagram of Fig. 2 consists in deriving from dictionary D, the best 1 One can note the results given in Table I are specic to TS-X images considered in the framework of the SWT and the FDCT. Thus, statistics given in this table may dier, depending on the data type and the transform used. distribution family and the parameters of this distribution for modeling the subband under consideration: this distribution follows from 1) a rst stage that aims at computing 3 maximum-likelihood estimates associated with GG, LOGN and WBL distributions respectively (see Table I for the motivation on selecting these families) and 2) a second stage for selecting the best distribution as the one admitting the smallest Kolmogorov stochasticity parameter among these maximum-likelihood estimates.

The following items pinpoint the steps involved in this model selection:

• Generate the set of parameters

S ML = {(α ML , β ML ), (µ ML , σ ML ), (a ML , b ML )}
of the \bests" GG, LOGN and WBL distributions from Maximum Likelihood (ML) estimation with respect to the set y of coecients analyzed,

• the Kolmogorov stochasticity parameter [START_REF] Kolmogorov | \Sulla determinazione empirica di una legge di distribuzione[END_REF] with respect to the above ML based distributions is derived from:

λ N (y, F ρ 1 ,ρ 2 ) = sup |F y,N (t) -F ρ 1 ,ρ 2 (t)| , (2)
where F y,N is the empirical cumulative distribution function (cdf ) of the N-sample dataset y and F ρ 1 ,ρ 2 is one among the ML based cdfs computed above.

• The relevant distribution is then the distribution

F ρ * 1 ,ρ * 2 with parameters derived from (ρ * 1 , ρ * 2 ) = arg min (ρ 1 ,ρ 2 )∈S ML λ N (y, F ρ 1 ,ρ 2 ).
The distribution type (GG, LOGN or WBL) is the index of the pair (ρ * 1 , ρ *

2 ) in S ML . The parsimonious description represented by the block diagram of Fig. 2 thus consists in a concise waveletcurvelet modeling of any I m by using P = 2Q + 6J + 8 parameters. Table II summarizes this parametric modeling. By taking into account that Q = 8 and J = 4 in the following, we have P = 48 essential parameters per image description. These parameters encompass multiscale and orientation based image features.

It is worth emphasizing that the P parameters associated with the description of I m from the multi-temporal sequence I can be computed and stored as long as the image I m is available, during the acquisition process. Furthermore, as I m is available, the similarity measurements between I m and the sequence {I , m -1} can be computed as well, for an operational purpose. These similarity measurements proposed are developped in the following section.

E. Similarity measurements

We use KL divergence as similarity measure between the dierent statistical distributions issued from dictionary D. The similarity measurements between image features (distribution parameters) are performed by using the approach proposed in [START_REF] Atto | \How to perform texture recognition from stochastic modeling in the wavelet domain[END_REF]. This approach is extended so as to make possible comparing features issued from two dierent parametric families. Furthermore, the approach proposed below also takes into account the modeling of SWT and FDCT approximations by using Edgeworth expansions.

We consider as similarity measure between two random variables X 1 and X 2 distributions f X 1 and f X 2 (being either Edgeworth expansions or parametric models from D), the symmetric KL similarity measure dened as

K(X 1 , X 2 ) = K(X 1 ||X 2 ) + K(X 2 ||X 1 ), where K(X 1 ||X 2 ) is the divergence K(X i ||X j ) = R f X i (x) log f X i (x) f X j (x) dx, i, j = 1, 2.
The similarity measures are applied to compare adjacent SWT subbands or FDCT subimages, the global similarity measure between I m and I being:

K W (I m , I ) = K(c J,0 [I m ], c J,0 [I ]) + j∈{1,2...,J} n∈{1,2,3} K(c j,n [I m ], c j,n [I ]) (3)
for the SWT, where c j,n [I m ], c j,n [I ] are the SWT coefcients of I m and I respectively and

K C (I m , I ) = K(I app C [m], I app C [ ]) + Q-1 q=0 K(I ang C [m][q], I ang C [ ][q]), (4)
when the FDCT is considered (

I ang C [m][q]
denotes the qth oriented curvelet subimage of the m-th image of the sequence I ).

In the following, K(X 1 ||X 2 ) is computed from parametric models associated with the distributions of X 1 and X 2 . When the detail and orientation coecients are under consideration, these distributions follow from a dictionary D composed with the 3 best distribution families given in Table I. Appendices B and C therefore provide the parametric forms of the KL divergences involved in comparing two arbitrary elements of dictionary D. We recall that approximation coecients are processed separately in the sense that they are modeled by using Edgeworth expansions: the close form of the KL divergence between Edgeworth expansions for p = 4 can be found [START_REF] Lin | \Edgeworth approximation of the Kullback-Leibler distance towards problems in image analysis[END_REF], [START_REF] Inglada | new statistical similarity measure for change detection in multitemporal SAR images and its extension to multiscale change analysis[END_REF]. We propose constructing SWT and FDCT MDDMs from:

K W = (K W (I m , I )) 1 m M,1 M and K C = (K C (I m , I )) 1 m M,1 M .
These matrices relate to the dierent cross similarities between pairwise observations (I m , I ) 1 m M,1 M . Matrices K W and K C are symmetric and have the form of the following matrix:

K =     0 • • 0 • • 0 • • 0    
where a star symbol is used to indicate the elements of the second diagonal.

The motivation beyond the use of K is to provide an intuitive and relevant tool for the analysis of image time series in change-image/epoch detection purpose. It is worth emphasizing that matrix K can be used in an operational context: when a new acquisition is provided, updating MDDM K only requires:

1) computing the P parameters associated with the parsimonious description of the current acquisition, 2) computing the KL divergences between the above parameters and those of the previous acquisitions and 3) adding a new row and a new column composed with the above KL divergences to K. Thus the MDDM matrix format suits for long size image time series.

B. Aggregation of SWT and FDCT MDDMs

The 2M 2 SWT and FDCT cross similarity measurements between elements of sequence I will hereafter be displayed by using a single MDDM. Indeed, since we have considered the symmetric Kullback-Leibler divergence, then the corresponding SWT and FDCT MDDMs are symmetric, so that only the upper triangular MDDM matrix is necessary and sucient for information processing, per transform. Aggregating these MDDMs can thus be performed by superimposing information from SWT and FDCT MDDMs: by replacing the redundant information issued from a single SWT MDDM with complementary Inputs

I[t k ] I[tm] ↓ Wavelet representations I W,k = SWT[I k ] I W,m = SWT[Im] ↓ ' & $ % Extraction of Description parameters                                κ [c J,0 [I k ]] {ρ [c 1,n [I k ]]} n=1,2,3 {ρ [c 2,n [I k ]]} n=1,2,3 . . . {ρ [c J,n [I k ]]} n=1,2,3                                                               κ [c J,0 [Im]] {ρ [c 1,n [Im]]} n=1,2,3 {ρ [c 2,n [Im]]} n=1,2,3 . . . {ρ [c J,n [Im]]} n=1,2,3                                ↓ ' & $ % Multiscaled KL measurements                                K (c J,0 [I k ] , c J,0 [Im]) {K (c 1,n [I k ] , c 1,n [Im])} n=1,2,3 {K (c 2,n [I k ] , c 2,n [Im])} n=1,2,3
. . .

{K (c J,n [I k ] , c J,n [Im])} n=1,2,3                                ↓ operator K W (I k , Im) K k,m ↓ Curvelet representations I C,k = FDCT[I k ] I C,m = FDCT[Im] ↓ ' & $ % Extraction of Description parameters                                κ[I app C,k ] ρ I ang C,k [0] ρ I ang C,k [1]
. . .

ρ I ang C,k [Q -1]                                                               κ I app C,m ρ I ang C,m [0] ρ I ang C,m [1]
. . .

ρ I ang C,m [Q -1]                                ↓ ' & $ % Oriented KL measurements                                K(I app C [t k ], I app C [tm]) K(I ang C [k][0], I ang C [m][0]) K(I ang C [k][0], I ang C [m][1])
. . . information issued from the FDCT MDDM (redundancy follows from the symmetry of the MDDM), we derive the joint SWT -FDCT MDDM K W,C , with

K(I ang C [k][0], I ang C [m][Q -1])                                ↓ operator K C (I k , Im) K m,k
K W,C = K W (I m , I )1l {m } + K C (I m , I )1l {m } 1 m M 1 M
In MDDM K W,C , the upper triangular divergence matrix represents SWT similarity measurements as expressed in (3) and the lower triangular divergence matrix provides FDCT similarity measurements from (4). Fig. 3 provides an overview of the dierent steps involved in the MDDM K W,C computation.

In matrix K W,C , the SWT and FDCT MDDMs are fused (lossless fusion) in order to have a complementary viewpoint that can be helpful for decision taking in many situations. In what follows, representation K W,C will be used for plotting MDDMs whereas analysis of the time series with respect to the MDDM properties will be addressed with respect to one among MDDMs K W and K C .

C. Change-image detection from MDDM analysis

Let us focus on an MDDM K being either K W or K C . The analysis of I derives from the construction of MDDM K:

• In absence of signicant changes from the observed scene and the acquisition system, matrix K does not exhibit any particular structure: it represents only acquisition noise. • In presence of a gradually maintained change (Example: deforestation), -the second diagonal of K, as well as, -each semi-row of K starting from a second diagonal element (that is {K m, : = m + 1, . . . , M}), are increasing sequences, if we assume that the different images are conform, i.e. the acquisitions have been achieved in the same conditions.

• An abrupt change occurs when every semi-row of K starting from a second diagonal element can be roughly associated with a step-like function in the following sense: the divergence values across the semirow form two clusters (or stationary subsequences), a rst cluster with small divergence values representing only acquisition noise (before changes) and a second cluster with large divergence values. The location in time of the changes lies in-between the two clusters.

• In most cases, the analysis of the neighborhood of a divergence values given in the MDDM characterizes and helps understanding the changes or the acquisition conditions at the corresponding acquisition time.

Fig. 4 presents an illustration of progressive changes from an MDDM K W,C computed from a sequence of Landsat images. Images given in Fig. 4 show a part of the Rondonia forest (Brazil) that has been subject to a drastic deforestation (more than 15000km 2 devastated per year). The images clearly highlight through a visual comparison, the progressive changes that have led to the forest destruction. These changes and their progressive nature appear in MDDM 2 K W,C : the second diagonals and any of the semi-column/row starting with a diagonal element are increasing sequences. The abrupt nature of the global change also appear by noting that the last column (respectively row) has very large K W,C values in comparison to the other columns (respectively rows). 2 The MDDM colors have been obtained from the following convention: We rst apply the normalization:

K W ≡ K W / max[K W ] and K C ≡ K C / max[K C ]
where max[M] denotes here the maximum over elements of the matrix M. We then construct the map K W,C from the above normalized similarity measurements. The color map then ranges from blue (associated with 0) to red (associated with 1) by passing through cyan, yellow, and orange. The same convention is used for all the MDDM color compositions given in the paper. A non-conform acquisition (image-outlier) is characterized by a column of K with large divergence values, among the other divergence values representing acquisition noise. Note that the -th column K[ ] of K represents dissimilarities of acquisitions in I with respect to the acquisition performed at time t . As a consequence, an index relating the non-conformity (or inconsistency) of a given acquisition is the cumulative dissimilarities from the columns of K: For illustrating the role of index D, let us consider the Rondonia images of Fig. 4: Table III provides index D measurements from the MDDM of the corresponding images. It appears clearly that the image acquired in 2001 diers signicantly from the other images of the Rondonia time series. Note that for this trivial example, the above remark can be observed directly from the Rondonia optical image time series or its MDDM K.

D[ ] = K[ ] 1 = M m=1 K m, . (5)
The analysis of SAR time series such as the one given in Fig. 5 is more intricate. Images in Fig. 5 are ERS 2 and ENVISAT SAR images showing the Arcachon Bay (France). In addition with the noise induced by acquisition conditions, other disturbances characterizing this scene are swell movement, sand displacement, vegetation variability and the presence of boats depending on the acquisition time. The MDDM K does not highlight any suspicious behavior, except that acquisition 5 (July 15, 2005) is not conform (the column 5 and row 5 have rather large values), in comparison with the other elements of the Arcachon time series. In particular, Table III provides the 3 acquisitions that are less conform from the analysis with respect to index D (colors red, blue and green in Table III).

The following section provides experimental results on the application of K and D for the analysis of image time series involving 13 and 11 acquisition dates respectively.

IV. Analysis of TS-X image time series by using

MDDMs

In this section, we consider two high resolution SAR image time series with respectively 13 and 10 images. The problem tackled is glacier monitoring in a change detection context. The sequences considered are TS-X images of Chamonix Mont-Blanc test site in the French Alps. Among these mountains, we focus hereafter on the Argenti ere glacier.

We rst consider a time series acquired from May 31, 2009 to October 21, 2009, in ascending orbit. The dates of acquisition are given in Fig. 6. The MDDM K W,C is provided in the same gure. As it can be seen from this MDDM, the upper right wedge of K has large divergence values whereas the second diagonal reects only acquisition noise: changes are progressive and follows from climatic seasonal variation that impact the consistency of the primary layer of the glacier. For change evaluation in the spatial domain, acquisitionsof-interest from the above MDDM are the rst and the last acquisitions from the image time series (dates 05-31 and 21-10).

The cumulative dissimilarity indices We now consider a time series of TS-X images of the Argenti ere glacier acquired in descending orbit (acquisition dates are given in Fig. 7). Images are displayed in the SAR slant range geometry.

This sequence extends to one year, with the specicity that no images are available between 2008-10-22 and 2009-05-28 (period including the winter). As it can be seen from the MDDM K W,C of Fig. 7, this \disconnection" have been well detected. Indeed, this \disconnection" involves two dierent seasons where the glacier state is such that the rst layers of the glacier have dierent backscattering properties. From the analysis of this MDDM, the acquisitions we will focus on are those performed at 2008: 09-29 (I 1 ), 05-29 (I 4 ) and 08-25 (I 10 ).

For this TS-X sequence, SWT and FDCT indices D show the same behavior in terms of the stationary sequences (see Table V), but have dierent ranking of cumulative deviations. The \mean" and the \max" fusion rules applied to SWT and FDCT indices D emphasize that FDCT results are more relevant than the SWT ones for analyzing this time series. This can be explained by refering to the suitability of FDCT for representing geometrically regular objects. The above results are instructive in the sense that it reveals a fact that is specic to the acquisition mode and the topography of this Argenti ere's area: glacier surface and geometry is more visible in descending acquisitions due to the glacier slope and valley main orientation. In this respect, we can recommend analyzing the FDCT based modeling for the TS-X descending mode when observing the Argenti ere glacier.

Note that, from Table V, one can distinguish mainly 2 changes in time which separate the acquisitions before 2009-05-29 and after 2009-06-20, the corresponding subsequences having rather stationary cumulative divergences.

V. Spatial change detection from dissimilarity maps between stationary subsequences

The analysis performed in Section IV above makes it possible to detect change-images in image time series. From this analysis, certain image time series can be seen, coarsely, as composed of a stationary subsequences inbetween the corresponding change-images.

The following addresses spatial change detection between the stationary subsequences identied from a time series. We begin by presenting the basics of the method through an example in Section V-A. We then apply this method for analyzing subsequences issued from the TS-X time series provided by Fig. 6 and Fig. 7 in Section V-B.

A. Spatial change detection: basics Specically, we consider images I 2 and I 5 of the sequence of the Argenti ere glacier described in Fig. 7. These images are TS-X images of the Argenti ere glacier (French Alps) acquired at dates 2008-10-10 and 2009-06-09, respectively. Images I 2 and I 5 are with sizes 3072 × 4864.

Change detection between I 2 and I 5 , i.e. in the observation space, involves evaluating local dissimilarities between these images. This problem is addressed hereafter by using the parsimonious descriptions and the similarity measures given in Section II. Indeed, for high resolution images such as TS-X images, the (local) information on a region of interest spreads over a large size pixel neighborhood. This makes relevant, the use of the statistical tools presented in Section II for retrieving and comparing local image features.

In particular, when considering the SWT and FDCT approximation modeling, this relevancy can be highlighted by the following remark.

In [START_REF] Inglada | \Change detection on SAR images by using a parametric estimation of the kullback-leibler divergence[END_REF], the information issued from the 4 rst cumulants3 is shown to be suitable for the description of 8bit RADARSAT images. In contrast, when 16-bit high resolution TS-X images are concerned, Fig. 8 shows that the information derived from the 4 rst cumulants 4 is not sucient to capture the image dynamic. These 4 rst cumulants (relating the mean, standard deviation, asymmetry and kurtosis) are rather suitable for modeling SWT and FDCT approximations, due to the fact that statistical distribution regularizes in the wavelet domain (iterative convolution properties).

M 1-r (I 2 , I 5 ) M K W (I 2 , I 5 )
Fig. 9. TS-X images of the Argenti ere glacier in descending orbit and their dissimilarity maps M 1-r , M K W , where r denotes the use of the correlation coecient and K W relates to the measure given by (3). For every method, the dissimilarity values are linearly scaled to form a colored image where the map of colors ranges from blue (minimum) to red (maximum) by passing through the colors green, yellow, and orange. Changes appear in red color whereas non-changed areas appear in blue. The bi-date correlation is approximately the same all along the Argenti ere glacier so that changes seems uniform on the glacier: only second order changes are captured by the correlations coecients. We need to take into account high order bi-date dependencies from M K W in order to discriminate the eective changes.

Furthermore, this relevancy is justied by noting that the parsimonious description under consideration encompasses both low and high order cumulants of the input random eld. In this respect, the characterization induced by the parameters of this description is more relevant than using a single variable such as the correlation coecient. Indeed, the latter parameter mainly reveals second order bidate dependencies, when applied locally on any pair (I m , I ) of observations, whereas second and higher order bidate dependencies are considered when using the parsimonious descriptions of local neighborhoods of (I m , I ).

For illustration, Fig. 9 provides dissimilarities between I 1 and I 2 by using measurements from M K W and from 1 -r, where M K W (I m , I ) (resp. r(I m , I )) denotes the MDDM obtained from measurements by using (3) and (resp. correlation coecient) between I m and I . These measurements are computed from sliding and overlapping windows with sizes 64 × 64, the overlap induced by the sliding window beginning at the center of the initial window.

As it can be seen in Fig. 9, the correlation between I 1 and I 2 is too weak all over the Argenti ere glacier, due to the speckle decorrelation. It follows that correlation is not relevant for assessing change signicance in this context. In contrast with the correlation coecient, K W highlights change amplitudes depending on the mountain's geometry: accumulation areas, seracs, . . . . This relevancy follows from the capability of the multiscaled and oriented SWT and FDCT similarity measurements in K W for discriminating dierent change levels which aect both correlation and higher order statistics.

B. Spatial change detection: similarities from acquisitions-of-interest identied in TS-X time series

From the MDDM of Fig. 6, we can derive mainly two acquisitions-of-interest; those performed at dates 2009-05-31 and 2009-10-21, respectively. The corresponding images are given in Fig. 10.

These acquisitions-of-interest reect the state of the glacier over two \extremal" behaviors, as it can be seen by comparing image I 1 (dusky) and image I 13 (with high luminance) provided in Fig. 10. In particular,

• the \dusky" state characterizing the glacier in May is conrmed by comparing image I 1 of Fig. 10 (acquired in ascending orbit at May 31) with I 5 given in Fig. 9 (acquired at 9 days interval, in descending orbit). • the \high luminance" characterizing the glacier state in October is assessed by comparing image I 13 of Fig. 10 and image I 2 given in Fig. 9: these images are acquired at the same time period (one year period and 11 days interval). Thus, the state of the glacier during short-time stationary sequences is characterized by images that look similar in terms of intensity variations.

Spatial inter-seasonal dissimilarities can then be evaluated from measurements M(•, •) computed over overlapping windows extracted from the acquisitions-of-interest I 1 , I 13 : these dissimilarities are given in Fig. 10. It follows from M(I 1 , I 13 ) that the main dierence between the glacier layer composition at May and October is concentrated at the upper part of the glacier in the acumulation area, where fresh snowfall modies the surface scattering (see [START_REF] Fallourd | \Monitoring temperate glacier displacement by multi-temporal TerraSAR-X images and continuous GPS measurements[END_REF]).

From the MDDM given in Fig. 7, we can derive 3 acquisitions-of-interest, those performed at dates 2008-09-29, 2009-05-29 and 2009-08-25. These acquisitions and their bidate dissimilarity maps are given in Fig. 11. One can conclude, by comparing the 3 dissimilarity maps M(I 1 , I 4 ), M(I 4 , I 10 ) and M(I 10 , I 1 ) that

• changes that do not appear in the 3 dissimilarity maps are occasional or transient. They follow (mainly) climatic conditions which are intrinsic to the dierent seasons (consequences of the winter and the summer in a high mountainous glacier area) and depend on 1) the glacier elevation and 2) the glacier slope orientation: in the change maps given in Fig. 11, the change distributions in amplitude relate to the form and the geometry of the glacier. • Permanent changes may refer to unstable zones like moraines or to continuous displacements of some glacier blocks, for instance in crevasses or serac fall areas.

VI. Conclusion

In this paper, we have proposed a new approach for analyzing SAR image time series. This approach relies on the computation of a multi-date divergence matrix (MDDM) with respect to some cross divergences issued from the sequences of SAR image description parameters.

Description parameters are issued from SWT and FDCT subband modeling. The modeling has been addressed by using Edgeworth expansions for approximation subbands and GG, LOGN and WBL statistical distribution families for detail and orientation subbands. The approach is multi-model (parametric modeling is not restricted to a single distribution family) and the parametric computation of Kullback-Leibler divergence has been extended to measure inter-family dissimilarities (between distribution families GG, LOGN and WBL).

This approach is computationally workable for long size image time series and makes their analysis straightforward from an MDDM display (that can support series with more than 1000 images). This analysis have shown relevancy in monitoring the temporal evolution of Alps glaciers at two dierent levels: rstly at the image (or sub-image) level, by analyzing the proposed MDDM and cumulative dissimilarity indices that reveal non-conform images and stationary subsequences; then at the pixel level, by computing dissimilarity maps between images of interest selected from the MDDM, for a spatial analysis of changing areas.

The paper opens some prospects relating the statistics of MDDMs for image time series analysis. Indeed, an MDDM makes possible both global and local analysis of a time series: encapsulating global and local MDDMs in the neighborhood of a critical scene will ease a multilevel analysis of the scene under consideration. In addition, beyond the multi-date visual scene observation provided by MDDMs, an MDDM representation makes also possible variance analysis in the time axis, searching for the median of the image time series, etc. Temporal segmentation tools could also be developed to analyze MDDMs and derive sub-matrices associated to seasonal or annual evolutions. This can be derived by applying matrix norms to the series of MDDMs obtained. This issue is challenging with respect to future SAR missions like Sentinel 1 to be launched in 2013, with more systematic acquisition strategies (every 12 days, then every 6 days with two satellites) and with free data access.

The GG, LOGN and the WBL distributions are used in a joint framework: this includes the case where a GG distribution has to be compared with a LOGN or a WBL distribution, and vice versa. In this respect (a LOGN and a WBL distribution operating on positive real values), we consider for comparison purpose, the magnitude of a random variable following the GG distribution: the magnitude of a GG distribution with scale α, shape β is given by:

f α,β (x) = β αΓ (1/β) e -(x/α) β 1l R + (x), (7)
where Γ is the standard Gamma function:

Γ (u) = R + x u-1 e -x dx
Note that this expression diers from the standard GG density by a factor 2.

The LOGN distribution with log-scale parameter µ > 0 and shape σ 2 is dened by:

f µ,σ (x) = 1 xσ √ 2π e -( (log x-µ) 2 2σ
2

) 1l R + (x). (8)
The WBL distribution with scale a > 0 and shape b > 0 is dened by:

f a,b (x) = b a x a b-1 e -( x a ) b 1l R + (x). (9) Appendix B
Intra-family KL divergences for GG, LOGN and

WBL distributions

The expression of the KL measure for distributions pertaining to the the same parametric family is available for many distribution families from the literature.

For GG distributions, the KL divergence is given in [START_REF] Do | \Wavelet-based texture retrieval using generalized gaussian density and kullback-leibler distance[END_REF]. The symmetric version of this divergence has the following form:

K(X 1 , X 2 ) = α 1 α 2 β 2 Γ 1+β 2 β 1 Γ (1/β 1 ) + α 2 α 1 β 1 Γ 1+β 1 β 2 Γ (1/β 2 ) - β 1 + β 2 β 1 β 2 . (10)
For Weibull distributions, the KL divergence is given in [START_REF] Kwitt | similarity measurement by Kullback-Leibler divergences between complex wavelet subband statistics for texture retrieval[END_REF]. We have that the symmetric version of this divergence is

K(X 1 , X 2 ) = λ 1 λ 2 k 2 Γ 1 + k 2 k 1 + λ 2 λ 1 k 1 Γ 1 + k 1 k 2 +e k 1 k 2 + k 2 k 1 -2 -2 +(k 1 -k 2 ) log λ 1 λ 2 . ( 11
)
where e is the Euler-Mascheroni constant.

Concerning LOGN distributions, the explicite form of the parameter-dependent KL divergence derives from:

Proposition 1: The symmetric KL divergence between LOGN random variables X 1 and X 2 is given by: K

(X 1 , X 2 ) = 1 2 (µ 1 -µ 2 ) 2 1 σ 2 1 + 1 σ 2 2 + 1 2 σ 2 1 σ 2 2 + σ 2 2 σ 2 1 -1. (12)
Proof: Proposition 1 derives from straightforward calculus, by expanding the KL integral form as a sum of Gaussian moments.

The following Appendix C derives close forms of the KL divergence for the joint framework involved in the use of a dictionary D. This joint framework requires comparing distributions issued from dierent parametric families (Inter-family similarity measurements).

Appendix C

Inter-family KL divergences between GG, LOGN and WBL distributions Consider two random variables X 1 , X 2 with probability density functions f X 1 , f X 2 respectively. The symmetric KL similarity measure between X 1 and X 2 is:

K(X 1 , X 2 ) = H(X 1 ||X 2 )-H(X 1 )+H(X 2 ||X 1 )-H(X 2 ) (13)
where the entropy H is dened for a random variable X with probability density function f X by H(X) = -R f X (x) log f X (x)dx [START_REF] Cand | \Fast discrete curvelet transforms[END_REF] and H(X i ||X j ) denotes the entropy of random variable X i conditional on the random variable X j . This conditional entropy is

H(X i ||X j ) = - R f X i (x) log f X j (x)dx (15)
when X i , X j admit probability density functions f X i , f X j respectively.

A. KL Divergence between a Generalized Gaussian Magnitude and a Weibull distribution Let X w , X g be random variables distributed as WBL and GG with probability density functions f a,b , f α,β respectively. Proposition 4 below provides the symmetric KL similarity measure between X w and X g . This proposition requires the following Propositions 2 and 3 relating the cross entropies between X w and X g . We have

where ψ is the digamma function satisfying d dx

Γ (x) = ψ(x)Γ (x). (18)
Proof: See Appendix E.

From Propositions 2, 3, and by taking into account that the entropy of the magnitude of a GG random variable is

H(X g ) = 1 β + log 1 β + log αΓ 1 β (19)
and the entropy of a WBL random variable is

H(X w ) = e 1 - 1 b + log a b + 1, (20)
we have:

Proposition 4: The symmetric KL divergence between a WBL and the magnitude of a GG random variables is given by: K(X w , X g ) = . [START_REF] Mccullagh | Tensor Methods in Statistics[END_REF] Proof: See Appendix G.

From Propositions 5, 6, Equation [START_REF] Inglada | new statistical similarity measure for change detection in multitemporal SAR images and its extension to multiscale change analysis[END_REF] and by taking into account that the entropy of a LOGN random variable X is

H(X ) = 1 2 + µ + 1 2
log(2πσ 2 ), [START_REF] Stuart | Kendall's Advanced Theory of Statistics[END_REF] we obtain:

Proposition 7: The symmetric KL divergence between an WBL and a LOGN random variables is given by: K(X w , X )= e bµ+b 2 σ 2 /2 a 2 + 1 2σ 2 × e 2 + π

C. KL Divergence between a Generalized Gaussian Magnitude and a Log-normal distribution Consider now a GG magnitude and a LOGN random variables X g , X with probability density functions f α,β , f µ,σ respectively. The following provides cross entropies and KL divergences between these random variables. From Propositions 8, 9, Equations ( 19) and ( 24), we have:

Stage 1 /Fig. 1 .

 11 Fig.1. Block diagram representing temporal and spatial image similarity measurements. In this diagram, transform T designates either the SWT or the FDCT. Specically, the term KL(T X/T Y) designates KL based similarity measurements from the SWT or the FDCT decompositions of X and Y, see for instance (3) and (4). Temporal analysis involves large size image (or subimage) comparison whereas spatial analysis applies only for specic images at the pixel level. The MDDM shown in the center of the gure has been computed upon a multi-date Gaussian random eld composed with 3 stationary subsequences.

  II. Feature extraction from SWT and FDCT parametric modeling -Similarity measurements A. Motivation Let I = {I m , m = 1, 2, . . . , M} be a sequence of M images of a given scene acquired at dierent dates:

Fig. 3 .

 3 Fig.3. Flowchart representing the steps involved in MDDM value computation. In this owchart, vector κ = (κ ) =1,2,3,4 designate the rst four cumulants and vector ρ = (ρ 1 , ρ 2 ) represents the scale and shape parameters of the best model obtained from D.

LandsatFig. 4 .

 4 Fig. 4. Multi-date image sequence from Rondonia forest, NASA/courtesy of nasaimages.org, and its MDDM K W,C . The upper (resp. lower) triangular matrix of the MDDM represents the SWT (resp. FDCT) contributions in K W,C . The rst row and the second diagonal of K W,C are increasing sequences: the changes are progressive and maintained. In addition, the changes from the last column of K W,C appears abrupt as corroborated by the visual image comparison. Images are with size 1920 × 2880.

  D. MDDM and the non-conformity of an acquisition with respect to the image time series ERS 2, June 06, 2003 ENVISAT, August 15, 2003 ENVISAT, June 25, 2004 ERS 2, June 10, 2005 ENVISAT, July 15, 2005 ENVISAT, June 30, 2006 K, Arcachon Bay

Fig. 5 .

 5 Fig. 5. Multitemporal ERS 2 and ENVISAT image sequence, CNES/courtesy. The MDDM K W,C only reects acquisition noise and an outlier (ENVISAT acquisition of 2005, July 15. Images are with size 3840 × 4480.

KFig. 6 .

 6 Fig.6. MMDM K W,C for a multi-date TS-X image sequence in ascending orbit. The changes from the MDDM are progressive from July to September and reect the glacier surface evolution.

KFig. 7 .

 7 Fig. 7. MDDM K W,C for a sequence of TS-X images acquired over Argenti ere glacier in descending orbit. Changes are abrupt (the rst row of the KL MDDM has a step located at the acquisition date 4 whereas the second diagonal has a unique outlier located at date 4) and progressive (decay of semi-row sequences after acquisition date 4). Images are with size 3072 × 4864.

  Argenti ere \descending", 2008-2009 2008 2009 Dates 09-29 10-10 10-21 05-29 06-09 06-20 07-01 07-12 08-14 08-25 D W [ ]

Fig. 8 .

 8 Fig. 8. Experimental setup: Edgeworth (parametric) expansion of order 4 is applied for modeling a 64 × 64 TS-X bloc (over the Argenti ere glacier serac area). The modeling performed in the image domain (top) is irrelevant whereas the modeling is relevant for the SWT (middle) and FDCT (down) approximations of the input TS-X bloc. The comparison of the Edgeworth expansion is made with respect to a smoothed histogram obtained by using (non-parametric) kernel density estimation (in red).

I 2 ,

 2 2008-10-10 I 5 , 2009-06-09

I 1 , 21 M(I 1 , I 13 )Fig. 10 .

 12111310 Fig.10. Dissimilarity map between two images issued from the TS-X sequence given in Fig.6. The acquisitions-of-interest are the rst and last acquisitions from the sequence. Images are with size 4096 × 3072.

I 1 ,Fig. 11 .

 111 Fig. 11. Change detection for some acquisitions-of-interest issued from the Argenti ere glacier TS-X sequence of Fig. 7. Images are with size 3072 × 4864.

Proposition 2 :

 2 The conditional entropy H(X w ||X g ) is given byH(X w ||X g ) =

Proposition 3 :

 3 The conditional entropy H(X g ||X w ) is given byH(X g ||X w ) =

2 log(2πσ 2 ) + log a + e b + e 2 + π 2 Proposition 6 :

 22226 Divergence between a Log-normal and a Weibull distributionTo conclude this section, the following provides cross entropies and KL divergences between a LOGN and a WBL random variables X , X w with probability density functions f µ,σ , f a,b respectively. Proposition 5: The conditional entropy H(X w ||X ) is given by H(X w ||X ) = 1 The conditional entropy H(X ||X w ) is given by H(X ||X w ) = b log a-log b+µ(1-b)+ e bµ+b 2 σ 2 /2 a 2

2 6 b 2 +

 2 µ 2 + (µ -log a) µ -log b -

Proposition 8 : 2 log 2 2σ 2 β 2 ( 26 )Proposition 9 :+ 1 α β e µβ+β 2 σ 2 / 2

 8222269122 The conditional entropy H(X g ||X ) is given byH(X g ||X ) = 1where ψ(2) denotes the standard Trigamma functionsatysfying d dz 2 Γ (z) = Γ (z) ψ (2) (z) + (ψ(z)) 2 . The conditional entropy H(X ||X g ) is given by H(X ||X g ) = log α-log β+Λ 1 β (28)where function Λ denotes the natural logarithm of the Gamma function.Proof: See Appendix I.

  2. Block diagram representing image feature extraction. Parameters (κ [X]), = 1, 2, 3, 4 designate the rst four cumulants of X: κ 1 and κ 2 are the mean and variance of X, κ 3 and κ 4 relate to the skewness and kurtuosis of X, respectively. These cumulants are associated with an Edgeworth expansion of X. Notation ρ 1 [X] and ρ 2 [X] designate the scale and shape parameters of the best model among the GG, LOGN and WBL distributions associated with |X| and issued from a model selection step.

TABLE II

 II Overview of the parametric modeling and the associated number of parameters.

	Parametric modeling	FDCT	SWT
	Edgeworth	κp

  III. KL based MDDM from SWT and FDCT modeling A. MDDM Consider the M size sequence I given in Section II-A. Section II-E provides 2M 2 SWT and FDCT similarity measurements (with a redundancy factor of 2 due to symmetry).

TABLE III Non

 III -conformity of acquisitions in Rondonia and Arcachon time series. Strong deviation from the normal behavior of the sequence occur for acquisitions with large D values (red, blue and green colors in the table).

			1		2	3	4
	Rondiana D W [ ] 11.54 7.27 13.03 25.32
		D C [ ]	7.41 4.94 8.49 16.31
			Arcachon
		1	2	3	4	5	6
	D W [ ]	5.56 2.89 3.23 3.85 8.82 2.36
	D C [ ]	9.97 4.58 4.81 6.68 15.67 3.94

  It appears from this table that neither image I 13 , nor image I 1 are representative of the glacier state on the acquisition period: D[START_REF] Bruzzone | advanced system for the automatic classication of multitemporal SAR images[END_REF] and D[START_REF] Cand | \New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] deviate signicantly from the rest of the cumulative dissimilarities. Table IV reveals a rather stationary state of the glacier from the period of July 03, 2009 to September 18, 2009. After this stationary period, we observe a signicant increase of the cumulative dissimilarity indices that can be explained by recent snowfalls and/or short-time continuous snow deposit on the glacier at the beginning of the autumn.

D[ ], = 1, 2, . . . , 13 of this time series are given in Table

IV

.

TABLE IV Non

 IV -conformity of acquisitions in the Argenti ere \ascending" time series. Strong deviation from the normal behavior of the sequence occur for acquisitions with large D values (red, blue and green colors in the table). Argenti ere \ascending", 2009 Dates 05-31 06-11 06-22 07-03 07-14 07-25 08-05 08-16 08-27 09-18 09-29 10-10 10-21 Non-conformity of acquisitions in the Argenti ere \descending" time series. Strong deviation from the normal behavior of the sequence occur for acquisitions with large D values (red, blue and green colors in the table).

	D W [ ]	6.76	4.03	4.52	2.07	2.02	2.60	2.10	2.23	2.15	2.43	4.28	2.03 14.14
	D C [ ]	6.94	3.50	4.73	2.09	2.06	2.75	2.20	2.21	2.17	2.52	3.84	2.04 15.93
							TABLE V

2008-09-29 2008-10-10 2008-10-21 2009-05-29 2009-06-09 2009-06-20 2009-07-01 2009-07-12 2009-08-14 2009-08-25

The information from the

rst cumulants is expressed in the Pearson system in[START_REF] Inglada | \Change detection on SAR images by using a parametric estimation of the kullback-leibler divergence[END_REF].[START_REF] Lhermitte | \Hierarchical image segmentation based on similarity of NDVI time series[END_REF] The information from the 4 rst cumulants is expressed in the Edgeworth expansion in this paper.

ACKNOWLEDGMENTS

The authors wish to thank the anonymous reviewers for their precise comments and useful advice on this paper and the German Space Agency (DLR) for providing the TerraSAR-X data (Project MTH0232).

Appendix A Analytic forms of the parametric distribution models

The Edgeworth expansion of order p of a random variable X, absolutely continuous, with mean µ and standard deviation σ is

where H r is the Chebyshev-Hermite polynomial of order r [START_REF] Mccullagh | Tensor Methods in Statistics[END_REF] and the coecient η r is a function of the r rsts cumulants of X. Note that when p = ∞, the right hand side expression of the above equation is exactly the probability density function of X, under regularity assumptions on this function [START_REF] Stuart | Kendall's Advanced Theory of Statistics[END_REF].

Proposition 10: The symmetric KL divergence between LOGN and GG magnitude random variables is given by:

(

We have:

where f α,β (x), f a,b (x) are given by ( 7), ( 9) respectively. From ( 7), ( 9) and ( 30), it follows that

By using integration by parts on the left hand side term of (31), we have

By applying a change of variable on the latter equality, we get

Proposition 2 then follows from (31), (33) and by taking into account the denition of the Γ function (see [START_REF] Park | \Similarity-based subsequence search in image sequence databases[END_REF]).

Appendix E Proof of Proposition 3

We have:

where f α,β (x), f a,b (x) are given by ( 7), ( 9) respectively. By taking into account (7), ( 9) and (34), the conditional entropy H(X g ||X w ) can be expressed in the following form

with

and

We have:

By noting taking into account [START_REF] Kolmogorov | \Sulla determinazione empirica di una legge di distribuzione[END_REF], we obtain

In addition, from a change of variable x -→ (x/α) β , we derive:

Proposition 3 is then a consequence of ( 35), ( 36), ( 37), ( 39) and (40).

Appendix F Proof of Proposition 5

We have:

where f µ,σ , f a,b are given by ( 8), [START_REF] Morio | \Information theory-based approach for contrast analysis in polarimetric and/or interferometric sar images[END_REF], respectively. The sketch of the proof is the following: rst, decompose log f µ,σ in (41) as a logarithmic polynomial. Then, identify the combination of derivatives of the gamma function by expanding the integrand of (41). Expressions ( 8), [START_REF] Kolmogorov | \Sulla determinazione empirica di una legge di distribuzione[END_REF] involving the Gamma and Digamma functions, as well as the Expression (27) involving the Trigamma function are useful in deriving [START_REF] Fallourd | \Monitoring temperate glacier displacement by multi-temporal TerraSAR-X images and continuous GPS measurements[END_REF] of Proposition 5.

Appendix G Proof of Proposition 6

We have:

where f µ,σ , f a,b are given by ( 8), [START_REF] Morio | \Information theory-based approach for contrast analysis in polarimetric and/or interferometric sar images[END_REF], respectively. The integrand in (42) can be expanded as a combination of Gaussian moments. This is performed via a change of variable t = log x. Proposition 6 then follows after some straightforward calculus.

Appendix H Proof of Proposition 8

We have:

where f α,β , f µ,σ are given by ( 7), [START_REF] Park | \Similarity-based subsequence search in image sequence databases[END_REF], respectively. The steps of the proof are similar to those given in Appendix F (proof of Proposition 5).

Appendix I Proof of Proposition 9

We have:

where f α,β , f µ,σ are given by ( 7), [START_REF] Park | \Similarity-based subsequence search in image sequence databases[END_REF], respectively. The steps of the proof are those given in Appendix 6 relating the proof of Proposition 6.


Khoa Truong

Tuan Thanh Phan

Duc Doan Tran

Nam Nguyen

Thoai

Khoa T Phan

Tuan T Tran

Duc D Nguyen

Nam Thoai

MaxNet and TCP Reno/RED on Mice Traffic

Keywords: Active queue management, TCP congestion control, MaxNet

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

TCP Reno [START_REF] Van Jacobson | Berkeley TCP Evolution from 4.3-Tahoe to 4.3-Reno[END_REF] uses AIMD mechanism [START_REF] Van Jacobson | Congestion Avoidance and Control[END_REF] in which the sending rate is increased until packet loss happens. To avoid buffer overflows at router, AQM RED (Active Queue Management Random Early Detection) [START_REF] Floyd | Random Early Detection Gateways for Congestion Avoidance[END_REF] can be used in conjunction with TCP Reno. The weakness of RED is that it does not take into account the number of incoming flows arrived at a bottleneck link to perform appropriate treatments to lighten down the heavy load. When there are a large number of sharing flows at a bottleneck link, the offered load will not be decreased, K.T. Phan T.T. Tran D.D.Nguyen N. Thoai Faculty of Computer Science and Engineering, Ho Chi Minh City University of Technology, Vietnam e-mail: khoaphan@cse.hcmut.edu.vn; nam@cse.hcmut.edu.vn also the queue length at router does not change. This results in longer response time for the mice traffic.

MaxNet [SAW08, [START_REF] Suchara | TCP MaxNet -Implementation and Experiments on the WAN in Lab[END_REF][START_REF] Wydrowski | MaxNet: A Congestion Control Architecture for Scalable Networks[END_REF] is a new congestion control mechanism using multi-bit signal instead of packet loss to control the sending rate. Besides, MaxNet router can control the magnitude of transient queues well regardless the number of new arrival flows. In other words, MaxNet can always keep a free space at routers for mice traffic to fly though. As a result, the response time to HTTP requests is much shorter with MaxNet than with TCP Reno/RED.

The rest of this paper is structured as follows. Section 2 is the theoretical analysis of queueing delay of RED and MaxNet routers. Section 3 shows the efficiency of MaxNet's quick control of the transfer rate of mice flows. We have some experiments and evaluations in Sect. 4. Finally, we present the conclusions and future work in Sect. 5.

Equilibrium Queueing Delay at RED and MaxNet Router

Queueing Delay at RED Router

RED routers calculate and compare the average queue length based on the two parameters: maximum threshold and minimum threshold. Base on this comparison, RED router operates in three modes [START_REF] Floyd | Random Early Detection Gateways for Congestion Avoidance[END_REF]:

-"No dropped": when the average queue length is less than the value of the minimum threshold, router assumes that its link is under-utilized. As a result, all packets are allowed to go through without marking or dropping. -"Probabilistic dropped": when the average queue length is between the minimum and maximum thresholds, router assumes that the network can be saturated, it then marks/drops packets with a probability corresponding to the traffic load. -"Forced dropped": when the average queue length is greater than the maximum threshold, all packets that go through the router will be marked/dropped to reduce heavy load on the link.

From the flow-level model of AIMD, the window size of TCP Reno is updated using the following equation [START_REF] Wei | FAST TCP: motivation, architecture, algorithms, performance[END_REF]: i .t/ D 0, then from (1), the end-to-end equilibrium mark/drop probability feedback to source i is derived as following: Equation (2) implies that, at the equilibrium point, end-to-end mark/drop probability of source must be greater than zero. As a result, from the marking scheme of RED, it can be asserted that each router on the end-to-end always keeps a backlog. This consequently causes inevitable queueing delay for mice traffic such as HTTP requests to fly through. Figure 1 illustrates two RED routers which always maintain backlog at equilibrium point.

w 0 i .t/ D 1 T i .t
q i D 3 2:w 2 i >0 (2)

Queueing Delay at MaxNet Router

The marking mechanism of MaxNet router uses an explicit multi-bit signal instead of marking/dropping packet as RED router's. Congestion price p l at MaxNet router is defined in [START_REF] Suchara | Implementation of Provably Stable MaxNet[END_REF]:

p l .t C dt/ D p l .t/ C dt y l .t/ l C l C l (3)
where y l .t/ is the aggregated rate at link l; C l is the link capacity and l is the target link utilization. In MaxNet router, at the equilibrium point, the link price adjustment (3) tries to match the aggregated input rate y i .t/ with l C l , leaving spare .1 l /C l capacity to absorb mice traffic and reduce the queueing delay. Figure 2 illustrates the queueing delay of two MaxNet bottleneck links at the equilibrium point. In contrast to RED routers, there is no backlog in both two MaxNet routers when the target link utilization is set to l where 0< l <1, hence, mice traffic can fly through links without being blocked.

Magnitude of Transient Queue of RED and MaxNet Routers

As pointed out in [START_REF] Feng | A Self-Configuring RED Gateway[END_REF][START_REF] Le | The Effects of Active Queue Management on Web Performance[END_REF], the weakness of RED is that it does not take into account the number of flows arriving a bottleneck link to have proper treatments to avoid heavy load. Assuming there are n flows sharing a bottleneck link. If a packet is

Mice traffic

MaxNet Router MaxNet Router

Elephant traffic

Fig. 2 Queueing delay of MaxNet routers marked or dropped then the offered load is reduced by a factor of .1 0:5n 1 /.When n is large then .1 0:5n 1 / ! 1, which means the offered load will not be decreased and the queue length doesn't change, either. That means if RED is not configured aggressively then marking a single packet could result in simple "droptail" packet drops. The packet loss then severely declines the throughput of Reno sources due to their AIMD mechanism. MaxNet router well controls the magnitude transient queue regardless the number of new arrival flows. At the equilibrium point, when the number of new arrival flows is small, transient queues exist, but the magnitude of these queues decreases rapidly as the number of flows increases [START_REF] Suchara | Implementation of Provably Stable MaxNet[END_REF]. This can be explained with the following simple case: assuming that there are N flows sharing a bottleneck link. At the equilibrium point, each flow transmits at the rate of i C i N . Thus, when a new flow joins, its advertised rate is at most i C i N . The aggregated arrivals at router are at most l C l C l C l N . Thus, this causes the overload:

0 overload ..1 C 1 N / l 1/C l (4)
Obviously, the larger the N is, the smaller the magnitude of transient queue becomes and eventually when N> 1 , the transient queue size drops to zero. As mice traffic is short-lived flows, an effective congestion control should quickly controls the rate of such flows to avoid uncontrolled overshoot and transient queue. Unlike SlowStart mechanism of Reno [START_REF] Van Jacobson | Congestion Avoidance and Control[END_REF], MaxNet source employs MaxStart mechanism [START_REF] Suchara | Implementation of Provably Stable MaxNet[END_REF] to seek for the target rate at initial state Fig. 3.

By adopting the multi-bit explicit signaling mechanism, MaxStart enables source to seek for its target rate within a significant short duration. New MaxNet flow is initiated at the minimum spare capacity of all links on its end-to-end path and then ramped up linearly to the advertised rate over two Round Trip Time [START_REF] Suchara | Implementation of Provably Stable MaxNet[END_REF]. Therefore MaxNet source converges to the target rate more quickly than Reno source.

Testbed Layout

In this testbed, Pentium IV PCs (CPU 1.8 GHz, 256 MB RAM, 100 Mbps network cards) are used. The MaxNet Router1 is configured with the output capacity 10 Mbps to make sure it is the bottleneck link. The target utilization of both MaxNet routers . l / is set to 0.95. Dummynet [DUM] is configured with 20 ms RTT delay.

The testbed of Reno/RED is same as MaxNet testbed where MaxNet routers are changed to RED routers. All of RED routers are configured with the RED parameters for web traffic [START_REF] Feng | A Self-Configuring RED Gateway[END_REF] as following:

w q D 0:002 : weighting factor for computing average queue size as suggested in [START_REF] Claypool | Active Queue Management for Web traffic[END_REF]forwebtraffic. -q avg D .1 w q /:q avg C w q :q with q is instantaneous queue size.

min th D 30: average queue length threshold for triggering probabilistic drops/marks. max p D 0:1: maximum mark/drop probability.

We simulated elephant traffic and mice traffic with iperf [IPE] and httperf [HTT] tools. For the both networks of Reno/RED and MaxNet, one long live FTP connection is generated at the "FTP source" in approximately of 60 s. Then 20 s after, HTTP connections are generated on the "HTTP source". Each HTTP connection sends one request of 62 bytes and HTTP response size is 4 KB. The HTTP response time is computed at the application layer by the duration from the first byte being sent out to the time when the first byte of response received.

Response Time of HTTP Connections

We adopted the cumulative probability for statistical analysis of the response time of HTTP requests. In Fig. 5, the response time of HTTP requests in MaxNet is significantly less than in TCP Reno/RED. Particularly, in the experiment with 200 HTTP requests spawn, 100% of MaxNet HTTP requests receive the reply at most at 135ms while in TCP Reno/RED, only 30% of the total HTTP requests receive the first reply less than 135ms.

Throughput of Elephant Flow

In Figs. 6 and7, packets drop occur at the RED bottleneck link, thus the throughput of Reno/RED is decreased. The greater the number of connections is, the more severity the drop becomes. In contrast, the throughput of elephant flow in MaxNet networks is not impacted regardless the number of arrival flows.

Conclusions and Future Work

At the equilibrium point, MaxNet can clear the buffer while Reno/RED always keeps a backlog in routers. Therefore, when elephant traffic is mixed with mice traffic, MaxNet has a shorter response time for mice traffic than TCP Reno/RED. If the number of arrival mice flows is large, Reno without proper treatment can cause packet loss which in turn eventually degrades the throughput of elephant traffic. In addition, MaxStart mechanism of MaxNet (using multi-bit signaling) can control mice flows to the target rate more quickly than Reno sources. By the experiments, we showed that the performance of mice and elephant traffic when using with MaxNet is better than with TCP Reno/RED's in terms of response time and network utilization. For future work, more experiments should be conducted with realistic web workload and other performance properties such as fairness, TCP friendly should be evaluated.

Fig. 1

 1 Fig. 1 Queueing delay of RED routers

  Fig. 3 SlowStart and MaxStart

Fig. 5

 5 Fig. 5 Response time of TCP Reno vs MaxNet with 200 HTTP connections

Fig. 6 Fig. 7

 67 Fig. 6 Throughput of elephant traffic when 50 HTTP flows join

Fig. 8 Fig. 9

 89 Fig. 8 Backlog at RED router

Acknowledgements

The authors would like to thank Bartek Wydrowski, Lachlan Andrew and MaxNet team [TEAM] for their advices and support.