Gabriel Vasile
email: gabriel.vasile@gipsa-lab.grenoble-inp.fr

Nikola Besic

Andrei Anghel

Circularity and Sphericity of Complex Stochastic Models in Multivariate High-Resolution SAR Images

This paper presents a new methodological framework to assess the conformity of multivariate high-resolution SAR data in terms of asymptotic statistics. Three important statistical properties are studied by applying statistical hypotheses testing, successively: circularity, sphericity and spherical symmetry. Starting from the classical tests designed for the multivariate Gaussian case, these tests are extended to the Spherically Invariant Random Vector (SIRV) stochastic model. A zero-mean test is proposed for both Gaussian and SIRV stochastic processes. The link between the spherical symmetry property and the conformity to the SIRV model is established asymptotically by the specific structure of the quadricovariance matrix. Two high and very high resolution datasets are used to illustrate departures from the standard model assumptions: TerraSAR-X multipass InSAR and ONERA RAMSES POLSAR images. As well, the derived tests are applied on the appropriate synthetic dataset.

The detection results are qualitatively and quantitatively analysed and some important inferences are drawn regarding these two datasets.

I. INTRODUCTION

Multidimensional Synthetic Aperture Radar (SAR) data, like POLSAR and multi-pass InSAR images, describe the interaction between the electromagnetic waves and the scatterers inside each resolution cell. Multivariate SAR images, acquired by either spaceborne or airborne sensors, are currently analyzed to extract useful information concerning the physical properties of the illuminated target. For distributed targets, this analysis relies on the stochastic properties of the SAR data. In general, multivariate SAR data can be locally modelled by a multivariate zero-mean circular Gaussian stochastic process, which is completely determined by its covariance matrix.

With the improved resolution of modern SAR platforms, the number of scatterers inside each resolution cell decreases considerably. The higher scene heterogeneity may eventually lead to non-Gaussian clutter modelling. More complex stochastic models, such as the SIRP (Spherically Invariant Random Process), are then required. Several special cases of univariate stochastic processes (K-compound, Weibull, etc.) have been extensively studied over the years (for example in coastal radar applications) before being reunited under the common umbrella of SIRP [START_REF] Picinbono | Spherically invariant and compound Gaussian stochastic processes[END_REF].

Polarimetry and multi-pass interferometry extend the dimensionality of the remotely sensed SAR data. Multivariate versions of SIRP distributions, namely SIRV (Spherically Invariant Random Vectors) [START_REF] Yao | A representation theorem and its applications to spherically-invariant random processes[END_REF], are employed for modelling high-resolution POLSAR [START_REF] Vasile | Coherency matrix estimation of heterogeneous clutter in high resolution polarimetric SAR images[END_REF][START_REF] Bombrun | Hierarchical segmentation of polarimetric SAR images using heterogeneous clutter models[END_REF][START_REF] Formont | Statistical classification for heterogeneous polarimetric SAR images[END_REF][START_REF] Vasile | Optimal parameter estimation in heterogeneous clutter for high resolution polarimetric SAR data[END_REF][START_REF] Besic | Stochastic approach in wet snow detection using multitemporal SAR data[END_REF][START_REF] Anghel | Short-range wideband FMCW radar for millimetric displacement measurements[END_REF] and InSAR [START_REF] Vasile | High resolution SAR interferometry: influence of local topography in the context of glacier monitoring[END_REF][START_REF] Julea | Simulation of SAR images and radar coding of georeferenced information for temperate glacier monitoring[END_REF][START_REF] Vasile | Stable scatterers detection and tracking in heterogeneous clutter by repeat-pass SAR interferometry[END_REF][START_REF] Vasile | Multidimensional very high resolution SAR interferometry for monitoring energetic structures[END_REF][START_REF] Vasile | Potential of multi-pass high-resolution SAR interferometry for dam monitoring[END_REF] data. This multiplicative model is expressed as a product between the square root of a scalar positive quantity (texture) and the description of an equivalent homogeneous surface (speckle) [START_REF] Ulaby | Textural information in SAR images[END_REF]. For an m-dimensional SAR system, the single channel model [START_REF] Fjortoft | Estimation of the mean radar reflectivity from a finite number of correlated samples[END_REF] has been extended as follows: in each azimuth / range location, k is the m × 1 complex target vector corresponding to the same area on the ground.

For distributed targets, the corresponding k vector is considered non-deterministic and may be written, under the SIRV assumption as k = √ τ • z, where τ is a positive random variable (texture) while z is a complex-valued, centred and Gaussian distributed, random vector with covariance matrix [M ] = E{zz H }. The two random variables are statistically independent.

Characterizing k reduces to the writing of the probability distribution p τ (τ ) of τ and the normalized covariance matrix [M ] of z, with trace{[M ]} = m. It is important to notice that in the SIRV definition, the texture probability density function is not explicitly specified. As a consequence, SIRVs describe a whole class of stochastic processes [START_REF] Rangaswamy | Computer generation of correlated non-Gaussian radar clutter[END_REF]. If p τ (τ ) is a Dirac pulse, then the Gaussian (multivariate) model is retrieved for k. In many field applications, p τ (τ ) is assumed to have a predefined analytical form. In this case, a specific SIRV is employed such as the multivariate K distribution (Gamma texture) [START_REF] Novak | Optimal speckle reduction in polarimetric SAR imagery[END_REF], the G 0 distribution (inverse Gamma texture) [START_REF] Nascimento | Hypothesis testing in speckled data with stochastic distances[END_REF] or the KummerU distribution (Fisher texture) [START_REF] Bombrun | Hierarchical segmentation of polarimetric SAR images using heterogeneous clutter models[END_REF]. Nevertheless, at fixed number of samples, as the number of parameters increases, the variance of their estimates will increase.

The four components of the POLSAR k target vector (in lexicographic order) describing the polarimetric characteristics of a given resolution cell are the elements of the Sinclair matrix: S hh , S hv , S vh and S vv . For mono-static configurations, where the reciprocity theorem applies, S hv = S vh , only three components remain: S hh , S hv and S vv . In this case, the dimension of the target vector k becomes m = 3.

In the case of conventional two-pass InSAR system [START_REF] Massonnet | Radar interferometry, limits and potential[END_REF], only two channels are involved and m = 2. By denoting with c = ρ exp(jφ) the complex correlation coefficient, the target relative displacement d 12 between the two acquisitions can be retrieved from the exact knowledge of SAR antenna phase center positions, terrain height, acquisition geometry, and an estimate of the differential interferometric phase φ 12 . ρ 12 is called interferometric coherence and it describes both the local phase stability and the amplitude decorrelation of the InSAR pair. The phase information φ 12 allows phase differences (interferograms) to be computed in order to measure topography or target displacements between repeated pass acquisitions. In the general case, the m-dimensional interferometric target vector k will contain information about the relative displacements between each combination of 2 passes. The main parameter to be estimated is the speckle covariance matrix, from which normalized correlation coefficients can be easily derived.

The question if the SIRV models, despite their flexibility, are appropriate for describing any multidimensional SAR dataset, still remains. More exactly, "what is to be gained, and through which strategies, in pursuit of higher-order statistics, and non-Gaussian models"1 ? One recommendation has been released by the POLinSAR scientific community, namely to "Increase reliance on quantitative objective norms for comparing the performance of alternative analysis" 2 . This fact motivated us to propose a general framework which allows quantitative evaluation of fitting SIRV stochastic models with respect to a given multidimensional SAR dataset. Moreover, the illustrated multivariate high-resolution SAR datasets show that in this case, it can be worth challenging two specific and very important stochastic properties: the circularity and the sphericity.

This paper is organized as follows. In Section II we introduce classical circularity and sphericity tests along with their extensions to the SIRV stochastic model. Section III contains in detail elaborated method of quantitative assessment of the SIRV conformity, relying on testing spherical symmetry. In Section IV, the results obtained using the proposed robust tests are presented and analyzed with synthetic, multi-pass InSAR and POLSAR datasets. Section V presents some general remarks regarding the proposed methodological framework. In Section VI, some conclusions and perspectives are presented. Eventually, Appendix I introduces a zero-mean test, which is proposed for both Gaussian and SIRV stochastic processes.

In this paper, we assume that the covariance matrix [M ] is of full rank m. In the specific case of a non-invertible matrix, the same considerations can be applied using only the non-zero signal subspace.

II. CIRCULARITY AND SPHERICITY

In the context of multivariate SAR data, we introduce two stochastic properties, often assumed in case of a Gaussian model: circularity and sphericity. After adequately defining the properties, we present the classical tests (Gaussian model) and propose the ones extended to the SIRV stochastic model.

A. Circularity

A complex-valued random variable (r 0 = x 0 + iy 0 ) is circular [START_REF] Picinbono | On circularity[END_REF] if its distribution remains invariant to multiplication with complex numbers on the unity circle i.e. if the real random vector r 0 = (x 0 , y 0 ) T is spherically symmetric with respect to the origin. When dealing with the circularity of complex random vectors, we ought to rely on the second order properties [START_REF]Second-order complex random vectors and normal distributions[END_REF] by introducing the complex extended target vector ג = [k T , k H ] T . Starting from ,ג the second-order statistical properties of the complex target vector k can be analysed using the extended covariance matrix :

[R] = גג{E H } = [C] [P ] [P ] * [C] * ∈ C 2m , (1)
where [C] = E{kk H } is the complex positive semi-definite Hermitian symmetric covariance matrix and [P ] = E{kk T } is the complex symmetric pseudo-covariance matrix3 of the target vector k [START_REF] Neeser | Proper complex random processes with applications to information theory[END_REF].

