Emmanuel Chevallier
email: emmanuelchevallier1@gmail.com

Frédéric Barbaresco

Jesus Angulo

Jesús Angulo

Kernel density estimation on the

Keywords: Kernel density estimation, Siegel space, Symmetric space, Radar signals

come

Introduction

Probability density estimation is a vast topic. There exists various standard techniques in the Euclidean context, such as histograms, kernel methods, or the characteristic function approach. These methods can sometimes be transposed to the case of Riemannian manifolds. However, the transposition often introduces additional computational efforts. This additional effort depends on the method used and the nature of the manifold. The Siegel space is a generalization of the hyperbolic space. It has a structure of symmetric Riemannian manifold, which enables to adapt different density estimation methods at a reasonable cost. Convergence rates of the density estimation using kernels and orthogonal series were progressively generalized to Riemannian manifolds, see [START_REF] Hendriks | Nonparametric estimation of a probability density on a Riemannian manifold using Fourier expansions[END_REF] [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF]. The Siegel space appears in radar processing in the study of Toeplitz block Toeplitz matrices, whose blocks represent covariance matrices of the signal, see [START_REF] Barbaresco | Robust statistical radar processing in Frchet metric space: OS-HDR-CFAR and OS-STAP processing in siegel homogeneous bounded domains[END_REF][START_REF] Barbaresco | Information Geometry of Covariance Matrix: Cartan-Siegel Homogeneous Bounded Domains, Mostow/Berger Fibration and Frchet Median[END_REF][START_REF] Barbaresco | Information geometry manifold of Toeplitz Hermitian positive definite covariance matrices: Mostow/Berger fibration and Berezin quantization of Cartan-Siegel domains[END_REF]. Information geometry is now a standard framework in radar processing, see [START_REF] Barbaresco | Robust Median-Based STAP in Inhomogeneous Secondary Data : Frechet Information Geometry of Covariance Matrices[END_REF][START_REF] Barbaresco | Robust statistical radar processing in Frchet metric space: OS-HDR-CFAR and OS-STAP processing in siegel homogeneous bounded domains[END_REF][START_REF] Barbaresco | Information Geometry of Covariance Matrix: Cartan-Siegel Homogeneous Bounded Domains, Mostow/Berger Fibration and Frchet Median[END_REF][START_REF] Barbaresco | Information geometry manifold of Toeplitz Hermitian positive definite covariance matrices: Mostow/Berger fibration and Berezin quantization of Cartan-Siegel domains[END_REF][START_REF] Degurse | A Riemannian Approach for Training Data Selection in Space-Time Adaptive Processing Applications[END_REF][START_REF] Degurse | Information Geometry for radar detection in heterogeneous environments[END_REF][START_REF] Barbaresco | Koszul Information Geometry and Souriau Geometric Temperature/Capacity of Lie Group Thermodynamics[END_REF]. The information geometry on positive definite Teoplitz block Teoplitz matrices is directly related to the metric on the Siegel space, see [START_REF] Jeuris | The Kahler mean of Block-Toeplitz matrices with Toeplitz structured block[END_REF]. Indeed these Toeplitz block Toeplitz matrices can be represented by a symmetric positive definite matrix and a set of coefficients laying in the Siegel disk. The metric considered on Toeplitz block Toeplitz matrices is a product metric between a metric on symmetric positive definite matrices and the Siegel disk metric, see [START_REF] Jeuris | The Kahler mean of Block-Toeplitz matrices with Toeplitz structured block[END_REF]. One already encounters the problem of density estimation in the hyperbolic space for electrical impedance [START_REF] Huckemann | Mobius deconvolution on the hyperbolic plan with application to impedance density estimation[END_REF], networks [START_REF] Asta | Geometric network comparison[END_REF] and radar signals [START_REF] Chevallier | Probability density estimation on the hyperbolic space applied to radar processing[END_REF]. In [START_REF] Said | New Riemannian Priors on the Univariate Normal Model[END_REF] was proposed a generalization of the Gaussian law on the hyperbolic space. Apart from [13], where authors propose a generalization of the Gaussian law, probability density estimation on the Siegel space has not yet been addressed. The contributions of the paper are the following. We review the main non parametric density estimation techniques on the Siegel disk. We provide some rather simple explicit expressions of the kernels proposed by Pelletier in [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF]. These expressions makes the kernel density estimation the most adapted method. We present visual results of estimated densities in the simple case where the Siegel disk correspond to the Poincaré disk. The paper begins with an introduction to the Siegel space in Section 2. Section 3 reviews the main non-parametric density estimation techniques on the Siegel space. Section 4 presents an application to radar data estimation.

The Siegel space

This section presents facts about the Siegel space. The interested reader can find more details in [START_REF] Siegel | Symplectic geometry[END_REF][START_REF] Terras | Harmonic analysis on symmetric spaces and applications II[END_REF]

The Siegel upper half space

The Siegel upper half space is a generalization of the Poincaré upper half space. Let Sym(n) be the space of real symmetric matrices of size n × n and Sym + (n) the set of real symmetric positive definite matrices of size n × n. The Siegel upper half space is defined by

H n = {Z = X + iY |X ∈ Sym(n), Y ∈ Sym + (n)}
H n is equipped with the following metric

ds = 2.tr(Y -1 dZY -1 dZ)
The set of real symplectic matrices Sp(n, R) is defined by

g ∈ Sp(n, R) ⇔ g t Jg = J where J = 0 I n -I n 0 and I n is the n × n identity matrix. Sp(n, R) is a subgroup of SL 2n (R), the set of 2n × 2n invertible matrices of determinant 1. Let g = A B C D ∈ Sp(n, R).
The metric ds is invariant under the following action of Sp(n, R),

g.Z = (AZ + B)(CZ + D) -1 .
This action is transitive, i.e.

∀Z ∈ H n , ∃g ∈ Sp(n, R), g.iI = Z.
The stabilizer K of iI is the set of elements g of Sp(n, R) whose action leaves iI fixed. K is a subgroup of Sp(n, R) called the isotropy group. We can verify that

K = A B -B A , A + iB ∈ SU (n) .
A symmetric space is a Riemannian manifold where the reversal of the geodesics is well defined and is an isometry. Formally, exp p (u) → exp p (-u) is an isometry for each p on the manifold, where u is a vector in the tangent space at p, and exp p the Riemannian exponential application at p. In other words, the symmetry around each point is an isometry. H n is a symmetric space, see [START_REF] Terras | Harmonic analysis on symmetric spaces and applications II[END_REF]. The structure of a symmetric space can be studied through its isometry group and the Lie algebra of its isometry group. The present work will make use of the Cartan and Iwasawa decompositions of the Lie algebra of Sp(n, R). Let sp(n, R) be the Lie algebra of Sp(n, R). Given A,

= ∪ k∈K Ad k (a) ( 2
)
where Ad is the adjoint representation of Sp(n, R).

The Siegel disk

The Siegel disk D n is the set of complex matrices {Z|I -Z * Z ≥ 0} where ≥ stands for the Loewner order. Recall that for A and B to Hermitian matrices, A ≥ B according to the Loewner order means that A -B is positive definite. The transformation where A stands for the conjugate of A. The point iI in H n is identified with the null matrix noted 0 in D n . Let Z ∈ D n . There exists P a diagonal matrix with decreasing positive real entries and U a unitary matrix such that Z = U P U t . Let τ 1 ≥ ... ≥ τ n be such that

Z ∈ H n → (Z -iI)(Z + iI) -1 ∈ D n
P =    th(τ 1 ) . . . th(τ n )    . Let A 0 =    ch(τ 1 ) . . . ch(τ n )    , B 0 =    sh(τ 1 ) . . . sh(τ n )    and g Z = U 0 0 U . A 0 B 0 A 0 B 0 .
It can be shown that

g Z ∈ Sp(n, C) and g Z .0 = Z (3)
We provide now a correspondence between the elements of D n and the elements of p defined in Eq. 1. Let

H Z =           τ 1 . . . τ n -τ 1 . . . -τ n           ∈ a (4)
and

a Z =           e τ 1 . . .
e τn e -τ 1 . . .

e -τn           ∈ A = exp(a).
It can be shown that there exists k ∈ K such that

Cexp(Ad k (H Z ))C -1 .0 = Z or equivalently Cka Z kC -1 .0 = Z.
Recall that Eq. 2 gives Ad k (H) ∈ p and kak ∈ exp(p). The distance between Z and 0 in D n is given by

d(0, Z) = 2 τ 2 i 1/2 , (5)
see [START_REF] Terras | Harmonic analysis on symmetric spaces and applications II[END_REF] page 292. !!Negative curvature of the Siegle space!!

Non parametric density estimation on the Siegel space

Let Ω be a space, endowed with a σ-algebra and a probability measure p.

Let X be a random variable Ω → D n . The Riemannian measure of D n is called vol and the measure on D n induced by X is noted µ X . We assume that µ X has a density, noted f , with respect to vol, and that the support of X is a compact set noted Supp. Let (x 1 , .., x k ) ∈ D k n be a set of draws of X. The Dirac measure in point a is defined as:

δ a (U ) : 1 if a ∈ U 0 if a ∈ U Let µ k = 1 k k i=1 δ x i
denotes the empirical measure of the set of draws. This section presents four non-parametric techniques of estimation of f from the set of draws (x 1 , .., x k ). The estimated density at x in D n is noted fk (x) = f (x, x 1 , ..., x k ). The relevance of a density estimation technique depends on several aspects. When the space allows it, the estimation technique should equally consider each directions and locations. This leads to an isotropy and a homogeneity condition. In the kernel method for instance, a kernel function K x i is placed at each observation x i . Firstly, in order to treat directions equally, the function K x i should be invariant under the isotropy group of x i . Secondly, for another observation x j , functions K x i and K x j should be similar up to the isometries that send x i on x j . These considerations strongly depend on the geometry of the space: if the space is not homogeneous and the isotropy group is empty, these indifference principles have no meaning. Since the Siegel space is symmetric, it is homogeneous and has a non empty isotropy group. Thus the density estimation technique should be chosen accordingly. The convergence of the different estimation techniques is widely studied. Results were first obtained in the Euclidean case, and are gradually extended to the probability densities on manifold, see [START_REF] Hendriks | Nonparametric estimation of a probability density on a Riemannian manifold using Fourier expansions[END_REF][START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF][START_REF] Huckemann | Mobius deconvolution on the hyperbolic plan with application to impedance density estimation[END_REF]. The last relevant aspect, is computational. Each estimation technique has its own computational framework, which presents pro and cons given the different applications. For instance, the estimation by orthogonal series presents an initial pre-processing, but provides a fast evaluation of the estimated density in compact manifolds.

Histograms

The histogram is the simplest density estimation method. Given a partition of the space D n = ∪ i A i , the estimated density is given by

f (x ∈ A i ) = 1 vol(A i ) k j=1 1 A i (x j )
where 1 A i stands for the indicator function of A i . Following the considerations of the previous sections, the elements of the partition should firstly be as isotropic as possible, and secondly as similar as possible to each other. Regarding the problem of histograms, the case of the Siegel space is similar to the case of the hyperbolic space. There exist various uniform polygonal tilings on the Siegel space that could be used to compute histograms. However, there are ratio λ ∈ R for which there is no homothety. Thus it is not always possible to scale the size of the bins to a given set of draws of the random variable. Modifying the scale of the density estimation requires then a change of the structure of the tiling. Thus the study of histograms has not been deepened.

Orthogonal series

The estimation of the density f can be made out of the estimation of the scalar product between f and a set of orthonormal functions {e j }. The most standard choice for {e j } is the eigenfunctions of the Laplacian. When the variable X takes its values in R n , this estimation technique becomes the characteristic function method. When the underlying space is compact, the spectrum of the Laplacian operator is countable, while when the space is non-compact, the spectrum is uncountable. In the first case, the estimation of the density f is made through the estimation of a sum, while in the second case is made through the estimation an integral. In practice, the second situation present a larger computational complexity. Unfortunately, eigenfunctions of the Laplacian operator are known on D n but not on compact sub-domains. For this reason the study of this method has not been deepened.

Kernels

Let K : R + → R + be a map which verifies the following properties: i)

R d K(||x||)dx = 1, ii) R d xK(||x||)dx = 0, iii) K(x > 1) = 0, iv) sup(K(x)) = K(0).
Let p ∈ D n . Generally, given a point p on a Riemannian manifold, exp p defines an injective application only on a neighborhood of 0. The Siegel space is a non-compact symmetric space and has thus only negative sectional curvatures. Thus, exp p is injective on the whole space. When the tangent space T p D n is endowed with the local scalar product,

θ p : q → θ p (q) = dvol dexp * p (Leb p ) (q), (6)
is the density of the Riemannian measure of D n with respect to the Lebesgue measure Leb p after identification of D n and T p D n induced by exp p , see Fig. 1.

Figure 1: M is a Riemannian manifold, T x M is its tangent space at x. The application induces a volume change θ x between T x M and M.

Given K and a positive radius r, the estimator of f proposed by [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF] is defined by: fk

= 1 k i 1 r n 1 θ x i (x) K d(x, x i ) r . (7)
The corrective factor θ x i (x) -1 is necessary since the kernel K originally integrates to one with respect to the Lebesgue measure, instead of the Riemannian measure. It can be noted that this estimator is the usual kernel estimator in the case of Euclidean space. When the curvature of the space is negative, which is the case of the Siegel space, the distribution placed over each sample x i has x i as intrinsic mean. The following theorem provides convergence rate of the estimator. It is a direct adaptation of theorem 3.1 of [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF]. Theorem 3.1. Let (M, g) be a Riemannian manifold of dimension n and µ its Riemannian volume measure. Let X be a random variable taking its values in a compact subset C of (M, g). Let 0 < r ≤ r inj , where r inj is the infimum of the injectivity radius on C. Assume the law of X has a twice differentiable density f with respect to the Riemannian volume measure. Let fk be the estimator defined in Eq. 7. The exist a constant C f such that

x∈M E x 1 ,...,x k [(f (x) -fk (x)) 2 ]dµ ≤ C f ( 1 kr n + r 4 ) (8) If r ∼ k -1 n+4 , x∈M E x 1 ,...,x k [(f (x) -fk (x)) 2 ]dµ = O(k -4 n+4 ) (9)
Proof. See appendix A.

Since the Siegel space has negative sectional curvatures, r inj = +∞. It can be easily verified that for an isometry α we have: fk (x, x 1 , ..., x k ) = fk (α(x), α(x 1 ), ..., α(x k )).

Each location and direction are processed as similarly as possible.

In order to obtain the explicit expression of the estimator one must first have the explicit expression of the Riemannian exponential, its inverse, and of the function θ p , see Eq. 6-7. These expressions are difficult and sometimes impossible to obtain on general Riemannian manifolds. In the case of the Siegel space, the symmetric structure makes the computation possible. Since the space is homogeneous, the computation can be made at the origin iI ∈ H n or 0 ∈ D n and transported to the whole space. In the present work, the random variable lays in D n . However in the literature the Cartan and Iwasawa decompositions are usually given for the isometry group of H n . Thus our computation starts in H n before moving to D n .