Schreier et al. proposed the circularity Generalized Likelihood Ratio Test (GLRT) by employing the previously defined extended covariance matrix [START_REF] Schreier | A generalized likelihood ratio test for impropriety of complex signals[END_REF]. The considered hypotheses are:

H 0 : [P ] = 0, k is circular, H 1 : [P ] = 0, k is not circular . ( 2
)
This circularity test is checking for the block-diagonality of the extended covariance matrix with respect to the covariance and pseudo-covariance matrix.

1) Gaussian random processes: The probability density function (PDF) of a zero-mean complex circular Gaussian target vector k can be generalized with respect to its associated extended target vector ג as [START_REF]Second-order complex random vectors and normal distributions[END_REF][START_REF] Van Den Bos | The multivariate complex normal distribution -A generalization[END_REF]:

p G (k) = π -m (det[R]) -1 2 e -ג H [R] -1 ג 2 . (3)
The maximum-likelihood (ML) estimator of the extended covariance matrix is the Sample Extended Covariance Matrix (SECM) obtained by replacing the statistical mean from Eq. 1 with the spatial average:

[R] SECM = 1 N N i=1 ג i ג H i = [C] SCM [P ] SP M [P ] * SP M [C] * SCM , (4)
where N is the number of samples, [C] SCM and [P ] SP M are the Sample Covariance Matrix (SCM) and the Sample Pseudocovariance Matrix (SPM) estimators, respectively. Under the constraint [P ] = [0] m imposed by H 0 , the ML estimator is [START_REF] Schreier | A generalized likelihood ratio test for impropriety of complex signals[END_REF]:

[R] SECM0 = [C] SCM [0] m [0] * m [C] * SCM . (5)
By introducing Eqs. 4 and 5 to the Likelihood Ratio Test (LRT) associated to Eq. 2, Schreier et al. derived in [START_REF] Schreier | A generalized likelihood ratio test for impropriety of complex signals[END_REF] the Generalized Likelihood Ratio Test (GLRT) for N i.i.d. observed samples:

Λ(k 1 , ...k N ) = det [R] SECM (det [C] SCM ) 2 H0 ≷ H1 λ. (6)
The asymptotic distribution of the decision statistic under the null hypothesis is H 0 : -N ln Λ → χ 2 m(m+1) (Theorem 2 in [START_REF] Ollila | Adjusting the generalized likelihood ratio test of circularity robust to non-normality[END_REF]). Schreier et al. showed that the GLRT from Eq. 6 is invariant with respect to invertible linear transforms.

2) Spherically Invariant Random Vectors: The GLRT makes use of the Gaussian assumption when inserting [C] SCM and

[P ] SP M estimators in the LRT. Therefore, it ought to be modified before being applied on SIRV data.

By directly applying Tyler's Corollary 1 from [START_REF] Tyler | Robustness and efficiency properties of scatter matrices[END_REF], Ollila and Koivunen [START_REF] Ollila | Adjusting the generalized likelihood ratio test of circularity robust to non-normality[END_REF] showed that the GLRT ( [C] M L , [P ] M L ) for SIRV can be derived by dividing the logarithm of the GLRT statistics from Eq. 6 (ln Λ) by the correction factor γ, estimated as:

γ = m i=1 E{|ki| 4 } E{|ki| 2 } 2 r i (k i ) 2 + 2 . ( 7
)
with respect to the circularity coefficients, representing the "amount of circularity":

r i (k i ) = |E{k 2 i }| E{|k i | 2 } (8)
Following the same procedure as Ollila and Koivunen, N is replaced by (N -m) from the Box approximation of the GLRT [START_REF] Walden | On testing for impropriety of complex-valued gaussian vectors[END_REF]. Finally, according to Tyler's Corollary 1, the decision statistic under the null hypothesis is

H 0 : -(N -m) ln Λ γ → χ 2 m(m+1) .

B. Sphericity

The sphericity of a m-dimensional complex circular random vector z 0 is defined with respect to the augmented real random vector. Let ζ 0 = [x 0 T , y 0 T ] T be the augmented 2m × 1 real random vector, where x 0 and y 0 are its real and imaginary parts, respectively. The augmented real target vector can be obtained from the complex extended random vector using the following transform [START_REF]Second-order complex random vectors and normal distributions[END_REF][START_REF] Van Den Bos | The multivariate complex normal distribution -A generalization[END_REF]:

ζ = 1 2 [I] m [I] m -j[I] m j[I] m .ג (9)
The sphericity GLRT was introduced by Mauchly for real-valued m-dimensional Gaussian random processes [START_REF] Mauchly | Significance test for sphericity of a normal n-variate distribution[END_REF]. If the random process is circular, it can be directly extended to complex-valued random processes also. The hypotheses are:

H 0 : [C] = ς[I] m , k is spherical, H 1 : [C] = ς[I] m , k is not spherical , (10)
where ς ∈ R is unknown.

For the sake of simplicity and without loss of generality, the test is build such that [M ] = [I] m . As a consequence, the derived sphericity test must present invariance to linear transforms, also.

Given a set of N i.i.d. observed samples, the Mauchly's sphericity LRT is:

Λ s (k 1 , ...k N ) = m (det[M ]) 1 m tr[M ] H1 ≷ H0 λ s . (11)
In [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF], Λ s is reported to be invariant with respect to scale and invertible linear transforms (see Chapter 10.7 on page 431):

Λ s (k 1 , ...k N ) = Λ s ([V ]k 1 , ...[V ]k N ), (12)
with

[V ][M ][V ] H = [I] m .
One can observe that this linear transform is not the same as the whitening transform changing the SIRV covariance matrix to [I] m . Consequently, a SIRV exhibiting sphericity is a particular case of a SIRV with equal variances of the differences between all combinations of target vector components.

In particular, the sphericity property needs to hold when processing conventional multi-pass interferometric SAR data by interferometric pair. If the sphericity is not respected, fully multivariate processing should be applied (for example when jointly using the InSAR coherence and phase for 3D phase unwrapping [START_REF] Hooper | Phase unwrapping in three dimensions with application to InSAR time series[END_REF] or relevant scatterer detection [START_REF] Hooper | A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches[END_REF]): the InSAR channels are not independent with equal variance (homoscedastic).

1) Gaussian random processes: When replacing the covariance matrix by the SCM (ML estimate under Gaussian assumption), the GLRT is obtained as:

Λ s (k 1 , ...k N ) = m (det [M ] SCM ) 1 m tr [M ] SCM H1 ≷ H0 λ s . (13)
Asymptotically, H 0 : -N ln Λ s → χ 2 (m+2)(m-1) [START_REF] Mauchly | Significance test for sphericity of a normal n-variate distribution[END_REF][START_REF] Muirhead | Asymptotic distributions in canonical correlation analysis and other multivariate procedures for nonnormal populations[END_REF].

2) Spherically Invariant Random Vectors: Following Muirhead and Waternaux original studies on the robustness of the GLRT from Eq. 11 when sampling from a SIRV [START_REF] Muirhead | Asymptotic distributions in canonical correlation analysis and other multivariate procedures for nonnormal populations[END_REF], Tyler proposed two different robust approximated GLRTs for sphericity [START_REF] Tyler | Robustness and efficiency properties of scatter matrices[END_REF][START_REF] Tyler | Radial estimates and the test for sphericity[END_REF].

The first method (Corollary 1 from [START_REF] Tyler | Robustness and efficiency properties of scatter matrices[END_REF]) is derived using the same reasoning as in the case of circularity. The derived approximated GLRT is reported to be inefficient for moderate departures from Gaussianity.

In this paper, we adopt the second method for constructing the approximated GLRT for sphericity (Corollary 4 from [START_REF] Tyler | Robustness and efficiency properties of scatter matrices[END_REF]):

a Λ s (k 1 , ...k N ) = (det [M ] F P ) 1 m < H0 λ s , (14)
where [M ] F P is the iterative Fixed-Point covariance matrix estimator (already employed with multivariate SAR data [START_REF] Vasile | Coherency matrix estimation of heterogeneous clutter in high resolution polarimetric SAR images[END_REF]):

[ M ] F P = m N N i=1 k i k H i k H i [ M ] -1 F P k i . ( 15
)
When originally introducing this estimator, Tyler showed in [START_REF]A distribution-free M-estimator of multivariate scatter[END_REF] that [M ] F P is an affine-invariant covariance matrix Mestimator [START_REF] Maronna | Robust M-estimators of multivariate location and scatter[END_REF]. It has been demonstrated by Tyler that the Fixed-Point estimator is an AML estimator for SIRVs [START_REF] Tyler | Statistical analysis for the angular central Gaussian distribution[END_REF]. Hence, all conditions required to apply Corollary 4 from [START_REF] Tyler | Robustness and efficiency properties of scatter matrices[END_REF] are now fulfilled and H 0 : -N m m+1 ln a Λ s → χ 2 (m+2)(m-1) in distribution. Due to the Theorem V.4 in [START_REF] Pascal | Performance analysis of covariance matrix estimates in impulsive noise[END_REF], N is replaced by N m m+1 .

Finally, the formerly introduced circularity test complements the sphericity test presented here: the test for circularity checks that there is no correlation between the real and the imaginary part of the complex random vector, while the sphericity test checks for equal independence between the target vector components (independence and homoscedasticity). Hence, circularity is a requirement before testing for sphericity, as k and k * cannot be independent.

III. SPHERICAL SYMMETRY

According to Vershik's definition [START_REF] Vershik | Some characteristic properties of Gaussian stochastic processes[END_REF], z 0 is spherically invariant if and only if the characteristic function of the augmented random vector ζ 0 can be written as:

Φ ζ0 (ζ 0 ) = φ ζ 0 T [C]ζ 0 2 , φ(v) = ∞ 0 e -τ v p τ (τ )dτ, (16)
where [C] is a positive definite characteristic matrix and p τ (τ ) is a probability measure on [0, ∞). More intuitively, by applying

ζ 0 = [C] -1/2
ζ 0 linear transform on the augmented random vector, the density generator function φ(v) in Eq. 16 takes the form:

φ(ζ 0 T ζ 0 /2) = f (||ζ 0 || 2 ), (17)
where ||.|| is the Euclidian norm. This implies that the SIRV ζ 0 is reduced to a new SIRV ζ 0 with its new covariance matrix equal to [I] m (identity matrix of order m). ζ 0 clearly exhibits spherical invariance.

Let now φ be defined in Eq. 16 so that the matrix [C] is the speckle normalized covariance matrix [START_REF] Vasile | Optimal parameter estimation in heterogeneous clutter for high resolution polarimetric SAR data[END_REF] and the conditions from the Yao's representation theorem are respected [START_REF] Yao | A representation theorem and its applications to spherically-invariant random processes[END_REF]. In consequence, for a given texture pdf p τ (τ ), the SIRV process describing the complex target vector k can be defined as :

p k (k) = 1 π m det[M ] h m k † [M ] -1 k , (18)
with

h m (q) = ∞ 0 τ -m exp -q τ p τ (τ )dτ . det[M ] denotes the determinant of the matrix [M ]. Yao (Lemma 4.1 in [2]
) also demonstrated SIRV closure under invertible linear transform (e.g. as suggested for Eq. 17). All target vectors k satisfying these conditions accept a product model stochastic representation [START_REF] Fang | Symmetric Multivariate and Related Distributions[END_REF].

Note also that a random target vector k E has an elliptically symmetric distribution [START_REF] Ollila | Complex elliptically symmetric distributions: survey, new results and applications[END_REF] if it is affinely equivalent in distribution to a spherically symmetric target vector k S :

k E = [A]k S + b. (19)
Therefore, spherical symmetry represents a particular case of elliptical symmetry when Hence, on one side, a spherically symmetric white SIRV vector, being characterized with unity covariance matrix, is necessarily spherical. On the other side, the previously introduced sphericity test is no longer valid in the absence of the elliptical symmetry, making the last property a mandatory constraint in testing the sphericity.

The spherical symmetry property in multivariate statistics is defined with respect to a SIRV with covariance matrix [I] m (white SIRV). The general SIRV case, with unknown covariance matrix, was studied under the "elliptical symmetry" property.

Since the family of spherical symmetric distributions -SSD can be considered as the standardized form of the family of elliptically contoured distributions -ECD (by employing the whitening transform from Eq. 17), we adopt in this paper the original Yao's notation by considering zero mean elliptically contoured distributions as SIRVs. Hence, we use the term "spherical symmetry" for both SSD and ECD goodness-of-fit testing. See [START_REF] Fang | Testing spherical and elliptical symmetry[END_REF] for an exhaustive presentation on SSD and ECD.

Spherical symmetry testing was firstly introduced by Kariya and Eaton [START_REF] Kariya | Robust tests for spherical symmetry[END_REF], using an alternative form of the Lehmann and

Stein lemma [START_REF] Lehmann | On the theory of some non-parametric hypotheses[END_REF] with known covariance matrix. In the common situation where the covariance matrix is estimated from data, several strategies for robust spherical symmetry testing have been proposed by Beran [START_REF] Beran | Testing for ellipsoidal symmetry of a multivariate density[END_REF], King [START_REF] King | Robust tests for spherical symmetry and their application to least squares regression[END_REF], Baringhaus [START_REF] Baringhaus | Testing for spherical symmetry of a multivariate distribution[END_REF], Fang et al. [START_REF] Fang | A necessary test of goodness of fit for sphericity[END_REF], Manzotti et al. [START_REF] Manzotti | A statistic for testing the null hypothesis of elliptical symmetry[END_REF] and Huffer et al. [START_REF] Huffer | A test for elliptical symmetry[END_REF] among the most recent publications. Li et al. proposed in [START_REF] Li | Some Q-Q probability plots to test spherical and elliptical symmetry[END_REF] a graphical method for spherical symmetry testing: the Q-Q probability plots. This method has been applied in [START_REF] Niu | Tests for the elliptical symmetry of hyperspectral imaging data[END_REF] to hyperspectral image analysis.

According to [START_REF] Huffer | A test for elliptical symmetry[END_REF], one of the most powerful spherical symmetry tests was proposed in [START_REF] Schott | Testing for elliptical symmetry in covariance-matrix-based analyses[END_REF] for real random vectors. In this section, we have adapted the Schott test for multidimensional complex SAR data analysis.

A. The Schott test for circular complex random vectors

Assuming the existence of the forth order moment (or quadricovariance) matrix:

[M ] 4 = k i k H i ⊗ k i k H i , (20)
Schott proposed the Wald test [START_REF] Moore | Generalized inverses, Wald's method, and the construction of chi-squared tests of fit[END_REF] for verifying that the structure of [M ] 4 corresponds to a SIRV (as originally given by Tyler in [START_REF] Tyler | Robustness and efficiency properties of scatter matrices[END_REF]). This structure holds for circular complex random vectors, also [START_REF] Ollila | Complex elliptically symmetric distributions: survey, new results and applications[END_REF].

With complex random vectors, the first modification is the sample quadricovariance estimator. According to [START_REF] Ferrari | Estimation of the frequencies of a complex sinusoidal noisy signal using fourth order statistics[END_REF], the sample complex quadricovariance estimator can be expressed in terms of the Kronecker product ⊗ as:

[M ] 4 = 1 N N i=1 k i k H i ⊗ k i k H i , (21)
where the transposed operator T is replaced by the conjugate and transpose operator H . Its corresponding standardized form is:

[M ] 4 * = [M ] -1 2 H ⊗ [M ] -1 2 H [M ] 4 [M ] -1 2 ⊗ [M ] -1 2 . ( 22
)
According to the Schott's theorem, the Wald test statistic for spherical symmetry can be expressed as:

T Schott = 2N β 1 tr [M ] 2 4 * + β 2 vec [I] m H [M ] 2 4 * vec [I] m -[3β 1 + (m + 2)β 2 ] m(m + 2)(1 + κ) 2 , (23)
where

β 1 = (1 + θ) -1 /24, (24)
β 2 = -3a[24(1 + θ) 2 + 12(m + 4)a(1 + θ)] -1 , (25)
a = (1 + θ) + (1 + κ) 3 -2(1 + κ)(1 + η), (26)
with the Mardia's kurtosis κ and the generalized higher order scalar moments θ, η given by:

(1 + κ) = 1 m(m + 2)N N i=1 k H i [M ] -1 k i 2 , (27)
(1 + θ) = 1 m(m + 2)(m + 4)N N i=1 k H i [M ] -1 k i 3 , (28)
(1 + η) = 1 m(m + 2)(m + 4)(m + 6)N N i=1 k H i [M ] -1 k i 4 . (29)
Asymptotically,

T Schott → χ 2 vm complex with v m complex = m 2 (m+1)(m+5)

12

-1. This represents the second modification with respect to the Wald test from [START_REF] Schott | Testing for elliptical symmetry in covariance-matrix-based analyses[END_REF]. According to Schott, the degrees of freedom is set according to the number of unknowns of the quadricovariance matrix:

v m complex = 2(v m real + 1) - m(m -1) 2 + m -1, (30)
where v m real = m 2 + m(m-1)(m 2 +7m-6)

24

-1 as in [START_REF] Schott | Testing for elliptical symmetry in covariance-matrix-based analyses[END_REF]. The second term in Eq. 30 comes from the number of real elements

of [M ] 4 : m(m-1)
2 elements of the form x 2 i x 2 j and m elements of the form x 4 i . It has been proven in [START_REF] Engle | Likelihood Ratio, and Lagrange Multiplier Tests in Econometrics[END_REF] that the Wald test and the LRT are asymptotically equivalent.

IV. RESULTS AND DISCUSSIONS

To illustrate the proposed tests, results obtained with the synthetic POLSAR dataset, high-resolution TerraSAR-X multi-pass InSAR and very high-resolution ONERA RAMSES [START_REF] Dreuillet | The ONERA RAMSES SAR: latest significant results and future developments[END_REF] POLSAR data are reported. These data sets are shown in Fig. 1.

[Figure 1

about here.]

This section is dedicated to the analysis of the three data sets in terms of circularity, sphericity and spherical symmetry. Since the mean equal to zero is a requirement for all the derived tests, we equally show results obtained by applying the zero-mean test, introduced in Appendix I, on the augmented real random vector ζ 0 . This test is valid for both Gaussian and SIRV stochastic models.

Finally, the discussions associated to the presented results concern these high-resolution multivariate SAR data sets, only.

The estimation neighbourhood is the 13 × 13 boxcar and the false alarm probability threshold is p f a = 0.01.

A. Synthetic data