The Killing form on the Lie algebra sp(n, R) of the isometry group of H n induces a scalar product on p. This scalar product can be transported on exp(p) by left multiplication. This operation gives exp(p) a Riemannian structure. It can be shown that on this Riemannian manifold, the Riemannian exponential at the identity coincides with the group exponential. Furthermore,

φ : exp(p) → H n g → g.iI (10)
is a bijective isometry, up to a scaling factor. Since the volume change θ p is invariant under rescaling of the metric, this scaling factor has no impact. Thus H n can be identified with exp(p) and exp iI∈Hn with exp |p . The expression of the Riemannian exponential is difficult to obtain in general, however it boils down to the group exponential in the case of symmetric spaces. This is the main element of the computation of θ p . The Riemannian volume measure on exp(p) is noted vol . Let

ψ : K × a → p (k, H) → Ad k (H)
Let a + be the diagonal matrices with strictly decreasing positive eigenvalues. Let Λ + be the set of positive roots of sp(n, R) as described in [START_REF] Terras | Harmonic analysis on symmetric spaces and applications II[END_REF] page 282, Λ + = {e i + e j , i ≤ j} ∪ {e i -e j , i < j} where e i (H) is the i-th diagonal term of the diagonal matrix H. Let C c (E) be the set of continuous compactly supported function on the space E. In [START_REF] Gangolli | Harmonic analysis of spherical functions on real reductive groups[END_REF] is given page 73 that for all t ∈ C c (p), there exists

c 1 > 0 such that p t(Y )dY = c 1 K a + t(ψ(k, H)) λ∈Λ + λ(H)dkdH, ( 11
)
where dY is a Lebesgue measure on the coefficients of Y . Let p = ψ(K ×a + ). λ ∈ Λ + never vanishes on a + and p \ p has a null measure. Thus

p t(Y ) λ∈Λ + 1 λ(H Y ) dY = c 1 K a + t(Ad k (H))dkdH, (12)
where H Y is the point in a + such that there exists k in K such that ψ(k, H Y ) = Y . Calculation in [START_REF] Gangolli | Harmonic analysis of spherical functions on real reductive groups[END_REF] page 73 also gives that for all t ∈ C c (p), there exists c 2 > 0 such that

Sp(n,R) t(g)dg = c 2 K a + K t(k 2 .exp(Ad k 1 (H)))J(H)dk 1 dHdk 2 , ( 13
)
where dg is the Haar measure on Sp(n, R) and

J(H) = λ∈Λ + e λ(H) -e -λ(H) = 2 |Λ + | λ∈Λ + sinh(λ(H)).
Thus for all t ∈ C c (Sp(n, R)/K),

Sp(n,R)/K t(x)dx = c 2 K a + t(exp(Ad k (H)))J(H)dkdH, ( 14
)
where dx is the invariant measure on Sp(n, R)/K. After identifying Sp(n, R)/K and exp(p) the Riemannian measure on exp(p) coincides with the invariant measure on Sp(n, R)/K. Thus for all t ∈ C c (exp(p)),

exp(p) t(x)dvol = c 2 K a + t(exp(Ad k (H)))J(H)dkdH. ( 15
)
Using the notation H Y of Eq. 12,

p t(exp(Y ))J(H Y ) λ∈Λ + 1 λ(H Y ) dY = c 1 K a + t(exp(Ad k (H)))J(H)dkdH.
(16) Combining Eq. 15 and Eq. 16 we obtain that exists c 3 such that p t(exp(Y ))

λ∈Λ + sinh(λ(H Y )) λ(H Y ) dY = c 3 exp(p) t(x)dvol . (17)
The term sinh(λ(H))

λ(H)
can be extended by continuity on a thus p t(exp(Y ))

λ∈Λ + sinh(λ(H Y )) λ(H Y ) dY = c 3 exp(p) t(x)dvol . (18)
Let dY be the Lebesgue measure corresponding to the metric. Then the exponential application does not introduce a volume change at 0 ∈ p. Since H 0 = 0 and sinh(λ(H)) λ(H) -→ H→0 1, we have c 3 = 1. Let log denote the inverse of the exponential application. We have

dlog * (vol ) dY = λ∈Λ + sinh(λ(H Y )) λ(H Y ) .
Since φ from Eq. 10 is an isometry up to a scaling factor, if Y ∈ p and

Cφ(exp(Y ))C -1 = exp 0 (u ∈ T 0 D n ), then dlog * (vol) dLeb 0 (u) = dlog * (vol ) dY (Y ),
where Leb 0 refers to the Lebesgue measure on the tangent space T 0 D n as in Eq. 6. Given Z ∈ D n , H Z from Eq. 4 verifies Cφ(exp

(Ad k (H Z )))C -1 = Z for some k in K. Thus θ 0 (Z) = dlog * (vol ) dY (Ad k (H Z )) = λ∈Λ + sinh(λ(H Z )) λ(H Z ) .
We have then

θ 0 (Z) = i<j sinh(τ i -τ j ) τ i -τ j i≤j sinh(τ i + τ j ) τ i + τ j
where the (τ i ) are described in section 2.2. Given

Z 1 , Z 2 ∈ D n , θ Z 1 (Z 2 ) = θ 0 (g -1 Z 1 .Z 2 )
where g -1 Z 1 is defined in Eq. 3. It is thus possible to use the density estimator defined in Eq. 7. Indeed,

1 θ Z 1 (Z 2 ) K(d(Z 1 , Z 2 )) = i<j τ i -τ j sinh(τ i -τ j ) i≤j τ i + τ j sinh(τ i + τ j ) K(( 2
τ 2 i ) 1/2 ) ( 19
)
where the (τ i ) are the diagonal elements of H g -1

Z 1 .Z 2 .
4 Application to radar processing

Radar data

In STAP radar processing, the signal is formed by a succession of matrix X representing the realization of a temporal and spatial process. Let B + n,m be the set of positive definite block Teoplitz matrices composed of n × n blocks of m × m matrices (PD BT). For a stationary signal the autocorrelation matrix R is PD BT, see [START_REF] Jeuris | The Kahler mean of Block-Toeplitz matrices with Toeplitz structured block[END_REF]. Authors of [START_REF] Jeuris | The Kahler mean of Block-Toeplitz matrices with Toeplitz structured block[END_REF] proposed a generalization of Verblunsky coefficients and defined a parametrization of PD BT matrices,

B + n,m → Sym + × D m-1 n R → (P 0 , Ω 1 , .., Ω m-1 ) (20)
in which the metric induced by the Kähler potential is the product metric of an affine invariant metric on Sym + and m -1 times the metric of the Siegel disk, up to a scaling factor. Among other references, positive definite block Teoplitz matrices have been studied in the context of STAP-radar processing in [START_REF] Barbaresco | Robust statistical radar processing in Frchet metric space: OS-HDR-CFAR and OS-STAP processing in siegel homogeneous bounded domains[END_REF][START_REF] Barbaresco | Information Geometry of Covariance Matrix: Cartan-Siegel Homogeneous Bounded Domains, Mostow/Berger Fibration and Frchet Median[END_REF][START_REF] Barbaresco | Information geometry manifold of Toeplitz Hermitian positive definite covariance matrices: Mostow/Berger fibration and Berezin quantization of Cartan-Siegel domains[END_REF].

Experimental results

In this section we show density estimation results of the marginal parameters Ω k . For the sake of visualization, only the Siegel disk D 1 is considered.