The synthetic dataset is composed out of nine different regions (Fig. 2a). Six of them are characterized with the SIRV PolSAR clutter, while for the remaining three the clutter is Gaussian. In four regions the additive thermal noise (circular or non-circular) is present. As well, we introduce a coherent scattering through the simulated elementary reflectors (two trihedrals, dihedral and dipole). In the derivation we relied on the ML deterministic texture estimator:

τ F P = k H [ M ] -1 F P k m . ( 31
)
[Figure 2 about here.]

Aside from the assumed detection of the coherent scattering sub-regions (elementary reflectors), circularity test rejects successfully the regions corrupted with additive non-circular thermal noise. The heterogeneity (borders between different regions)

does not appear to influence this test significantly.

Sphericity test detects properly the Gaussian clutter with additive circular thermal noise as well as the Gaussian clutter without noise. As it is the case with the circularity test, the influence of non-stationarity is negligible.

However, the spherical symmetry test does not seem to be immune to the heterogeneity, but still quite successfully rejects the coherent scattering (deterministic target).

B. High-resolution multi-pass InSAR data

The 3-pass interferometric stripmap HH images were acquired in 2009, at 11-day interval, over the Argentire village, France with a mean incidence angle of 5 0 , an azimuth-resolution of 3.3 m and a slant-range resolution of 1.8 m. This data set has been used for SAR tomography over this area as it exhibits a high coherence level over the main buildings from [Figure 3 about here.]

[Table 1 about here.]

The pixels illustrated in magenta on Fig. 3-(a) indicate where the zero-mean test from Appendix I is rejected. It can be observed from Table I that this dataset is zero-mean: the percentage of rejected pixels is much less than the significance level a priori set (1% in all cases). Fig. 3-(b) shows in red the pixels where the adjusted circularity test is rejected. These pixels should be processed as 2m real random vectors. The percentage of rejected non circular pixels (cf. Table I) is, although larger than the significance level, still small enough (3.14%) for us to conclude that this dataset is circular.

Fig. 3-(c) and Fig. 3-(d) illustrate in blue the pixels where the sphericity test from Eq. 14 is accepted, and in green the pixels where the spherical symmetry test from Eq. 23 is rejected. The quantitative results summed up in Table I indicate that both sphericity and spherical symmetry properties are significant for this dataset. Qualitatively, it can be observed that:

• localized mainly in regions with high density of strong scatterers, the nonspherical pixels should be treated as a fully multivariate process inside the local neighborhood when estimating the InSAR coherence and phase parameters;

• mostly located in the same areas of the image, the nonspherically symmetric pixels indicate where the SIRV (and consequently the compound Gaussian) model fails to properly describe the multivariate clutter.

This specific behavior may be linked with the presence of strong deterministic scattering in urban areas: the analyzed target is not distributed but deterministic. Nevertheless, the proposed tests build a methodological framework to study this effect with respect to the SIRV model.

Spherical symmetry is a prerequisite for sphericity. For quantitatively validating the obtained results, the percentage pixels detected as spherical and rejected for spherical symmetry was computed in Table I: the value is quite small considering the significance level. This condition is met for non circular and non-zero mean pixels, also.

Additionally, about 18% of non circular pixels are rejected by the Schott test for spherical symmetric circular complex random vectors. This illustrates that circularity should be tested before testing for spherical symmetry.

C. Very high resolution POLSAR data

Illustrated in Fig. 1-(b), this data set was acquired over Toulouse, France with a mean incidence angle of 50 0 . It represents a fully polarimetric (monostatic mode) X-band acquisition with a spatial resolution of approximately 0.5 m in range and azimuth.

[Figure 4 about here.]

[Table 2 about here.]

The background image from Fig.

is the fixed point Polarimetric Whitening Filter (PWF) span obtained by applying Eq. 31 with polarimetric target vectors. We have used here the same representation and the same color coding as in the previous section, while Table II sums up the POLSAR results. It can be observed that this dataset is globally zero-mean and not spherical: the corresponding percentages are less than the imposed significance level. The SIRV model (spherical symmetry) holds in about three quarters of all realizations, while the non sphericity confirms that multivariate statistical modeling is correctly employed since the clutter is not spherical.

In other words, Fig. 4-(d) underlines the same phenomena as in Fig. 3-(d): the spherical symmetry is rejected over the urban areas exhibiting strong deterministic scattering. Quantitatively, one can notice in Table II that the Schott test for spherical symmetric circular complex random vectors is rejecting non circular pixels in about 31% of all cases, only. Hence circularity testing is mandatory prior to testing the SIRV model conformity.

Finally, both Fig. 4-(c) and Table II indicate a relatively high percentage of noncircular pixels in this POLSAR dataset. These pixels are mainly located in the weak backscattering image areas (shadowing). Since the data are zero-mean, this effect is not induced by any non-centered thermal noise additive component. However, it may be introduced by the specific calibration of this airborne very high resolution POLSAR data: the motion compensation module tracks and rectifies the signature of specific calibrated point targets on the ground.

V. GENERAL REMARKS

When dealing with multivariate high-resolution SAR data, it is crucial to decide if a specific stochastic model is properly fitting the experimental dataset inside the estimation neighborhood. The accepted stochastic process with the smallest number of parameters should be selected. In other words, if both the SIRV model and the Gaussian model are fitting the data, the latter will have better estimation performances with a finite number of samples. Based on the results presented in this paper, we provide the methodological framework to asses multivariate SAR data conformity:

1) Zero-mean test from Appendix I. STOP if the test is rejected: the mean parameter should be included in the SIRV stochastic model.

2) Circularity test: STOP if the test is rejected: the augmented (or extended) random vectors should be used instead of the complex target vectors.

3) Sphericity test: STOP if the test is accepted: it is better to use one complex circular random variable for each channel, separately.

4) Spherical symmetry test: STOP if the test is rejected: SIRVs are not stationary inside the neighborhood and improved statistical (or simply deterministic) modeling should be addressed.

5) Testing for non-Gaussianity (for example the one proposed in [START_REF] Vasile | Optimal parameter estimation in heterogeneous clutter for high resolution polarimetric SAR data[END_REF]). STOP if the test is rejected: the dataset is more likely to be locally Gaussian and multivariate zero-mean complex circular Gaussian statistics should be employed.

By successively performing the proposed tests for a specific multivariate SAR dataset, it is possible to asymptotically evaluate the pertinence of various model-based statistical processing schemes (filtering, segmentation or detection).

In terms of theoretical performance analysis, the adjusted generalized LRT is asymptotically uniformly most powerful according to the Neyman-Pearson lemma. This "optimality" holds provided the estimators plugged into the LRT (or the Wald test) are consistent and unbiased, which is the case for our study.

Special care must be taken when applying the tests with multivariate SAR data. In theory, the proposed conformity testing holds as long as the observed number of samples (estimation neighborhood) is large enough with respect to the dimension of the target vector, especially for the spherical symmetry test (based on the specific structure of the quadricovariance).

Finally, it is important to stress that no predefined analytical form was imposed on the texture probability function when establishing the conformity tests. Therefore, they can be directly applied for a wide class of stochastic processes currently used for describing multivariate high-resolution SAR data.

VI. CONCLUSIONS

This paper has presented a new methodological framework to asymptotically asses the conformity of multivariate highresolution SAR data. The proposed approach consist of applying successively three statistical hypotheses tests for verifying three important statistical properties: circularity, non sphericity and spherical symmetry. The latter is asymptotically equivalent, under certain hypotheses, to the conformity of the experimental data with respect to the SIRV product model. In addition, the zero-mean (from Appendix I) and the non-Gausiannity [START_REF] Vasile | Optimal parameter estimation in heterogeneous clutter for high resolution polarimetric SAR data[END_REF] tests can be used to decide which model is better suited to asymptotically fit the experimental data.

The proposed framework, aside from the fact that it is gathering the most notable advances in the field of signal processing, is introducing the extension of both the sphericity and the spherical symmetry tests with zero-mean complex circular SIRV assumption.

The effectiveness of the proposed detection schemes was illustrated by high-resolution TerraSAR-X multi-pass InSAR and very high resolution ONERA RAMSES POLSAR data. The conclusions driven from the analysis of the obtained results are important with respect to the two tested datasets: non-sphericity can be an important issue for spaceborne multi-pass InSAR, while circularity is important for airborne POLSAR data.

It has been illustrated that in strong heterogeneous clutter, such as the urban environment, the SIRV model can fail. The bottom line is that characterization of urban regions is much more complex (and difficult) -since a more complex model, with more parameters, may be required. In the light of the results shown in this paper, SIRV models may be less appropriate for urban areas characterization. However, alternative explanations are possible. As an example, the root of this inappropriateness might as well be the assumed ergodicity / stationarity (in spatial sense) for the backscattered signal and, also, even in the hypothesis of randomness: targets exhibit a deterministic behavior.

First, the very use of a sliding analysis window for estimating the stochastic parameters of the scattered signal may be questioned, as it implicitly assume that the considered signal is ergodic / stationary. While this hypothesis holds for distributed and uniform targets, where the physical parameters (and, thus, the electromagnetic scattering behavior) differs very little from one resolution cell to another, in urban areas the physical structure (and, as such, its electromagnetic behavior) may change considerably from one resolution cell to the next. This makes the hypothesis of ergodicity / stationarity less applicable.

Second, one should note that even the randomness of the radar echo is not given, but assumed. This is mainly a way to deal with the inherent complexity of the signal. Anyway, for identical measuring conditions, the recorded radar data is perfectly identical. Even if small differences in measuring conditions lead to strong discrepancies in the recorded data, this is not an evidence for randomness, as such behavior can be fully explain under a deterministic paradigm -the chaotic models. Various parameters, such as meteorological conditions and, even more important, the changes that the target suffers in time (between two succeeding acquisitions, for example), account for the observed randomness of the recorded data. However, these changes of the target are more significant for green targets (such as forests and agricultural fields), where humidity and wind modify both their physical structure and their electromagnetic behavior. On the other hand, those changes are less significant for urban targets and, as such, randomness is less likely for the latter.

In perspective, applying chaotic (or pseudo-chaotic) [START_REF] Gao | Direct dynamical test for deterministic chaos and optimal embedding of a chaotic time series[END_REF][START_REF] Gao | Power-law sensitivity to initial conditions in a time series with applications to epileptic seizure detection[END_REF][START_REF] Haykin | Chaotic dynamics of sea clutter[END_REF] models to POLSAR / InSAR data from urban areas can be a possible solution. These models should be able to take into account the deterministic features of those areas (presence of dihedral angles, straight edges, cavities, etc.), while still leaving room for some unpredictability (orientation of those elements).

Using chaotic models in POLSAR and multi-pass InSAR data will make the object of our future work.

APPENDIX I: ZERO-MEAN TEST FOR BOTH GAUSSIAN AND SIRV STOCHASTIC MODELS

Under circularity, the zero-mean testing can be performed by using the T 2 -statistic. When testing the hypothesis H 0 that a mean vector k is equal to zero, the generalized likelihood ratio criterion for the circular multivariate Gaussian model is:

T 2 = N kH [ M ] -1 SCM k, (32)
with T 2 → χ 2 m under H 0 (see Theorem 5.2.3 from [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF]). In case circularity is not assumed, this test can be directly applied on the augmented real random vector ζ 0 . This results holds asymptotically for the SIRV models according to Theorem 5.7.1 from [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF]. More details can be found in [START_REF] Anderson | On the theory of multivariate elliptically contoured distributions and their applications[END_REF][START_REF]Maximum likelihood estimators and likelihood ratio criteria for multivariate elliptically contoured distributions[END_REF][START_REF] Anderson | Maximum-likelihood estimates and likelihood-ratio criteria for multivariate elliptically contoured distributions[END_REF].

The GLRT from Eq. 32 is not scale invariant like the ones from Sections II-A.2 and II-B.2. Thus, Tyler's Corollaries 1 and 4 from [START_REF] Tyler | Robustness and efficiency properties of scatter matrices[END_REF] cannot be applied directly.

  [A] = [I] m ,b = [0] m×1 . In our case, by identifying in the SIRV whitening transform from Eq. 17, it yields [A] = [C] -1/2 and b = [0] m×1 (due to the zero mean).

Fig. 1 -

 1 (a)[START_REF] Ilie | Spaceborne SAR tomography: application in urban environment[END_REF] . The background image from Fig. 3-(a), (b), (c), (d) is the fixed point Interferometric Whitening Filter (IWF) span (Eq. 31).

Fig. 1 :Fig. 2 :Fig. 3 :Fig. 4 :

 1234 Fig. 1: Data sets: (a) Synthetic POLSAR data, intensity color composition of the target vector elements k 1 -k 3 -k 2 in Pauli basis, (b) TerraSAR-X 3-pass InSAR data, amplitude color composition of the complex random vector elements k 1 -k 2 -k 3 , (c) Toulouse, RAMSES POLSAR data, intensity color composition of the target vector elements k 1 -k 3 -k 2 in Pauli basis. figure

TABLE I :

 I Argèntiere, TerraSAR-X 3-pass InSAR data, X-band: detection results (percentage computed with respect to the total number of pixels).

	table
		Non-zero-mean Non-circularity	Sphericity	Non-spherical-symmetric
		(3751 pixels)	(32882 pixels) (424856 pixels)	(328125 pixels)
	Non-zero-mean	0.36 %	0.01%	0.02%	0.18%
	Non-circularity	0.01%	3.14%	1.27%	0.55%
	Sphericity	0.02%	1.23%	40.52 %	3.14%
	Non-spherical-symmetric	0.18%	0.6%	3.14%	31.3%

TABLE II :

 II Toulouse, RAMSES POLSAR data, X-band: detection results (percentage computed with respect to the total number of pixels).

	table
		Non-zero-mean Non-circularity	Sphericity	Non-spherical-symmetric
		(45249 pixels) (230104 pixels) (1283 pixels)	(373983 pixels)
	Non-zero-mean	4.52%	0.91%	0%	2.3%
	Non-circularity	0.91%	23.01%	0.02%	7.2%
	Sphericity	0%	0.02%	0.13%	0%
	Non-spherical-symmetric	2.3%	7.18%	0%	37.4%

POLinSAR