Recall that D 1 coincides with the Poincaré disk. The results are partly extracted from the conference paper [START_REF] Chevallier | Probability density estimation on the hyperbolic space applied to radar processing[END_REF]. Data used in the experimental tests are radar observations from THALES X-band Radar, recorded during 2014 field trials campaign at Toulouse Blagnac Airport for European FP7 UFO study (Ultra-Fast wind sensOrs for wake-vortex hazards mitigation). Data are representative of Turbulent atmosphere monitored by radar. Fig. 23illustrate the density estimation of six coefficients on the Poincaré unit disk. Fig. 2 correspond to a clear environment and Fig. 3 to a rainy environment. The densities are individually re-scaled for visualization purposes. For each environment the dataset is composed of 120 draws. The densities of the coefficients Ω k are representative of different backgrounds. These information on the background are expected to ease the detection of interesting targets.

Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 Ω 6
Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 Ω 6

Conclusion

Three non parametric density estimation techniques have been considered. The main advantage of histograms in the Euclidean context is there simplicity of use. This make histograms an interesting tool despite the fact that they do not present optimal convergence rate. On the Siegel space, histograms loses their simplicity advantage. They were thus not deeply studied. The orthogonal series density estimation also present technical disadvantages on the Siegel space. Indeed, the series become integrals which make the numerical computation of the estimator more difficult than in the Euclidean case. On the other hand, the use of the kernel density estimator does not present major differences with the Euclidean case. The convergence rate obtained in [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF] can be extended to compactly supported random variables on non compact Riemannian manifolds. Furthermore the corrective term whose computation is required to use Euclidean kernels on Riemannian manifolds turns out to have a reasonably simple expression. Our future efforts will concentrate on the use of kernel density estimation on the Siegel space in radar signal processing. We strongly believe that the estimation of the densities of the Ω k will provide a interesting description of the different backgrounds.

The authors would like to thank Salem Said, Michal Zidor and Dmitry Gourevitch for the help they provided in the understanding of symmetric spaces and the Siegel space.

A Demonstration of theorem 3.1 Lemma A.1. Let (M, g) be a Riemannian manifold, let C be a compact subset of M and let U be a relatively compact open subset of M containing C. Then there is a compact Riemannian manifold (M , g ) such that U is an open subset of M , the inclusion i : U → M is a diffeomorphism onto its image and g = g on U .

Proof. We can assume that M is not compact. Let f : M → R be smooth function M which tends to +∞ at infinity. Since U is compact, f -1 (] -∞, a[) contains U for a large enough. By Sard Theorem there exists a value a ∈ R such that f -1 (a) contains no critical point of f and such that f -1 (]-∞, a[) contains U . It follows that N = f -1 (]-∞, a]) is a submanifold with boundary of M . Since f tends to +∞ at infinity N is compact as well as its boundary ∂N = f -1 ({a}). Call M the double of N . It is a compact manifold which contains N such that the inclusion i : N → M is a diffeomorphism onto its image (see [START_REF] Munkres | Elementary Differential Topology[END_REF], Theorem 5.9 and Definition 5.10 ). Choose any metric g 0 on M . Consider two open subsets W 1 and W 2 in M and two smooth functions f 1 , f 2 : M → [0, 1] such that

U ⊂ W 1 ⊂ W 1 ⊂ W 2 ⊂ W 2 ⊂ int N the interior of N , f 1 (x) = 1
on W 1 , vanishes outside of W 2 , and

f 2 (x) = 1
outside W 1 and vanishes in U . Define g on M by g = f 1 g + f 2 g 0 on N and g = f 2 g 0 outside of N . Since f 1 + f 2 > 0, g is positive definite everywhere on M . Since f 1 vanishes outside of W 2 , g is smooth on M . Finally, since f 1 = 1 and f 2 = 0 on U , g = g on U .

B

  and C three real n by n matrices, let us write A B C -A t = (A, B, C). We have that sp(n, R) = {(A, B, C)|B and C symmetric} The Cartan decomposition of sp(n, R) is given by sp(n, R) = t ⊕ p where t = {(A, B, -B)|B symmetric and A skew-symmetric} p = {(A, B, B)|A, B, symmetric} (1) The Iwasawa decomposition is given by sp(n, R) = t ⊕ a ⊕ n where a = {(H, 0, 0)|H diagonal} n = {(A, B, 0)|A upper triangular with 0 on the diagonal , B symmetric} It can be shown that p

  is an isometry between the Siegel upper half space and the Siegel disk. Let C = I -iI I iI . The application g ∈ Sp(n, R) → CgC -1 identifies the set of isometries of H n and of D n . Thus, it can be shown that a matrix g = A B A B ∈ Sp(n, C) acts isometrically on D n by g.Z = (AZ + B)(AZ + B) -1

  ||u|| = d(p, exp p (u)) where ||.|| is the Euclidean distance associated to the local scalar product and d(., .) is the Riemannian distance. The corresponding Lebesgue measure of T p D n is noted Leb p . Let exp * p (Leb p ) denote the push-forward measure of Lep p by exp p . The function θ p defined by:

Figure 3 :

 3 Figure 3: Estimation of the density of 6 coefficients Ω k under rainy conditions. The expression of the used kernel is K(x) = 3 π (1-x 2 ) 2 1 x<1 . Densities are rescaled for visual purpose.

Figure 2 :

 2 Figure 2: Estimation of the density of 6 coefficients Ω k under clear conditions. The expression of the used kernel is K(x) = 3 π (1-x 2 ) 2 1 x<1 . Densities are rescaled for visual purpose.

We can now prove theorem 3.1. Let X be a random variable as in theorem 3.1. Following the notations of the theorem and the lemma, let

U is open, relatively compact and contains C. Let (M , g ) be as in the lemma. Let f and f be the kernel density estimators defined on M and M respectively. Theorem 3.1 of [1] provides the desired results for f . for r ≤ r inj , the support and the values on the support of f and f coincide. Thus the desired result also hold for f .


J F Semblat

Modeling Seismic Wave Propagation and Amplification in 1D/2D/3D Linear and Nonlinear Unbounded Media

Keywords: earthquake engineering, engineering seismology, wave propagation, numerical modeling

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

Modeling seismic wave propagation