Workshop, ESA ESRIN, Frascati, Italy: seed question no. 3 for the "Methods & Theoretical Modelling" session. 2 POLinSAR 2013 Workshop, ESA ESRIN, Frascati, Italy: see "Methods & Theoretical Modelling" session summary at (https : //earth.esa.int/web/guest/polinsar-2013).

Also known as "relation matrix" or "complementary covariance matrix".

ACKNOWLEDGMENT

The authors would also like to thank Dr. C. Tison (CNES, France) for providing the high-resolution POLSAR images over Toulouse and the German Aerospace Center (DLR) for providing the TerraSAR-X stripmap SAR images through the MTH0232 and the MTH0828 projects.


Nikola Besic
email: nikola.besic@gipsa-lab.grenoble-inp.fr

Gabriel Vasile

Jocelyn Chanussot

Srdjan Stankovic
email: srdjan@ac.me

Jean-Pierre Dedieu

Stochastically Based Wet Snow Detection with Multitemporal SAR Data

This paper proposes a new method for wet snow detection using multitemporal SAR data. The proposed change detection method is primarily based on the comparison between two X band SAR images acquired during the accumulation (winter) and the melting (spring) seasons, in the French Alps. The new membership decision method is relying on the local intensity statistics of the SAR images by considering the backscattering ratio as a stochastic process: the probability that "the intensity ratio is smaller than the corresponding locally varying dry/wet snow threshold" is larger than a predefined confidence level.

The proposed snow backscattering simulations indicate a far more complex relation between the backscattering properties of the two snow types, with respect to the conventional assumption of the augmented electromagnetic absorption associated to the wet snow. The most conclusive results are confirmed at X band by in situ measurements. We show that this complexity is mostly caused by the dominance of a different backscattering component for each type of snow, leading to a significant angular sensitivity. Therefore, we introduce a variable threshold matrix instead of an unique threshold: the threshold is expressed as a function of the local incidence angle (LIA) for each pixel. The proposed snow backscattering simulator (SBS) allows deriving the locally varying threshold as the ratio between dry and wet snow backscattering coefficients computed for the SAR system parameters (polarization, frequency), the local topography parameters (LIA) and for the a priori derived and spatially averaged target parameters (underlying ground and snow cover surface properties, snow grain size, underlying ground dielectric permittivity, snow wetness). The threshold varying with respect to the frequency and the choice of the winter season reference image, allow the proposed algorithm applicability in wider frequency range.

I. INTRODUCTION

Snow cover represents a significant water resource in mountainous regions such as the French Alpes. The exact amount of liquid water contained in the snow pack is given through the physical quantity -Snow Water Equivalent (SWE). Due to numerous difficulties in performing "in situ" measurements (high altitudes, global coverage), the estimation of this quantity turns out to be suitable for remote sensing applications. Generally, the SWE is a function of two independent snow physical properties: density and depth [START_REF] Shi | Estimation of snow water equivalence using sir-c/x-sar, part i: Inferring snow density and subsurface properties[END_REF]. The significant snow pack penetration, along with the high spatial resolution and the clouds penetration capability, candidates Synthetic Aperture Radar (SAR) remote sensing for this kind of application [START_REF] Ulaby | Microwave response of snow[END_REF], [START_REF] Rott | Possibilities and limits of synthetic aperture radar for snow and glacier surveying[END_REF].

Depending on the liquid water content, snow can be considered as dry or wet. Dry snow is defined as snow consisted just of ice crystals embedded in air, without any liquid water. Its presence is characteristic for temperatures below 0 • C. Wet snow can be found during the melting season, when the temperature exceeds 0 • C. It contains a certain amount of liquid water introduced quantitatively through the wetness (volumetric liquid water content) [START_REF] Rees | Remote Sensing of Snow and Ice[END_REF]. These two snow types behave like two completely different materials with respect to their dielectric properties: the wet snow contains liquid water with a dielectric constant differing significantly from the one of the ice. This fact is pointing out to a difference in backscattering mechanisms [START_REF] Ulaby | Snowcover influence on backscattering from terrain[END_REF], and consequently that different feature extraction methods should be applied. Therefore, the estimation of any snow pack parameter defining the SWE, by means of SAR remote sensing, requires firstly the proper identification of the snow cover type.

Although studies on snow mapping by polarimetric SAR existed already [START_REF] Rott | Thematic studies in alpine areas by means of polarimetric sar and optical imagery[END_REF], it was not before the ERS-1 started providing repeat pass images that the idea of snow mapping based on SAR multitemporal data appeared [START_REF] Rott | Capabilities of ers-1 sar for snow and glacier monitoring in alpine areas[END_REF]. Further refinement on the initial change detection method leaded eventually to the compact algorithm based on the ratio of two C band SAR images introduced in [START_REF] Nagler | Retrieval of wet snow by means of multitemporal sar data[END_REF]. The Nagler and Rott method requires as inputs the SAR wet snow image and the reference SAR image of the dry snow (or the snow free terrain). After corregistration, either mutilooking or speckle filtering is applied on the two SAR intensities before constructing the backscattering ratio image. This ratio image is then georeferenced and an unique threshold of -3 dB is used to discriminate the wet snow from other surfaces. The expected difference in backscattering is justified by the increased electromagnetic absorption of the wet snow. The resulting maps are successfully validated by using snow terrain optical images.

We propose in this article a novel snow change detection method with X-band SAR data for wet snow detection [START_REF] Besic | Stochastically based wet snow mapping with sar data[END_REF]. Our goal is to introduce an algorithm which can be applied to a wider frequency range and is eventually more suited to the presence of the speckle noise. The first novelty consists in the restriction in choosing the reference image. Due to the demonstrated difference in backscattering between the snow free terrain and the dry snow, the latter is considered as reference.

Single snow layer backscattering is simulated by applying fundamental scattering theories: Integral Equation Model, Quasi Crystalline Approximation and Quasi Crystalline Approximation with Coherent Potentials. Qualitative conclusions concerning the interaction between the two snow types and the electromagnetic waves are drawn [START_REF] Besic | Dry snow backscattering sensitivity on density change for swe estimation[END_REF]. Additionally, some conclusions referring to the wet snow backscattering are confirmed through in situ measurements, validating the proposed simulator. Consequently, it is shown that the ratio between the wet and dry snow backscattering coefficients is significantly sensitive to the value of the local incidence angle (LIA), therefore requiring the use of a variable threshold (function of the LIA).

Further, the speckle noise statistics [START_REF] Vasile | High resolution SAR interferometry: influence of local topography in the context of glacier monitoring[END_REF], [START_REF] Julea | Simulation of SAR images and radar coding of georeferenced information for temperate glacier monitoring[END_REF], [START_REF] Vasile | Optimal parameter estimation in heterogeneous clutter for high resolution polarimetric SAR data[END_REF], [START_REF] Anghel | Short-range wideband FMCW radar for millimetric displacement measurements[END_REF], [START_REF] Besic | Stochastic approach in wet snow detection using multitemporal SAR data[END_REF] is introduced through the local estimation of the backscattering coefficients ratio probability. This allows additional enhancement of the discrimination accuracy, which is illustrated through the matching of the independently obtained HH and VV maps. The validation is performed by comparing snow maps obtained using TerraSAR-X data with interpolated temperature map.

The paper is organized as follows. In Section II, we briefly introduce the proposed snow backscattering simulator (SBS) followed by the simulation results and by the comparison with in situ ground based radar measurements. This section contains the qualitative conclusions concerning dry and wet snow backscattering. Section III is the core of the article, containing the description of the wet snow detection algorithm. In Section IV, we present the results obtained with two TerraSAR-X dual-pol stripmap images acquired in the French Alps. This section as well provides the method validation, while Section VI concludes this article.

II. SNOW BACKSCATTERING ANALYSIS

On one side, dry snow is characterized by a dielectric constant depending on the ice volume fraction, only. In this case, the dielectric constant is real, indicating the absence of medium dielectric losses [START_REF] Rees | Remote Sensing of Snow and Ice[END_REF], [START_REF] Hallikainen | Dielectric properties of snow in the 3 to 37 ghz range[END_REF]. One can notice that a general consensus

has not yet been achieved in the scientific community: the work of Matzler et al. [START_REF] Matzler | Solar System Ices[END_REF], [START_REF] Matzler | Dielectric properties of freshwater ice at microwave frequencies[END_REF] may indicate a complex dielectric constant. We have adopted the former assumption [START_REF] Rees | Remote Sensing of Snow and Ice[END_REF], [START_REF] Hallikainen | Dielectric properties of snow in the 3 to 37 ghz range[END_REF] by neglecting the imaginary part.

On the other side, the wet snow dielectric constant is a complex quantity with the imaginary part reflecting the presence of medium absorption losses. It is also function of the ice volume fraction, as well as of the wetness and the operating frequency [START_REF] Rees | Remote Sensing of Snow and Ice[END_REF], [START_REF] Hallikainen | Microwave Remote Sensing of Sea Ice[END_REF].

In terms of scattering theory, the extinction (absorption and scattering) properties of the two snow volume types require different approaches in volume modeling. The model we use for dry snow is based on the multiple scattering assumption: Quasi Crystalline Approximation (QCA) through the Dense Media Radiative Transfer (DMRT), taking thus into account the scattered waves influence on the other particles. Furthermore, it assumes waves coherence and it stays valid even in the case of a media more dense than the dry snow [START_REF] Tsang | Scattering of Electromagnetic Waves[END_REF], [START_REF] Tsang | Scattering of Electromagnetic Waves[END_REF]. Note that the extinction losses are reduced to the scattering phenomena only by neglecting the electromagnetic absorption of the ice particles (real ice dielectric constant). We assume the effective shape of the particles to be spherical.

The wet snow medium reveals a significant spatial variation of dielectric permittivity [START_REF] Tsang | Scattering of Electromagnetic Waves[END_REF], [START_REF] Tsang | Scattering of Electromagnetic Waves[END_REF]. Such an important difference in dielectric properties (between the host medium and the scattering particles) requires applying a modified QCA theory. Hence the selected QCA with Coherent Potentials (QCA-CP) assumes that the coherent wave interacts with the new average medium rather than the host one. The water inclusions are introduced by means of mixed ice-water spherical particles [START_REF] Longepe | Snowpack characterization in mountainous regions using c-band sar data and a meteorological model[END_REF] with a radius determined by the wetness. Due to non negligible wet snow absorption we introduce the albedo coefficient, also [START_REF] Longepe | Snow remote sensing using spaceborne sar imagery at l-and c-bands[END_REF].

[Figure 1 about here.] The total backscattering coefficient σ 0 of a snow pack (distributed target) is expressed as the sum of four backscattering coefficients, corresponding to the following components [START_REF] Fung | Microwave Scattering and Emission Models for Users[END_REF], [START_REF] Shi | Estimation of snow water equivalence using sir-c/x-sar, part ii: Inferring snow depth and particle size[END_REF] (Fig. 1):

1) Snow pack surface component (σ s ), 2) Underlying ground surface component (σ gr ), 3) Snow volume component (σ v ), 4) Ground-volume interaction component (σ gv ),

σ 0 = σ s + σ gr + σ v + σ gv . ( 1
)
Surface backscattering is simulated using the Integral Equation Model (IEM-B) by assuming an isotropically exponentially correlated surface (introduced through the surface spectra) [START_REF] Fung | Microwave Scattering and Emission Models for Users[END_REF].

The output result is the angular distribution of the backscattering coefficient for different values of snow pack parameters.

The incidence angle ranges from 0 • up to 70 • .

[Table 1 about here.]

Table 1 shows the input parameters for the snow backscattering simulations. The selected values are set either according to [START_REF] Rees | Remote Sensing of Snow and Ice[END_REF], [START_REF] Jiang | Ice and water permittivities for millimeter and sub-millimeter remote sensing applications[END_REF] (for the ice and water dielectric permittivities), or according to the local in situ measurements 12 . The surface Root Mean Square (RMS) height represents a potential exception: a 12 mm ground RMS height is not actually corresponding to the terrain roughness. The value of this parameter is strictly imposed by the IEM-B model limitations and by the constrain of using same surface properties for each frequency band (limitting value at 10 GHz) [START_REF] Fung | Microwave Scattering and Emission Models for Users[END_REF].

We summarize in the sequel our conclusions with respect to the dry/wet snow backscattering behavior. Despite the fact that the detection strategy introduced in Section III is designed with X band SAR images, the necessity to demonstrate the particularities of the snow target backscattering in this band, and the potential of broader frequency coverage, imposed the comparative study with respect to C and L bands.

A. Dry snow

According to previous simulation results [START_REF] Besic | Stochastically based wet snow mapping with sar data[END_REF], the dominant component of dry snow backscattering is the underlying ground component for all three frequency bands. This is particularly true at L band, because both snow pack surface and volume component are increasing with the increase in frequency. In general, dry snow backscattering can be considered as not significantly sensitive to the depth of the layer, but fairly sensitive to the snow density.

[Figure 2 about here.] Fig. 2 illustrates the rise in backscattering (underlying ground component) at L band with the increase of snow density.

Firstly noticed by Shi and Dozier, this phenomenon was explained by the fact that the density dominantly affects backscattering through the change of the air-snow interface refraction angle and the change of the wavelength in the snow medium [START_REF] Shi | Estimation of snow water equivalence using sir-c/x-sar, part i: Inferring snow density and subsurface properties[END_REF]. The near correspondence between [START_REF] Shi | Estimation of snow water equivalence using sir-c/x-sar, part i: Inferring snow density and subsurface properties[END_REF] and our qualitative conclusions (with slightly different input parameters) validates the simulation framework proposed in this paper.

At C band, and by assuming the same roughness parameters, the former effect is not present any more: density increase causes decrease of backscattering (Fig. 2). This is because density dominantly affects backscattering through the change in ground-snow dielectric contrast.

The underlying ground component remains dominant at X band. Fig. 2 shows the same backscattering behavior with respect to the density change as for the C band. However, the sensitivity is significantly increased. On one side, the impact of the

1 c Électricité de France (EDF). 2 c Météo-France.
change in the ground-snow dielectric contrast is independent of frequency. On the other side, the ground incident wavelength and the change in ground incident angle are inducing the effect opposite to the one characteristic for the L band.

B. Wet snow

The dominant backscattering component of wet snow is the snow pack surface component [START_REF] Besic | Stochastically based wet snow mapping with sar data[END_REF]. This causes high sensitivity of the wet snow backscattering on wetness change. Therefore, the increased wetness, reflected by an increase in wet snow dielectric constant, leads to stronger total backscattering.

Nevertheless, the initial assumption which justifies the change detection approach is not compromised. The difference in backscattering intensity with respect to the corresponding dry snow is maintained, due to the difference in both the dielectric constant and the roughness properties, between the upper wet snow layer and the underlying ground. This remains true for a wide range of wetness values [START_REF] Hallikainen | Microwave Remote Sensing of Sea Ice[END_REF].

[Figure 3 about here.]

On one side, higher operating frequency would suggest stronger snow-pack surface backscattering. On the other side, the increase in frequency is causing the decrease of the wet snow dielectric constant [START_REF] Rees | Remote Sensing of Snow and Ice[END_REF]. According to our simulator, the backscattering decrease due to the later effect exceeds the former increasing phenomena [START_REF] Besic | Wet snow backscattering sensitivity on density change for swe estimation[END_REF]. This indicates the reduced share of snow pack surface component in total backscattering at higher frequencies.

[Figure 4 2).

[Table 2 about here.]

The backscattering coefficients for different local incidence angles are derived from the corresponding ground range profiles.

Fig. 3 illustrates the comparison between the measured and the scaled simulated backscattering sensitivity on wetness of the upper layer. One can notice that the backscattering intensity is directly proportional to the wetness of the upper layer.

This result, although not qualitatively consistent with the analysis in [START_REF] Arslan | Investigating of snow wetness parameter using a two-phase backscattering model[END_REF], confirms the wet snow pack surface backscattering component is indeed dominant.

In summary, the conclusions derived in this section point out unequivocally that the entire complexity of the snow backscattering should not be neglected in distinguishing between wet and dry snow covers. Moreover, simply assuming that less intense wet snow backscattering is only due to the electromagnetic absorption induced by water inclusions, can be misleading. This is especially true at X band: the dominance of both the underlying ground component of dry snow and the snow pack surface component of wet snow backscattering, indicates that the expected backscattering difference should be treated as a function of the parameters determining surface backscattering. Therefore, it is necessary to introduce a variable threshold (rather than unique) for change detection.

III. THE DETECTION ALGORITHM

In the proposed algorithm, a variable change detection threshold is introduced as a local function of the Local Incidence Angle (LIA) and a global function of the operating frequency. This step is justified by the significant sensitivity of the ratio between the wet and dry backscattering coefficients to the value of local incidence angles (Fig. 6) and by the fact that the operating frequency affects the change detection approach through the threshold value, exclusively.

Aside from those, the speckle noise statistics is used in the estimation procedure, making this method stochastic. It provides the wet snow probability map, which is transformed into the binary map based on the reliability of the variable threshold.

The wet snow detection algorithm (Fig. 4) is consisting of four principal parts.

A. Input data

The following data are required as input:

1) Winter SAR image (W), slant range geometry -image acquired during the winter season, when the dry snow assumption is fairly valid due to the air temperature at ground level;

2) Melting season SAR image (M), slant range geometry -image acquired at the end of the winter season, when the increase in the air temperature causes melting to occur;

3) Digital Elevation Model (DEM) and Sensor & Orbit parameters; 4) Terrain characteristics -approximate information about the terrain roughness and its dielectric properties, 5) Optical snow free terrain image -image acquired during the period when the type of the underlying ground can be observed.

According to our backscattering simulator results (averaged over different underlying ground dielectric constant and grain size values), the assumption that bare ground and dry snow backscattering difference is always smaller than -3 dB is not valid for the X band (Fig. 5). This implies that, at higher frequencies, the difference between bare ground and melting season backscattering coefficients, which is greater than 3 dB, does not unambiguously impose the wet snow presence. The results from Fig. 5 are questioning this assumption for C band. Therefore, the comparison is made directly between the dry snow cover image and the mixed dry/wet snow cover one in order to assure applicability in a wider range of frequencies.

[Figure 5 about here.]

B. SAR image processing

In this part, the input images are calibrated and corregistrated using the resampled SAR geometry intensity simulation. The slant range LIA is also derived.

[Figure 6 about here.]

C. Probability map

The backscattering simulator provides the variable threshold matrix T (Fig. 6), function of the local incidence angle for a given terrain a priori characteristics. These values are derived as the ratio of the simulated wet and dry snow backscattering results, averaged over different wet snow volumetric liquid water content, particle size, and underlying ground dielectric permitivity values (Table 1).

Assuming the gaussianity of the SAR clutter, the intensity over homogeneous regions can be modeled by Gamma probability density function, according to the fully developed speckle model [START_REF] Goodman | Some fundamental properties of speckle[END_REF]:

G(τ |ν, µ) = 1 Γ(ν) ν µ ν τ ν-1 e -νx µ , (2)
with µ being the texture intensity mean, ν -shape factor providing deviation with respect to the corresponding Gaussian distribution and Γ -the Gamma function. If the intensity is expressed as µ ν X, the random variable X follows chi-squared distribution χ 2 (ν) with ν degrees of freedom. The ratio of two chi-squared random variables, normalized with respect to the degrees of freedom, follows Fisher-Snedecor distribution F(ν 1 , ν 2 ) [START_REF] Bombrun | Hierarchical segmentation of polarimetric SAR images using heterogeneous clutter models[END_REF]. This implies that the ratio of two Gamma random variables, having different shape factors but the same mean value is modeled using the Fisher-Snedecor distribution:

F(r|k, ν 1 , ν 2 ) = Γ(ν 1 + ν 2 ) Γ(ν 1 )Γ(ν 2 ) ν 1 kν 2 ( ν1r kν2 ) ν1-1 (1 + ν1r kν2 ) ν1+ν2 . (3)
The proposed stochastic approach is exactly based on the probability estimation relaying on the previously elaborated assumption. The algorithm uses the boxcar neighborhood, coupled with approximate maximum likelihood estimator, in order to obtain local statistics for each of the areas in the image. Due to the poor performances of the Fisher-Snedecor MLE, the estimation is performed rather on the normalized ratio intensity (ξ) than the ratio intensity itself (r = µξ).

If the ratio intensity (r) follows the Fisher-Snedecor distribution, the normalized ratio intensity (ξ) is modeled by Beta prime distribution:

B (ξ|ν 1 , ν 2 ) = Γ(ν 1 + ν 2 ) Γ(ν 1 )Γ(ν 2 ) ξ ν1-1 (1 + ξ) ν1+ν2 . (4)
For each local neighborhood, we are deriving the mean value (µ), normalizing the texture locally and estimating ν 1 and ν 2 parameters. This way, we are able to define the probability density function (B (i,j) ) for every possible region in the image. By integrating the obtained PDF with respect to the normalized texture, we are getting the cumulative distribution function (B (i,j) )

for the normalized threshold of the central pixel:

i out (i, j) = B (i,j) T (i, j) µ(i, j) = T (i,j) µ(i,j) -∞ B (i,j) (ξ)dξ. ( 5
)
This value is the probability that the ratio is smaller or equal to the local threshold or exactly the wet snow probability.

[Figure 7 about here.]

D. Confidence level

The obtained wet snow probability map is transformed into the wet snow binary map by defining the confidence level C:

I out (i, j) = 1 if i out (i, j) ≥ C, 0 if i out (i, j) < C. ( 6
)
The confidence level value is related to the underlying terrain properties. Namely, the simulator generating the variable threshold is assuming the bare sol as the underlying ground. Any kind of model mismatch is decreasing the reliability of the computed variable threshold. In the final step, this effect can be compensated through the confidence level value: globally or locally, depending on the degree of the a priori know terrain characteristics.

Having the optical image of the targeted snow-free terrain among the input data, we are able to classify different distributed targets. We are locally assigning the confidence level which deviates from the full confidence (one) proportionally to the distance between the given class and the bare ground in the orthogonalized space of the given image (Fig. 8). This way we are accounting for the presence of low vegetation, forests, lakes and glaciers in the region of interest, compromising the validity of the derived threshold.

Otherwise, the global confidence level is estimated based on the overall validity of the assumed terrain characteristics, the ground dielectric constant particularly.

[Figure 8 about here.]

IV. RESULTS AND DISCUSSION

The proposed algorithm is illustrated by the results obtained using two X-band TerraSAR-X stripmap images (Fig. 7) acquired over the Grand Rousse massif near Grenoble, France:

• the winter image acquired on the 08th of February 2009 holding for the dry snow assumption according to the local meteorological data,

• the melting season image acquired on the 02nd of March 2009 in the presence of wet snow, according to the same source.