Many various numerical methods are available to model seismic wave propagation and we will discuss them first. Afterwards, the issue of seismic wave amplification (site effects) in both linear (weak motion) and nonlinear (strong motion) ranges will then be examined. To analyze seismic wave propagation in 2D or 3D geological structures, various numerical methods are available (Fig. 1):

 the finite difference method is accurate in elastodynamics but is mainly adapted to simple geometries [START_REF] Bohlen | Accuracy of heterogeneous staggered-grid finite-difference modeling of rayleigh waves[END_REF][START_REF] Frankel | A Three-Dimensional Simulation of Seismic Waves in the Santa Clara Valley, California, from a Loma Prieta Aftershock[END_REF][START_REF] Moczo | 3D Heterogeneous Staggered-Grid Finite-Difference Modeling of Seismic Motion with Volume Harmonic and Arithmetic Averaging of Elastic Moduli and Densities[END_REF][START_REF] Virieux | P-SV Wave Propagation in Heterogeneous Media: Velocity-Stress Finite-Difference Method[END_REF]),  the finite element method is efficient to deal with complex geometries and numerous heterogeneities (even for inelastic constitutive models [START_REF] Bonilla | Computation of Linear and Nonlinear Site Response for Near Field Ground Motion[END_REF]) but has several drawbacks such as numerical dispersion and numerical damping [START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF][START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF][START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF], Semblat, 2000a[START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF] and (consequently) numerical cost in 3D elastodynamics,  the spectral element method has been increasingly considered to analyse 2D/3D wave propagation in linear media with a good accuracy due to its spectral convergence properties [START_REF] Chaljub | Spectral Element Analysis in Seismology[END_REF][START_REF] Faccioli | Spectral Domain Decomposition Methods for the Solution of Acoustic and Elastic Wave Equations[END_REF][START_REF] Komatitsch | The Spectral Element Method: An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures[END_REF],  the boundary element method allows a very good description of the radiation conditions but is preferably dedicated to weak heterogeneities and linear constitutive models [START_REF] Beskos | Boundary elements methods in dynamic analysis: Part II (1986-1996[END_REF][START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF][START_REF] Dangla | A plane strain soil-structure interaction model[END_REF], 2005[START_REF] Sanchez-Sesma | Seismic Response of Three-Dimensional Alluvial Valleys for Incident P, S and Rayleigh Waves[END_REF], Semblat 2008, 2000b). Recent developments have been proposed to reduce the computational cost of the method especially in the high frequency range [START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF], 2009[START_REF] Fujiwara | The fast multipole method for solving integral equations of three-dimensional topography and basin problems[END_REF],  the finite volume method was recently developed in the field of elastodynamics [START_REF] Glinsky-Olivier | A finite-volume method for the 2D seismic wave propagation[END_REF],  the Aki-Larner method which takes advantage of the frequency-wavenumber decomposition but is limited to simple geometries [START_REF] Aki | Surface Motion of a Layered Medium Having an Irregular Interface due to Incident Plane SH Waves[END_REF][START_REF] Bouchon | A Boundary Integral Equation-Discrete Wavenumber Representation Method to Study Wave Propagation in Multilayered Media Having Irregular Interfaces[END_REF]),  the scaled boundary finite element method which is a kind of solution-less boundary element method [START_REF] Wolf | The scaled boundary finite element method[END_REF],  other methods such as series expansions of wave functions [START_REF] Liao | A series solution and numerical technique for wave diffraction by a three-dimensional canyon[END_REF][START_REF] Sanchez-Sesma | Diffraction of elastic waves by three-dimensional surface irregularities[END_REF]. Furthermore, when dealing with wave propagation in unbounded domains, many of these numerical methods raise the need for absorbing boundary conditions to avoid spurious reflections. Since each method has specific advantages and shortcomings (Table I), it is consequently often more interesting to combine two methods to take advantage of their peculiarities. It is for instance possible to couple FEM and BEM [START_REF] Aochi | A complete BIEM-FDM-FEM simulation of an earthquake scenario[END_REF][START_REF] Dangla | A plane strain soil-structure interaction model[END_REF][START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF]) allowing an accurate description of the near field (FEM model including complex geometries, numerous heterogeneities and nonlinear constitutive laws) and a reliable estimation of the far-field (BEM involving accurate radiation conditions). The main advantage of the method in elastodynamics is the accurate description of the infinite extension of the medium through exact radiation conditions [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF], Semblat 2008). In contrast to other discretization methods, there is no need to consider absorbing boundary conditions. The Boundary Element Method involves singular integrals [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF]. There are three main kinds of singularities: (i) weak singularity (i.e. in the ordinary Riemann sense), (ii) strong singularity (i.e. in the Cauchy principal value sense) or (iii) hyper-singularity (i.e. in the Hadamard finite part sense). Strong and hyper-singular integrals have to be converted to regular ones in the regularization of the BEM formulations [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF][START_REF] Sladek | Singular integrals and boundary elements[END_REF]). Indeed, weak singularity is not treated by regularization. However, from the point of view of numerical integrations, one should devote a great attention to the evaluation of these integrals because standard integration quadratures fail in accuracy [START_REF] Beskos | Boundary elements methods in dynamic analysis: Part II (1986-1996[END_REF][START_REF] Dangla | A plane strain soil-structure interaction model[END_REF]. Therefore each type of singularity has to be treated by appropriate techniques [START_REF] Bonnet | Regularization of the displacement and traction BIE for 3D elastodynamics using indirect methods[END_REF][START_REF] Dangla | A simple and efficient regularization method for 3D BEM: application to frequency-domain elastodynamics[END_REF][START_REF] Sladek | Singular integrals and boundary elements[END_REF]). It is noticed that the regularization can be performed either before or after the discretization, i.e. in the global or local (intrinsic) coordinate space, as observed in some papers mentioned above. A comprehensive review of BEM in dynamic analysis has been proposed by [START_REF] Beskos | Boundary elements methods in dynamic analysis: Part II (1986-1996[END_REF]. An example (pressurized cavity in a full space) showing the influence of the regularization process on the numerical solution is proposed in Fig. 2: the uncorrected solution is very far from the regularized one which matches well the analytical solution [START_REF] Dangla | A simple and efficient regularization method for 3D BEM: application to frequency-domain elastodynamics[END_REF].

The main drawback of the BEM is that it leads to full non symmetric matrix systems. For large models, the numerical cost and memory storage requirements may then by huge. Recent researches considered symmetric Galerkin boundary element methods more efficient for large numbers of unknowns [START_REF] Bonnet | Symmetric Galerkin boundary element methods[END_REF]. Current researches try to reduced the number of computations through Fast Multipole formulations initially developed in the field of physics [START_REF] Fujiwara | The fast multipole method for solving integral equations of three-dimensional topography and basin problems[END_REF][START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF], 2009). Such alternative approaches are discussed in the following. [START_REF] Dangla | A simple and efficient regularization method for 3D BEM: application to frequency-domain elastodynamics[END_REF].

Fast Multipole formulation for the Boundary Element Method

In the field of physics (Maxwell or Laplace equations), recent advances in boundary element methods lead to a very important decrease of the computational cost. Instead of point to point interaction as in classical boundary element methods, the fast multipole method [START_REF] Greengard | Accelerating fast multipole methods for the Helmholtz equation at low frequencies[END_REF][START_REF] Fujiwara | The fast multipole method for solving integral equations of three-dimensional topography and basin problems[END_REF] considers interactions between groups of points (cells centered on a multipole, Fig. 3) and hence avoids multiple computations of nearly identical terms corresponding to very close points. It is possible to apply this method at a single scale or even at various scales through a multilevel approach. The size of each cell around its multipole depends on the distance to the other cells: the larger the distance between cells, the larger the cells (Fig. 3). The computations of the singular integrals are then performed through this approach. They are split in an integral on the surface around the singularity and another integral on the complementary distant surface. The first one is estimated using classical regularization techniques, whereas the latter is computed with a fast multipole algorithm [START_REF] Greengard | Accelerating fast multipole methods for the Helmholtz equation at low frequencies[END_REF]. The advantage for the computational cost is very significant since it depends on N 2 for classical methods and on N or NlogN for the fast multipole method [START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF]. Furthermore, the computational cost is reduced for both the memory storage and the calculation time. This method allows the analysis of very large problems involving millions of unknowns on a singleprocessor PC (to compare to several tens of thousands previously). Current researches also investigate a fast multipole method well-adapted to the Helmholtz equation and even to elastodynamics [START_REF] Fujiwara | The fast multipole method for solving integral equations of three-dimensional topography and basin problems[END_REF][START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF], 2009). The fast multipole method then allows the computation of very large models considering a larger number of heterogeneities, a more realistic geometrical representation of geological structures (especially in 3D) as well as higher frequency values (detailed modeling of short wavelengthes amplification).

Modeling seismic wave amplification by the Boundary Element Method

Amplification of seismic waves in geological structures

The amplification of the seismic motion mainly occurs in alluvial deposits the characteristics of which (geometry, wave velocities) control the amplification process [START_REF] Bard | The two dimensional resonance of sediment filled valleys[END_REF], Dobry 1976, Duval 1998[START_REF] Sanchez-Sesma | Diffraction of elastic waves by three-dimensional surface irregularities[END_REF], 2000, Semblat 2000, 2003). Considering the Boundary Element Method, seismic wave amplification (i.e. site effects) can be analyzed numerically taking into account the infinite extension of the medium as well as complex geometrical features. The main phenomena that can be recovered by the BEM models are the following:

 amplification due to the (1D) velocity contrast between soil layers,  focusing effects due to complex layers or free-surface geometries,  basins effects (2D/3D) due the trapped surfaces waves in surficial layers. Significant wave scattering may also be found around strong topographic irregularities such as crests or hills slopes [START_REF] Bouckovalas | Numerical evaluation of slope topography effects on seismic ground motion[END_REF][START_REF] Paolucci | Amplification of earthquake ground motion by steep topographic irregularities[END_REF][START_REF] Reinoso | Three-dimensional scattering of seismic waves from topographical structures[END_REF][START_REF] Sanchez-Sesma | Diffraction of elastic waves by three-dimensional surface irregularities[END_REF], Semblat 2009). Such phenomena are called topographic site effects and may aso be studied by the Boundary Element Method. In both case (alluvial basins and topographic irregularities), BEM models in 2D or 3D allow the estimation of the amplification factor of the seismic motion. The BEM formulation can be considered in frequency-domain for long durations and narrow bandwidths [START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF][START_REF] Dangla | A simple and efficient regularization method for 3D BEM: application to frequency-domain elastodynamics[END_REF] or time-domain for shorter durations and larger bandwidths [START_REF] Gaul | A comparative study of three boundary element approaches to calculate the transient response of viscoelastic solids with unbounded domains[END_REF][START_REF] Jin | A 2D time-domain boundary element method with damping[END_REF][START_REF] Manolis | Boundary Element Methods in Elastodynamics[END_REF]). Other numerical methods may also be considered as discussed in previous sections [START_REF] Faccioli | Spectral Domain Decomposition Methods for the Solution of Acoustic and Elastic Wave Equations[END_REF][START_REF] Glinsky-Olivier | A finite-volume method for the 2D seismic wave propagation[END_REF][START_REF] Komatitsch | The Spectral Element Method: An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures[END_REF][START_REF] Moczo | 3D Heterogeneous Staggered-Grid Finite-Difference Modeling of Seismic Motion with Volume Harmonic and Arithmetic Averaging of Elastic Moduli and Densities[END_REF][START_REF] Virieux | P-SV Wave Propagation in Heterogeneous Media: Velocity-Stress Finite-Difference Method[END_REF]).

Topographic amplification of seismic waves

The analysis of wave scattering around topographic irregularities has been performed by various authors for canonical configurations [START_REF] Bouckovalas | Numerical evaluation of slope topography effects on seismic ground motion[END_REF][START_REF] Reinoso | Three-dimensional scattering of seismic waves from topographical structures[END_REF][START_REF] Sanchez-Sesma | Diffraction of elastic waves by three-dimensional surface irregularities[END_REF] or for actual 2D/3D topographies [START_REF] Paolucci | Amplification of earthquake ground motion by steep topographic irregularities[END_REF][START_REF] Semblat | Seismic Site Effects in a Deep Alluvial Basin: Numerical Analysis by the boundary element method[END_REF]. As shown in Fig. 4 in the case of Caracas [START_REF] Semblat | Seismic Site Effects in a Deep Alluvial Basin: Numerical Analysis by the boundary element method[END_REF], the maximum amplification Amax along a complex topography may lead to significant changes in the seismic wavefield. When comparing the top plot at frequency 0.25 Hz to the bottom one (0.35 Hz), these topographic site effects are nevertheless found to be less than those found in alluvial basins. Detailed analyses of such effects are proposed in [START_REF] Bouckovalas | Numerical evaluation of slope topography effects on seismic ground motion[END_REF][START_REF] Paolucci | Amplification of earthquake ground motion by steep topographic irregularities[END_REF], Semblat 2009). [START_REF] Bouckovalas | Numerical evaluation of slope topography effects on seismic ground motion[END_REF] proposed empirical laws to estimate the topographic effects for 2D steep crests. For 3D models, efficient numerical methods are needed. [START_REF] Paolucci | Amplification of earthquake ground motion by steep topographic irregularities[END_REF] studied an actual steep italian site using the Spectral Element Method. Since the development of Fast Multipole BEMs in elastodynamics allows to handle large scale 3D models, several large canonical 3D topographies have been recently studied [START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF]. As shown in Figure 5 for an ellipsoidal canyon and an oblique incident plane P-wave [START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF], the FM-BEM solution is in good agreement with previous solutions obtained by [START_REF] Reinoso | Three-dimensional scattering of seismic waves from topographical structures[END_REF] for both the vertical and horizontal motion components. Several other canonical configurations were defined in the framework of the QSHA research project (Quantitative Seismic Hazard Assessment) and are used to compare the efficiency and accuracy of various numerical methods (data available at http://qsha.unice.fr/). Current researches aim at investigating actual 3D topographies by the Fast Multipole Method. example in Caracas, Venezuela [START_REF] Semblat | Seismic Site Effects in a Deep Alluvial Basin: Numerical Analysis by the boundary element method[END_REF]. (data for various canonical configurations available at http://qsha.unice.fr/)

Amplification of seismic waves in alluvial basins

The amplification of seismic waves in alluvial basins is generally larger than that observed/computed around topographic irregularities [START_REF] Baise | The effect of shallow San Francisco bay sediments on waveforms recorded during the Mw 4.6 Bolinas California Earthquake[END_REF][START_REF] Beauval | Quantification of frequency-dependent lengthening of seismic groundmotion duration due to local geology: applications to the Volvi area (Greece)[END_REF][START_REF] Frankel | A Three-Dimensional Simulation of Seismic Waves in the Santa Clara Valley, California, from a Loma Prieta Aftershock[END_REF][START_REF] Sanchez-Sesma | Seismic Response of Three-Dimensional Alluvial Valleys for Incident P, S and Rayleigh Waves[END_REF], Semblat 2009[START_REF] Sommerville | Emerging art: earthquake ground motion[END_REF]). Since extensive field measurements were performed on Volvi EuroSeisTest, Greece [START_REF] Bard | Wave propagation in complex geological structures and their effects on strong ground motion[END_REF][START_REF] Beauval | Quantification of frequency-dependent lengthening of seismic groundmotion duration due to local geology: applications to the Volvi area (Greece)[END_REF], Chavez-Garcia 2000[START_REF] Pitilakis | Site effects : recent considerations and design provisions[END_REF]), BEM models were performed for this site (Semblat 2005). The numerical results were compared to numerous seismological field data from various earthquake recordings.

The Volvi basin is 6km wide and 250m deep. As shown in Figure 6, the amplification level changes from one frequency to the other. The influence on the amplification process of both the basin geometry and the soil layering were studied in (Semblat 2005). Since the geometry of the basin is complex, strong focusing effects occur leading to large amplification levels (around 10) at some peculiar frequencies (Figure 6). However, depending on the accuracy of the geological model, the description of the soil layering has a strong influence at higher frequencies (Figure 6, bottom right). It has been quantified for a simplified (2 layers) as well as detailed (6 layers) model of the Volvi profile in (Semblat 2005). The detailed analysis leads to larger amplitude trapped surface waves at higher frequencies.

F1=0.6Hz, A1=2.9 F2=0.8Hz, A2=9.5 F4=1.2Hz, A4=7.5 F6=2.4Hz, A6=7.3 Figure 6. Seismic wave amplification in the Volvi-EuroSeisTest basin, Greece (Semblat 2005): BEM simulations at different frequencies Fi and maximum related amplifications.

3 Modeling seismic wave propagation by the Finite/Spectral Element Method

Numerical wave dispersion

Considering the finite/spectral element method to model seismic wave propagation, as for other numerical methods, two types of numerical errors can be studied [START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF][START_REF] Chaljub | Spectral Element Analysis in Seismology[END_REF][START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF] : the relative period error and the algorithmic damping. The relative period error is different from one integration scheme to another [START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF]. For wave propagation problems, the relative period error appears in the estimation of wave-velocities and is called numerical dispersion. Seismic wave propagation in a numerical scheme then depends on the element size, the integration scheme, the element type, etc. This phenomenon is called numerical dispersion as a reference to physical dispersion leading to a dependence of the wave-velocity on frequency [START_REF] Deraemaeker | Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions[END_REF][START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF][START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF], Semblat 2000b). To model wave propagation phenomena by the finite element method or the finite difference method, one has to control the accuracy of the numerical scheme since the numerical error tends to increase during the propagation process.

Various theoretical works concern the analysis of the numerical error made in the estimation of the approximated wave number [START_REF] Deraemaeker | Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions[END_REF][START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF], Semblat 2009). For instance, it was shown that there is a cutoff frequency above which there is no propagation at all [START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF]. Depending on the frequency of the excitation, the numerical wave propagates slower or faster than the theoretical solution. It is then necessary to analyze numerical dispersion of seismic waves and precisely quantify the numerical error.

Higher order finite elements vs spectral elements

Numerical wave dispersion is also influenced by the order of the shape functions of the finite elements. Higher order finite elements are for instance known to have a very good precision for elastic-plastic computations. In acoustics and elastodynamics, several theoretical works propose analytical expressions to estimate the numerical dispersion [START_REF] Deraemaeker | Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions[END_REF][START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[END_REF]. Some comparisons between lower order and higher order finite elements to model wave propagation were recently proposed by [START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF]. In Figure 7, the results from 1D FEM simulations (Semblat 2000b) are plotted for various shape functions orders.

The same number of degrees of freedom is considered in each case. The linear finite elements (Figure 7, top) are shown to have strong numerical dispersion (numerical wave velocity artificially increased). For quadratic finite elements (Figure 7, middle), the accuracy is satisfactory when compared to the theoretical delays (vertical dashed lines). Finally, for the higher order finite elements (Figure 7, bottom), the numerical wave dispersion is very low. These results show that, for the same number of DOFs, the accuracy of higher order finite element for wave propagation simulations is much better. In addition to "classical" higher order finite elements, spectral finite elements have increasingly been studied since their spectral convergence is of great interest for wave propagation simulations [START_REF] Chaljub | Spectral Element Analysis in Seismology[END_REF][START_REF] Faccioli | Spectral Domain Decomposition Methods for the Solution of Acoustic and Elastic Wave Equations[END_REF][START_REF] Komatitsch | The Spectral Element Method: An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures[END_REF]. The spectral elements are generally considered for high order only (4 to 8) since the difference with classical finite elements is not significant below. An evaluation of the cost effectiveness (and accuracy) of both Chebyshev spectral and p-finite elements has been proposed by [START_REF] Dauksher | An evaluation of the cost effectiveness of Chebyshev spectral and p-finite element solutions to the scalar wave equation[END_REF].

Modeling linear attenuation

The analysis of seismic wave propagation in attenuating media, such as soils and rocks, raises the need for suitable attenuation models. One of the most common damping model in structural dynamics is the Rayleigh formulation giving a damping matrix as a linear combination of stiffness and mass matrices [START_REF] Hughes | Duality and Unified Analysis of Discrete Approximations in Structural Dynamics and Wave Propagation: Comparison of p-method Finite Elements with k-method NURBS[END_REF]. The damping ratio is then minimum for a given frequency and infinite for zero and infinite frequencies (Semblat 2009). In the field of structural dynamics, the main advantage of the Rayleigh formulation is that it leads to a diagonal damping marix in the real eigenmodes base [START_REF] Chopra | Dynamics of Structures, Theory and Applications to Earthquake Eng[END_REF].

For wave propagation phenomena, the Rayleigh damping formulation is equivalent to a particular rheological model: generalized Maxwell model [START_REF] Semblat | Rheological interpretation of Rayleigh damping[END_REF]. For weak to moderate damping values, this equivalence allows an easier and explicit estimation of the two coefficients of the Rayleigh formulation from experimental behaviour parameters. Other types of "matrix based" damping formulations are available such as the Caughey formulation [START_REF] Chopra | Dynamics of Structures, Theory and Applications to Earthquake Eng[END_REF], Semblat 2009) corresponding to a generalization of the Rayleigh formulation and allowing complex damping-frequency dependences.

Various types of mechanical or rheological viscoelastic formulations were also proposed. [START_REF] Kjartansson | Constant Q wave propagation and attenuation[END_REF] investigated frequency constant attenuation ("CQ") models for wave propagation problems. [START_REF] Day | Numerical simulation of wavefields using a Padé approximant method[END_REF] considered a Padé approximant to simulate wave propagation in attenuatingmedia. [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF] proposed a rheological model leading to a nearly constant attenuation ("NCQ") in a given frequency range. [START_REF] Carcione | Time-domain Modeling of Constant-Q Seismic Waves Using Fractional Derivatives[END_REF] studied a damping formulation based on fractional derivatives and [START_REF] Moczo | On the rheological models used for time-domain methods of seismic wave propagation[END_REF] discussed time domain simulations using rheological models. All these approaches are formulated in the linear range and need, for time domain computations, to deal with memory variables which may be costly. As discussed in the following, nonlinear constitutive laws are often necessary in the case of strong earthquakes.

Radiation conditions at infinity

For finite differences, finite or spectral element methods, it is necessary to take into account the radiation conditions through efficient techniques to avoid spurious wave reflections at the mesh boundaries. It can be done considering absorbing boundary conditions or infinite elements for purely solid [START_REF] Chadwick | Diffraction of short waves modelled using new mapped wave envelope finite and infinite elements[END_REF] as well as multiphase media [START_REF] Modaressi | Paraxial Approximation for Poroelastic Media[END_REF]. It is nevertheless difficult for very heterogeneous media [START_REF] Chammas | Characterization of heterogeneous soils using surface waves : homogenization and numerical modeling[END_REF].

Absorbing layer methods rather try to model the wave attenuation in a thick layer at the medium boundaries (Semblat 2009). They have been increasingly developed in recent years since their efficiency appears good in various configurations. They are generally known as "Perfectly Matched Layer" methods. Various types of PML formulations have recently been proposed:  Classical PML [START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF]: the wavefield propagating in the layer is damped according to a onedimensional amplitude decrease law similar to that of an attenuating medium:

  i x x x ) ( ~  (1)
This complex valued law leads to an amplitude decrease of the numerical wave along the x-axis. As illustrated in Figure 8, grazing incidences are related to small horizontal wavenumbers and a limited amplitude decrease will be obtained. This uniaxial law is thus not sufficient to have an optimal absorbing effect in all configurations.

 Filtering PML [START_REF] Festa | PML absorbing boundaries[END_REF][START_REF] Festa | The Newmark scheme as a Velocity-Stress Time staggering: An efficient PML for Spectral Element simulations of elastodynamics[END_REF]: for shallow numerical models, the classical PML formulation may also amplify surface waves (Figure 8). A filtering PML formulation (F-PML) has been proposed by [START_REF] Festa | The Newmark scheme as a Velocity-Stress Time staggering: An efficient PML for Spectral Element simulations of elastodynamics[END_REF] to avoid this problem:

c i x x x       ) ( ~(2)
 Multidirectional PML [START_REF] Meza-Fajardo | A Nonconvolutional, Split-Field, Perfectly Matched Layer for Wave Propagation in Isotropic and Anisotropic Elastic Media: Stability Analysis[END_REF]: another alternative to ensure the numerical stability of PMLs (e.g. grazing incidences) is the multidirectional PML formulation (M-PML) recently proposed by Meza Fajardo and Papageorgiou [START_REF] Meza-Fajardo | A Nonconvolutional, Split-Field, Perfectly Matched Layer for Wave Propagation in Isotropic and Anisotropic Elastic Media: Stability Analysis[END_REF]. The original uniaxial law is generalized in the multi-axial case through the following system:

               z i z z y i y y x i x x x z x y x x ) ( ) ( ) ( 1 ~1 ~1 ~      (3)
It is thus possible to choose the attenuation vector  in order to optimize the absorbing effect depending on the wave type and incidence. Another interest of this formulation is its numerical stability for strongly anisotropic media (Figure 8).

Modeling strong ground motion and non linear effects

To model strong ground motion, it is necessary to combine seismological tools to geotechnical methods involving nonlinear constitutive laws [START_REF] Heuze | Estimating site-specific strong earthquake motion[END_REF]. The nonlinear soil models should reproduce one of the main features of soil dynamic/cyclic behaviour: the shear modulus decreases and the hysteretic damping increases with strain amplitude [START_REF] Seed | Moduli and damping factors for dynamic analyses of cohesionless soils[END_REF]). Such models are generally considered in the framework of the finite element method and often involve Masing type cyclic behaviour [START_REF] Bonilla | Computation of Linear and Nonlinear Site Response for Near Field Ground Motion[END_REF]. The influence of plasticity, confining pressure and pore-pressure build up can also be analyzed with more complex models [START_REF] Hujeux | A double memory model with multiple mechanisms for cyclic soil behaviour[END_REF][START_REF] Gyebi | Finite element analysis of viscoplastic soils with Q-factor[END_REF][START_REF] Lade | Elasto-plastic stress-strain theory for cohesionless soil with curved yield surfaces[END_REF][START_REF] Loret | Growth and decay of acceleration waves in non-associative elastic-plastic fluidstaturated porous media[END_REF][START_REF] Mellal | A simplified numerical approach for nonlinear dynamic analysis of multilayered media[END_REF][START_REF] Park | Soil damping formulation in nonlinear time domain site response analysis[END_REF][START_REF] Prevost | A simple plasticity theory for frictional cohesionless soils[END_REF]. It is then possible to have a deeper insight into the liquefaction processes. More and more strong motion data from well-controlled bore-holes are available and allow comparisons with nonlinear computations.

In the last decades, linear (viscoelastic) equivalent models were also extensively used to have a simplified description (i.e. few parameters) of the shear modulus decrease and the hysteretic damping increase [START_REF] Schnabel | SHAKE-A computer program for equation response analysis of horizontally layered sites[END_REF]. Simplified models are interesting since they may allow combined seismological/geotechnical computations to simultaneously take into account basin effects at large scales (2D/3D) and nonlinear local effects. Recent researches proposed new simplified models avoiding some drawbacks of equivalent models such as frequency independence [START_REF] Assimaki | Model for dynamic shear modulus and damping in granular soils[END_REF][START_REF] Kausel | Seismic simulation of inelastic soils via frequency-dependent moduli and damping[END_REF]). As shown in Fig. 9 (left), this nonlinear frequency dependent model combines a strain spectrum with the classical G() and () curves. [START_REF] Delépine | A simplified non linear model to analyze site effects in alluvial deposits[END_REF][START_REF] Delépine | Nonlinear viscoelastic wave propagation: an extension of Nearly Constant Attenuation (NCQ) models[END_REF] proposed another alternative: a simple constitutive model generalizing the classical NCQ (Nearly Constant) attenuation model [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF] in the nonlinear range. This extended NCQ model leads to decreasing shear modulus and increasing damping/attenuation. This model allows the analysis of seismic waves propagation from strong earthquakes and may be used at large scales (e.g. alluvial basins). Finally, a crucial issue is the availability of the field/laboratory data to feed and validate such models.

Conclusion

This paper presents recent results in various research fields related to seismic wave propagation in geological structures. New methods or major improvements of existing methods (staggered grid FD, Finite Volumes, Fast Multipole-BEM) are promising ways to better model seismic wave propagation, amplification and interactions in 3D surficial geological structures. The major advances of current works also concern improved treatments of the radiation conditions for spectral elements/finite elements (absorbing layer methods such as F-PML or M-PML).

Another key point is strong seismic motion: recent simplified nonlinear constitutive models allow the analysis of seismic wave propagation from strong earthquakes at large scales. Nevertheless, several topics have not been discussed in this paper: fault mechanics and the influence of fluids [START_REF] Aochi | A complete BIEM-FDM-FEM simulation of an earthquake scenario[END_REF], soil-structure interaction [START_REF] Kham | Seismic Site-City Interaction: Main Governing Phenomena through Simplified Numerical Models[END_REF][START_REF] Sextos | Inelastic dynamic analysis of RC bridges accounting for spatial variability of ground motion, site effects and soil-structure interaction phenomena. Part 1: Methodology and analytical tools[END_REF]) and site-city interaction [START_REF] Bielak | Earthquake ground motion and structural response in alluvial valleys[END_REF][START_REF] Clouteau | Modifications of the ground motion in dense urban areas[END_REF][START_REF] Groby | Simulation of seismic response in a city-like environment[END_REF], Semblat 2008[START_REF] Tsogka | Seismic response of a set of blocks partially embedded in soft soil[END_REF][START_REF] Wirgin | Effects of buildings on the duration and amplitude of ground motion in Mexico city[END_REF], soil-pile interaction and liquefaction [START_REF] Brennan | Effectiveness of vertical drains in mitigation of liquefaction[END_REF][START_REF] Takahashi | Liquefaction-induced large displacement of pile-supported wharf[END_REF], etc.

Acknowledgements

The author would like to thank all his colleagues and students who contributed to some of the results presented herein and the related papers listed in the References.

Figure 1 .

 1 Figure 1. Various numerical methods to model seismic wave propagation: (a) the Finite Difference Method, (b) the Finite Element Method, (c), the Spectral Element Method, (d) the Boundary Element Method, (e) the Discrete Element Method, (f) the Finite Volume Method.

Figure 3 .

 3 Figure 3. Comparison of the principles of the standard BEM and the Fast Multipole BEM.

Figure 4 .

 4 Figure 4. Topographic (top) vs statigraphic (bottom) seismic wave amplification:example in Caracas, Venezuela[START_REF] Semblat | Seismic Site Effects in a Deep Alluvial Basin: Numerical Analysis by the boundary element method[END_REF].

Figure 5 .

 5 Figure5. Validation of the Fast Multipole BEM from a 3D configuration studied by[START_REF] Reinoso | Three-dimensional scattering of seismic waves from topographical structures[END_REF]:diffraction of an oblique plane P-wave by a semi-ellipsoidal canyon[START_REF] Chaillat | A multi-level Fast Multipole BEM for 3-D elastodynamics in the frequency domain[END_REF].(data for various canonical configurations available at http://qsha.unice.fr/)

Figure 7 .

 7 Figure 7. Numerical wave dispersion for various finite element orders: linear (top), quadratic (middle) and fourth-order (bottom) (Semblat 2000b).

Figure 8 .

 8 Figure 8. Various critical configurations for the efficiency/stability of PML absorbing conditions (from Semblat, 2009).

Figure 9 .

 9 Figure 9. Recent improvements of the equivalent linear model for the analysis of strong seismic motion in soils: (a) Frequency dependent nonlinear model (Kausel 2002); (b) Nonlinear "extended NCQ" model (Delépine 2009).

Table I .

 I Features of various numerical methods for wave propagation modeling.

	Method	FDM	FEM	SEM	BEM	DEM	FVM
	Continuous	yes	yes	yes	yes	no(/yes)	yes
	Constitutive	lin./NL	lin./NL	lin./NL	linear	lin./NL	lin.
	Radiation	abs. bound.	abs. bound.	abs. bound.	intrinsic	abs. bound.	abs. bound.

2 Modeling seismic wave propagation by the Boundary Element Method

2.1 Interest of the Boundary Element Method ("BEM")