The local incidence angle map is computed using the Digital Terrain Model (Datum: WGS-84, UL Geo: For comparison, we provide in Fig. 8 several pairs of binary maps (HH and VV) derived using different change detection methods. The first set of images (a,e) present the results obtained by applying the variable threshold on the ratio of original input images. The following set (b,f) is derived by performing spatial speckle filtering on the input images before calculating the ratio of SAR intensities. The last two sets are the final outputs of our algorithm (before georeferencing): firstly by thresholding the probability maps with a global value of confidence level, and secondly by using a predefined confidence level map for locally thresholding the probability. This map has been derived from an optical image (SPOT) acquired during the snow-free season.

The quantitative correspondence between the HH and VV maps, as it can be seen in Fig. 8, augments in case of the proposed stochastic approach.

The qualitative verification of the proposed method validity is performed using the local temperature measurements at the ground level. The measurements data, provided by the Électricité de France (EDF), are related to the area of interest (Grand Rousse massif) at the same date (03 March 2009), meaning that comparison with the available SAR images was possible.

The procedure is based on the stated fact that wet snow presence is characteristic for the local ground temperature above 0 • C, while dry snow can be found below 0 • C. The temperature measurements are compared to the obtained wet snow probability maps (Fig. 9). The approximate matching between the high temperature regions and the high probability regions, two information sources which could be considered as independent, is pointing to the validity of the proposed method.

[Figure 9 about here.]

Using the matching percentage between the HH and the VV binary maps as a quantitative indicator (Fig. 8), the proposed method can be considered as an improvement in the wet snow detection. This is particularly true when using the locally predefined confidence level. The implicitly introduced spatial correlation between the wet snow areas, on one side, and the possibility to subtly account for the presence of different distributed targets underlying the snow layer by thresholding the obtained probability map, on the other side, bring the introduced stochastic approach superior with respect to conventional change detection approaches.

Note that the illustrated relation between the bare ground and dry snow backscattering (Fig. 5) indicates that the eventual presence of the bare ground in the melting season image cannot cause a false wet snow detection.

V. CONCLUSIONS

Based on fundamental scattering theories, the proposed backscattering simulator allowed us to demonstrate the different backscattering mechanisms of the dry and the wet snow cover. The detailed sensitivity analysis pointed out the dominance of the underlying ground backscattering component in case of dry snow and the dominance of snow pack surface component in case of wet snow. The later conclusion is confirmed by the performed in situ measurements. Rather than assuming the cause of the backscattering difference to be the absorption phenomena of the wet snow due to the water inclusions, we provided a detailed comparative description of the backscattering. It resulted in the conclusions that the difference occurs to be significantly dependent on the local incidence angle and the operating frequency. Consequently, we have proposed an improvement of the conventionally used wet snow detection method, by introducing a variable threshold given as function of the local incidence angle for each pixel and function of the operating frequency, also. This fact, as well as the choice of the winter image as the reference, allows application in wider frequency range. Instead of directly thresholding the ratio of multitemporal images, by considering the spatial correlation, we rather estimate the probability of the wet snow occurrence, making the algorithm stochastic.

Further refinements of the proposed detection method are going in two main directions. In the first, the other parameters significantly affecting surface backscattering (e.g. roughness and dielectric constant) would be also included as spatial variables in the threshold derivation. These parameters will be either retrieved by means of spaceborne remote sensing or by interpolating the local in situ measurements. Concerning the second direction, further work will mostly consist in adjusting the proposed method to polarimetric SAR input data [START_REF] Lessard-Fontaine | Multi-temporal dry and wet snow mapping in alpine context using polarimetric radarsat-2 time-series[END_REF]. The idea is to exploit the dual-pol images in the probability derivation, rather than just using them in the verification of the results.

Input parameter

  about here.] 1) "In situ" measurements: This paragraph deals with the analysis of X-band ground based radar signals (VV polarization channel) acquired on the 23rd of April 2013, in the Grand Rousse region of the French Alpes. Two sets of measurements were performed over the same spot (same surface parameters) at 4 hour interval. The first set of range profiles was acquired in the morning, while the second acquisition took place in the afternoon, allowing snow wetness variation in the upper snow layer due to the sunny weather conditions. The parameters of the upper snow layer which characterize the snow pack surface backscattering component were carefully measured (see Table

Fig. 2 :

 2 Fig. 2: Snow backscattering sensitivity on density change (VV polarisation (a) dry snow, L band; (b) wet snow, L band; (c) dry snow, C band; (d) wet snow, C band; (e) dry snow, X band; (f) wet snow, X band. figure

Fig. 3 :FigFig. 5 :

 35 Fig. 3: Ground based X-band radar in situ measurements compared to simulation results. figure

Fig. 6 :Fig. 7 :Fig. 8 :

 678 Fig. 6: Variable thresholds function of the LIA. figure

TABLE I :

 I Input parameters for the snow backscattering simulations: values in brackets are used for computing the variable threshold (median values for particle equivalent radius and ground dielectric constant are employed in Section II for the snow backscattering analysis).

					Backscattering coefficient for different ice volume fractions									Backscattering coefficient for different volumetric liquid water content (density)
		-5													-10
												180kg/m3										3.6% (220kg/m3)
		-6										190kg/m3 200kg/m3			-12							5.6% (240kg/m3) 7.6% (260kg/m3)
	Backscattering coefficient [dB]	-12 -11 -10 -9 -8 -7										210kg/m3 220kg/m3	Backscattering coefficient [dB]	-24 -22 -20 -18 -16 -14							9.6% (280kg/m3) 11.6% (300kg/m3)
		-13													-26
		10 -14	15	20	25	30	35 LIA [degree]	40	45	50	55	Dry snow 60 10 -28	15	20	25	30	35 Wet snow LIA [degree]	40	45	50	55	60
				Snow density		(a)			180,190,200,210,220 kg/m 3	220,240,260,280,300 kg/m 3 (b)
				Volumetric liquid water content			0 %					3.6,5.6,7.6,9.6,11.6 (7,9,11,13) %
				Snow depth												2 m
				Frequency Particle equivalent radius Parameter Ice dielectric constant Acquisition time Water dielectric constant Snow dielectric constant Snow RMS height Volumetric liquid water content Snow correlation length Snow density Ground dielectric constant Snow depth	1.5 (L), 6 (C), 10 (X) GHz Sample I Sample II 225 (121,169,225,289,361) µm 11h00 15h00 3.15 2.94 4.39 80 + j50 06.8 % 16.00 % 10 mm 175 kg/m 3 250 kg/m 3 50 mm 6 (5,5.5,6,6.5,7) 2.15 m 2.15 m
	table		Ground RMS height Ground correlation length TABLE II: In situ measurements: snow parameters. 12 mm 50 mm
	table						Fig. 1: Snow pack backscattering mechanism
	figure

ACKNOWLEDGMENT

The authors would like to thank to the Électricité de France (EDF) and the Météo-France for providing us the ground truth temperature data and other input parameters of interest, and to the Institut national de l'information géographique et forestière (IGN France) for the DEM. The SPOT optical image was provided by the ISIS 2011-593 contract, copyright CNES/Spot-Image.


Sébastien Marmin

David Ginsbourger

Jean Baccou

Jacques Liandrat

Frédéric Perales

Adaptive design experiment and non-stationary kriging : applications to uncertainty analysis in nuclear mechanical studies

Introduction

Three propagation behaviour

Around the inclusion

Through the inclusion

Intermediate

It is related to the study of the heterogeneous material cracking (typically concrete) and to the simulation of the cracking energy by the code XPER. Two uncertained input parameters are considered : the length of the inclusions, L and the ratio of interface energies, W. An adaptive design of experiments is here relevant since we observe a high variation of the response E corresponding to the transition between two cracking propagation modes.

This work aims at integrating uncertainties in mechanical analysis with complex computer codes from nuclear safety studies. More precisely, it is intended to analyse if the variations of the uncertain code inputs can drastically change the behaviour of the system or even turn the system into a non-secure state. To do so, numerical design of experiments methods based on spatial statistics can be used for launching new simulations in the most relevant areas (high variation areas, over sets, ...). These methods are based on a modeling of the response with an interpolating Gaussian process. The covariance of the process is often assumed stationary, whereas this strong hypothesis leads sometime to poor predictions.

Context A f i rst example of mechanical application

Objectives build new non-stationary Gaussian process models to describe the heterogeneous behaviour of the response.

integrate the new models in a adaptive design of experiment. validate the methods on mechanical applications related to the cracking of heterogeneous materials.

  process Y(x) of mean µ(x) and of covariance C(x,x')

+ Data conditioning

Averaging

Predictions m(x) with a variance s(x)

Gaussian process modeling

Non-stationary model

Results

Predictive error

We have tested seven Gaussian process models. We show that the models with warped Gaussian processes adapted to the data structure is more effective.

Criteria comparison

Stationary model

Perspectives

We will compare new derivative-based criteria with further norm and second order derivatives. Then we will aim at improving the non-stationary Gaussian process modeling by matching them with wavelets methods from the image processing f i eld. Finally, new mechanical applications involving more uncertain code inputs and different types of non stationarity will be considered We can interpret the idea of warping a Gaussian process to make it non-stationary as a transformation of the input space. The nonstationary covariance is the composition of a stationary covariance with a transformation function T locally determined by a density function g.

Sample of a stationary process

-1,0 0,0 1,0

Sample of a nonstationary Gaussian process > 1 contraction < 1 dilatation

Warping Inf i ll criteria

Predictive error

To compare inf i ll criteria, we evaluate the prediction errors in the high variation region after adding 10 points. We repeat the experiment 10 times with differents initial designs.

On the considered test case, we observed good performance of G-IMSE for the prediction in high-variation zone.

We use two variance-based criteria:

We def i ne two gradient-based criteria focusing more on variation region: