Mohamad Jouni

Mauro Dalla Mura

Pierre Comon

Classification of Hyperspectral Images as Tensors Using Nonnegative CP Decomposition

Keywords: Hyperspectral Imagery, Morphological Profiles, Tensor Decomposition, Scene Classification

A Hyperspectral Image (HSI) is an image that is acquired by means of spatial and spectral acquisitions, over an almost continuous spectrum. Pixelwise classification is an important application in HSI due to the natural spectral diversity that the latter brings. There are many works where spatial information (e.g., contextual relations in a spatial neighborhood) is exploited performing a so-called spectral-spatial classification. In this paper, the problem of spectral-spatial classification is addressed in a different manner. First a transformation based on morphological operators is used with an example on additive morphological decomposition (AMD), resulting in a 4-way block of data. The resulting model is identified using tensor decomposition. We take advantage of the compact form of the tensor decomposition to represent the data in order to finally perform a pixelwise classification. Experimental results show that the proposed method provides better performance in comparison to other state-of-the-art methods.

Introduction

Hyperspectral Imaging is one of the most important tools in remote sensing. A hyperspectral image (HSI) is typically a three-way 1 block of data acquired when many (two-way) images are taken over hundreds of almost-continuous spectral bands and stacked together forming a so called hypercube. HSI is employed in several fields of applications such as astronomy, Earth and planetary observation, monitoring of natural resources, precision agriculture and biomedical imaging. One of its most common usages is classification, i.e., a thematic discrimination of different types of objects present in the image. In other words, the goal of classification is to assign a label to each pixel belonging to the same thematic class according to some of its characteristics (e.g., its spectrum, the homogeneity of its spatial neighborhood or the shape of a region it belongs to). In the case of HSI, the image can be seen as a mapping from a couple of positive integers (x, y) that correspond to horizontal and vertical positioning, to a vector of positive real values that correspond to the radiance measured by the sensor or the reflectance of the scene. Hence, a direct way to classify a HSI is to take the spectral bands as features and the pixels as samples. However, this approach raises some problems. First, in terms of spectral resolution, the acquisition of HSI induces high dimensionality in the spectra, which means that the samples are put in a high dimensional feature space. Second, by only considering the spectral information can be a limiting factor since spatial information (e.g., contextual relations, size and shape of the objects) is not exploited. For this reason, many works aimed at incorporating spatial information based on the pixels' neighborhood, so that each pixel in the HSI has features that include both spectral and spatial information, and this is known in the literature as spectral-spatial classification; such previous works could be found in [START_REF] Bruzzone | A multilevel context-based system for classification of very high spatial resolution images[END_REF][START_REF] Duarte-Carvajalino | Multiscale representation and segmentation of hyperspectral imagery using geometric partial differential equations and algebraic multigrid methods[END_REF][START_REF] Li | Spectral-spatial hyperspectral image segmentation using subspace multinomial logistic regression and markov random fields[END_REF][START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF][START_REF] Fauvel | A spatial-spectral kernel-based approach for the classification of remote-sensing images[END_REF][START_REF] Fauvel | Advances in spectral-spatial classification of hyperspectral images[END_REF][START_REF] Ghamisi | A survey on spectral-spatial classification techniques based on attribute profiles[END_REF][START_REF] Ghamisi | Frontiers in spectral-spatial classification of hyperspectral images[END_REF] with different approaches.

One way to incorporate spatial information in HSI is by using tools defined in the framework of Mathematical Morphology (MM). A well-known example is the morphological profile (MP) [START_REF] Atli Benediktsson | Classification and feature extraction for remote sensing images from urban areas based on morphological transformations[END_REF], which is built by stacking the results of a set of openings and closings based on geodesic reconstruction applied on each band of the HSI (or on each component obtained after dimensionality reduction as in [START_REF] Atli Benediktsson | Classification of hyperspectral data from urban areas based on extended morphological profiles[END_REF]) and has been successfully employed in remote sensing [START_REF] Dalla Mura | The evolution of the morphological profile: From panchromatic to hyperspectral images[END_REF][START_REF] Ghamisi | A survey on spectral-spatial classification techniques based on attribute profiles[END_REF]. Attribute Profiles (AP) [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF] are another example, in which a scalar image (e.g., a band of a HSI) is filtered with attribute thinnings and thickenings computed on component trees [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF] or self-dual attribute filters based on the the tree-of-shapes [START_REF] Cavallaro | Remote sensing image classification using attribute filters defined over the tree of shapes[END_REF].

Since HSI is by nature a three-way block of data, applying transformations based on MM results in a four-way block of data, which would be of high rank. Since we care about having low rank data, we merge the first two ways that correspond to pixel locations. As a result we obtain a three-way block of data where each dimension could be described by pixels, spectral bands, and spatial features (e.g., by means of morphological transformations), respectively.

In this paper we aim at modeling the 3-way hyperspectral data using tensors, and address the problem of pixelwise classification using tensor tools. Our contribution is inspired by [START_REF] Velasco | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF], where an Additive Morphological Decomposition (AMD) is applied on HSI, resulting in a 4-way block of data, which is then modeled by tensor means. In [START_REF] Velasco | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF], the corresponding tensor is dealt with using Tensor-PCA (TPCA) in order to apply dimensionality reduction of the data, which is then passed down to classification. In this paper, our motivation comes from the fact that most tensors enjoy a unique Canonical Polyadic (CP) decomposition [START_REF] Comon | Tensors: a brief introduction[END_REF], and we take advantage of how it represents the data for classification purposes. We sum up our contribution as follows:

-CP decomposition is used as a way to represent the pixel data in a compact form, which is more direct and easy to deal with and works thanks to a low dimensionality of the feature space. Moreover, CP decomposition can be adapted to incorporate the nonnegativity of its entries, which leads to a decomposition that can have a physical meaning.

-The proposed method shows promising results, when comparing it to TPCA approach.

The rest of the paper is organized as follows. In Section 2 we briefly talk about adding a new dimension to the data set based on MM. In Section 3 we describe CP decomposition, and explain how it can be a suitable way to directly represent the pixel data in a low dimensional feature space ready for classification.

Experiments and results are reported in Section 4. Finally, a conclusion is drawn in Section 5 and some future work is outlined.

Mathematical Morphology

Let us consider { Φi , Φ i } i=1...m , a set of extensive and anti-extensive operators, respectively, where m refers to the size of the filter and the greater m the coarser the filter effect. The MP of a grayscale image (i.e., with scalar pixels value) I using the latter set of operators is defined as: In the case of HSI, or more generally in the case of multivariate images, a direct extension to MP is the Extended Morphological Profiles (EMP) [START_REF] Atli Benediktsson | Classification of hyperspectral data from urban areas based on extended morphological profiles[END_REF]. The concept is the same: one simply concatenates the MP of each image along a fourth way and the result is provided as an input to the classifier. Along the same lines, works in the literature extended this concept to other possible rearrangements or derivations, and an example that we adopt in our experiments is the AMD [START_REF] Velasco | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF].

As explained in [START_REF] Velasco | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF], first we suppose that we have any image I decomposed through AMD such that:

I = Φm ( Φm-1 (I)) + Φ m (Φ m-1 (I)) 2 + m i=1 R - i -R + i 2 = S + m i=1 R i , (2)
The term S is the structure and contains the unfiltered objects of the decomposition, while the terms {R i } i=1...m are the residuals and they contain elements that are filtered out according to the parameters used. Each element of the set {S, R i } i=1...m is a 3-way data set, and at the end they are stacked to form a 4-way data set.

Unlike [START_REF] Velasco | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF], we prefer to build the morphological way using the set {S, R - i , R + i } i=1...m in order to preserve the nonnegativity of the data, as a result we have a 4-way hyperspectral data set

D = [S, R - 1 , R + 1 , . . . , R - m , R + m ],
which is represented by an I 1 × I 2 × J × (2m + 1) array, I 1 × I 2 is the number of pixels, J is the number of spectral bands, and 2m + 1 = K is the number of terms in the morphological decomposition.

We note that the case study is not restricted to this kind of application, one could also build a 4th way using attribute thinning and thickenings as done in attribute profiles [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF] for example.

Tensor Decomposition

In this section, we explain how to process by tensor decomposition the data obtained with AMD. D is of high rank, this can be explained starting with the fact that usually a 2D image is full rank which is high, and a tensor built upon this image will have a rank that is at least as much as that of the image. So we group the first two ways so that the first dimension becomes I 1 × I 2 = I, and D becomes an I × J × K three-way array.

CP Decomposition

A third order tensor D r is referred to as decomposable if it can be written as

D r def = a r ⊗ b r ⊗ c r , (3)
A tensor Y can always be decomposed into a sum of decomposable tensors D r :

Y = R r=1 λ r D r , (4)
The tensor rank, usually denoted by R, is the least number of terms required such that the CP decomposition (4) is exact. Suppose we have a 3-way tensor Y ∈ R I×J×K , each decomposable term D r in the CP decomposition can be expressed as in [START_REF] Li | Spectral-spatial hyperspectral image segmentation using subspace multinomial logistic regression and markov random fields[END_REF], where the set of vectors {a r } r=1...R ∈ R I , {b r } r=1...R ∈ R J , and {c r } r=1...R ∈ R K form the columns of the so-called factor matrices, A ∈ R I×R , B ∈ R J×R , and C ∈ R K×R . Hence another way to write equation 4 is, with an obvious notation:

Y = (A, B, C) • Λ (5)
The diagonal tensor Λ governs the interaction between the components (columns) of the factor matrices such that components of different r indices do not interact. R can then be seen as directly related to the degrees of freedom in the decomposition.

When Y is nonnegative, it is often desirable to constrain factor matrices to be also nonnegative. In this case, the rank is called nonnegative rank and sometimes denoted R + . Having nonnegative factors is essential for the sake of physical interpretation of the results where negative values have no meaning. Consequently, the goal is to minimize the cost function:

1 2 Y -(A, B, C).Λ 2 (6) s.t. A 0, B 0, C 0
An interesting aspect is that each of the factor matrices represents one way of the tensor. For example, an element y i,j,k is decomposed into R + components and the decomposed information is mapped to the i-th row of A, the j-th row of B, and the k-th row of C.

In the case of HSI, A corresponds to the way of pixels, B corresponds to the way of spectral bands, and C corresponds to the way based on MM. Accordingly each row of A represents a pixel, that of B represents a spectral band, and that of C is defined by a scale or a structure. Thanks to the nonnegativity constraints and especially after normalizing the columns of the factor matrices, one could assume that columns of B represent spectral signatures and those of C represent a combination of scale domination along the set of structuring elements being chosen, and finally each column of A could be seen as an abundance map, that could be reshaped to show a grayscale image, with respect to information received at B and C, very similar to the case of spectral unmixing (but for two dimensions only) -except that unmixing is presently not our concern.

In this paper, we take advantage of the low-rank aspect of our CP decomposition, i.e. how it represents the data in a more compact form with almost no loss of information (depending on how exact the decomposition is). First, instead of dealing with a multidimensional data set of large dimensions, we are able to deal with matrices of much smaller size, as explained earlier (since a small rank permits compression). Second, since our approach is pixelwise classification, we mainly think of A as the matrix to deal with. Theoretically speaking, similar pixels should have the same composition of coefficients with respect to the components in A, so we take the components (columns) as features and the rows as samples in the classification.

High-Order Singular Value Decomposition (HOSVD)

Because of the huge dimensions of our tensor, directly applying CP decomposition would be computationally very demanding. One solution to make the decomposition possible is to pass by a pre-processing step by compressing the original tensor into one with smaller dimensions without losing any information. This can be done by providing the minimal compression size for each way such that the information is preserved. In that sense, let's suppose a tensor Y with dimensions I × J × K, where I JK (which is our case), the compressed tensor can have a size of JK × J × K whenever its rank satisfies R ≤ JK.

G = Y •1 V T 1 •2 . . . •N V T N ; return G and V 1, . . . , V N
This can be explained as follows. First let us denote the d-mode unfolding of Y by Y (d) , which is the matrix flattening of the tensor obtained by reordering the elements with respect to the d-th index. Each of these unfoldings has a matrix rank that essentially appears in its economic SVD. More precisely, after applying the SVD on each unfolding, we get three right singular (orthogonal) matrices that form the basis of the HOSVD factor matrices:

Y (d) = U d Σ d V T d (7) ∀d = {1, 2, 3}
These matrices, denoted by V 1 , V 2 , and V 3 have dimensions I × R 1 , J × R 2 and K × R 3 respectively, with R 1 ≤ JK, R 2 ≤ J, and R 3 ≤ K. Hence they can be used to compress the original tensor into smaller dimensions without any loss of information. In fact, this is how HOSVD is computed by finding orthogonal matrices as the basis of the decomposition, and the compressed tensor G is found by projecting the original one onto the basis of factor matrices. HOSVD of Y can be written as follows:

Y = G • 1 V 1 • 2 V 2 • 3 V 3 = (V 1 , V 2 , V 3 ) • G (8)
where • d indicates a summation over the dth tensor index. Subsequently, we shall compress only the first mode, so that R 1 < JK, but R 2 = J and R 3 = K.

As a consequence, we may choose V 2 and V 3 both as Identity matrices. The algorithm that we implement to compress an N -way tensor using HOSVD is described in Algorithm 1; way c is a vector containing the ways that we want to compress and dim c is a vector containing the compressed dimension size.

Alternating Optimization -Alternating Direction Method of Multipliers (AO-ADMM)

In order to compute the CP decomposition and cope with the constraints of nonnegativitity and compression together, the algorithm AO-ADMM [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF] is adopted for its efficiency and flexibility especially with the different kinds of constraints.

The CP decomposition can be computed via the Alternating Least Squares algorithm (ALS), where one factor matrix is updated alternatingly at a time by fixing the others and minimizing the convex least squares cost function:

H d = argmin H d 1 2 Y (d) -W H T d 2 (9) ∀d = 1 . . . N,
where H d is the factor matrix corresponding to the d-th way of the tensor, and W is the Khatri-Rao product of all the factor matrices excluding H d . The optimization problem [START_REF] Atli Benediktsson | Classification and feature extraction for remote sensing images from urban areas based on morphological transformations[END_REF] needs to be modified when constraints on the factor matrices are imposed, say through regularization functions r(H d ):

H d = argmin H d 1 2 Y (d) -W H T d 2 + r(H d ) (10) ∀d = 1 . . . N,
As the name suggests, AO-ADMM works by alternating between the optimization problems of each factor matrix alone while fixing the others. For example in our case, the three factor matrices (H 1 , H 2 , H 3 ) are recommended to have nonnegative entries, while only the first one is concerned with compression. AO-ADMM allows to pass different kinds of constraints and parameters for each mode-decomposition and alternates accordingly.

Algorithm 2 Alternating Optimization (Least Squares loss)

Require: Y, H1, . . . , HN , U [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF] for the update of µ) until some termination criterion is reached (number of iterations) Normalize the columns of the factor matrices and store the weights in Λ; return H1, . . . , HN , Λ It is important to note that the tensor fed in AO-ADMM is the compressed one that is computed as a result of Section 3.2, which allows to have negative entries, and this presents a certain inconsistency with the fact of having nonnegativity constraints on the decomposition. For this reason we present a solution with ADMM inspired by [START_REF] Cohen | Fast decomposition of large nonnegative tensors[END_REF]. So our minimization problem goes as follows: arg min [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF] for µ and Hf )

H1,H2,H3,Λ 1 2 G -(H 1 , H 2 , H 3 ).Λ 2 (11)
s.t. V 1 H 1 0, H 2 0, H 3 0, Λ 0 Algorithm 3 ADMM of a mode d Require: Y , W , H, U , k, µ, , constraint, imax, V H and U are already initialized; G = W T W ; ρ = trace(G)/k; Calculate L from Cholesky decomposition such that G + (ρ + µ)Ik = LL T ; F = W T Y ; Hf = H; repeat H ← (L T ) -1 L -1 (F + ρ(H + U ) T + µHf T ); (See
H ← proximity(constraint, HT ,U ,V ); (refer to Algorithm 4)

U ← U + H -HT ; Update r and s; (refer to [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF] for the updates of r and s) until (r < and (s is undefined or s < )) or (imax > 0 and i ≥ imax) return H and U Algorithm 4 Proximity Update of H in ADMM Require: constraint, Ht, U , V switch (constraint) case Nonnegativity: H ← (Ht -U ) + ; case Compression and Nonnegativity:

Hu ← V (Ht -U ); Hu ← H + u ; H ← V T Hu; end switch return H
This minimization is executed by ADMM as follows, for d ∈ {2, 3}:

H1 ← (W T W + ρI) -1 (W T G (1) + ρ(U + H1) T ) H1 ← V T 1 max(0, V 1( HT 1 -U )) U ← U + H1 - HT 1      (12)
Hd

← (W T W + ρI) -1 (W T G (d) + ρ(U + H d ) T ) H d ← max(0, HT d -U ) U ← U + H d - HT d      (13)
Of course, the calculation of (W T W + ρI) -1 should be done once for all outside the loops, as explained in the pseudo-codes Algs. 2, 3, and 4.

Results

Description of the dataset

Our dataset is described by a real hyperspectral image, University of Pavia, as shown in Fig. 2 in true colors (by choosing the bands 53, 31, and 8 as Red, Green, and Blue channels respectively). It is an image acquired by the ROSIS sensor with geometric resolution of 1.3 meters. The image has 610 × 340 pixels and 103 spectral bands. The dataset contains a groundtruth image that consists of 9 classes: trees, asphalt, bitumen, gravel, metal sheets, shadows, self-blocking bricks, meadows, and bare soil. There are also 42776 labeled pixels available as test set and 3921 pixels available as training set, the latter is fixed at that.

Classification

For the classification part we use Support Vector Machines (SVM), which has proved useful in the application of pixelwise classification of HSI. Practically SVM is carried out through the open source machine learning library provided by Libsvm [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF]. After decomposing the tensor using CP decomposition, we use the first factor matrix as the set of data to be classified by SVM.

Results and discussion

In the following experiments, we fix the morphological transformations on structuring elements whose neighborhoods are defined arbitrarily by disk shapes with varying sizes of radii: 1, 6, 11, 16, 21, and 26 pixels, this accounts to m = 6 structuring elements, thus a morphological decomposition of K = 13 components. We note that considering the whole set of {S, R - i , R + i } i=1...6 without dimensionality reduction yielded a classification with better accuracy than that of the set of {S, R i } i=1...6 (96.60% to 93.06% respectively) under the same AMD parameters.

AMD. The HSI of University of Pavia is loaded in its 3rd-order form, I ∈ R 610×340×103 . I is then decomposed into a morphological structure, going from a dataset of one tensor to a set of tensors (S, R -'s and R + 's, forming a 4th-order dataset) representing information of the scene with different scales:

I → D = {S, R - 1 , R + 1 , . . . , R - 6 , R + 6 } such that I = S + 6 i=1 R - i -R + i 2
. As we mentioned in Section 3, D is vectorized by merging the first two modes resulting in a 3-way tensor of dimensions 207400× 103 × 13.

Tensors. Now that the data is modeled as a 3rd-order tensor, we seek to decompose it using CP decomposition. The latter is carried out by taking the compressed tensor as input, the factor matrices are initialized randomly (absolute value of i.i.d. standard Gaussian distribution), and only the first factor matrix is compressed after that. The rank is set to different values followed by different initializations of factor matrices.

Classification. The accuracy of the classification is related to the reconstruction error of the CP decomposition since it is important that factor matrix H 1 represents the data in a shape as good as possible. Two factors play major roles in the reconstruction error: the number of iterations and the input rank of the decomposition. Hence we try different values of the two variables and record both the reconstruction error and the result of the overall accuracy as seen in Table 1. Per-class and average accuracy is also recorded for two methods in Table 2. As discussed in Section 3.1, H 1 is passed to the classifier since it is considered as the representative of pixel data in compact low-dimensional form.

Method

Rec. Error OA AMD + TPCA + SVM (PC= [START_REF] Atli Benediktsson | Classification of hyperspectral data from urban areas based on extended morphological profiles[END_REF][START_REF] Fauvel | A spatial-spectral kernel-based approach for the classification of remote-sensing images[END_REF] Results. We compare our results to the method of [START_REF] Velasco | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF]. Looking at Table 1, we can see that the better the reconstruction error in the CP decomposition, the better the overall accuracy of the classification. The accuracy shows promising results for our proposed method. It is worth to stress the influence of rank and/or number of iterations on the reconstruction error. But at some point the error decreases very slowly as more iterations are run, towards what seems to be a limit, which is mainly influenced by the value of the input rank. Figure 2 shows images of the training set, the test set, and the results of various classifications; parameters are mentioned in the caption.

Conclusion

In the framework of pixelwise classification of hyperspectral images, an extraneous diversity is built from hyperspectral data thanks to Mathematical Morphology. Then the use of tensor decomposition allows to take into account both spectral and spatial information of the scene. Experimental results using a real image show that the CP decomposition is a promising way to represent the data in a low dimensional feature space and improve on classification accuracy. In future works, we plan to investigate other potential interests of the CP decomposition, like spectral unmixing of HSI.

  MP(I) := {Φ m (I), . . . , Φ 1 (I), I, Φ1 (I), . . . , Φm (I)}, (1) in which the following ordering relation holds MP(I) i ≤ MP(I) i+1 for i = 1, . . . , 2m and with MP(I) i denoting the i-th component of the MP. A MP using Opening and Closing by reconstruction is shown in Figure 1.

Fig. 1 .

 1 Fig. 1. Morphological profile of a portion of the 8th spectral band of the HSI of University of Pavia using 3 different SE; disks with sizes [1,6,11].

Algorithm 1 COMPRESS

 1 Require: Y, wayc, dimc for d ∈ wayc do Unfold Y into Y (d) such that the mode d takes the second way of the matrix; Compute the right singular matrix from the SVD of Y (d) , denoted by V d ; Truncate the columns of V d by dimc(d); end for for d / ∈ wayc do V d = I dimc(d) ; (Identity matrix) end for

Fig. 2 [

 2 A] image in true colors, [B] available training set, [C] test set, [D-F] classification results: [D] AMD+TPCA(15)+SVM (91.09%), [E] AMD+CPD (itr=50,R=20)+SVM (91.17%), [F] AMD+CPD (itr=50,R=30)+SVM (94.08%)

  1, . . . , U

N H1, . . . , HN are initialized; U 1, . . . , U N are initialized to zero matrices; for d = 1, . . . , N do Store the different unfoldings Y d ; end for repeat for d = 1, . . . , N do W = j =d Hj; update H d and U d using Algorithm 3; end for update µ if necessary; (refer to

Table 1 .

 1 Records of various tests in terms of Reconstruction Error (Rec. Error) (in case of CPD) and Overall Accuracy (OA).

	)	-	91.09

Table 2 .

 2 Records of Per-class, Overall and Average Accuracies.

The number of ways of an array refers to the number of its indices. A HSI is typically a three-way array of dimensions I1 × I2 × J, where I1 and I2 are space dimensions (i.e. pixels) and J denotes the number of spectral bands.


Mohamad Jouni

Mauro Dalla Mura

Pierre Comon

HYPERSPECTRAL IMAGE CLASSIFICATION USING TENSOR CP DECOMPOSITION

Keywords: Remote Sensing Image, Mathematical Morphology, Attribute Profiles, Tensor Decomposition, Scene Classification

Image classification has been at the core of remote sensing applications. Optical remote sensing imaging systems naturally acquire images with spectral features corresponding to pixels. Spectral classification ignores the spatial distribution of the data which is becoming more relevant with the development of spatial resolution sensors, and many works aim to incorporate spatial features based on neighborhood through for example, Mathematical Morphology (MM). Additionally, one could stack multiple morphological transformations of the image resulting in a highly complex block of data. Since classification is a tool that requires a matrix of samples and features, and simply stacking the different sets of features can lead to the problem of high dimensionality, we propose a way to create a matrix of low dimensional feature space by modeling the data as tensors and thanks to Canonical Polyadic (CP) decomposition. Experiments on real image show the effectiveness of the proposed method.

INTRODUCTION

Classification of pixels has become a core subject in remote sensing applications allowing each pixel to be associated to a semantic label. For this procedure to take place in a pixelwise approach, it is necessary to have a set of samples and features, the former representing the pixels of a scene and the latter representing distinctive characteristics acquired from the scene.

Thanks to optical remote sensing systems, spectral features of the scene can naturally be acquired so that each pixel is associated to a discrete spectral signature. When the spectral acquisition spans a very high number of narrow bands, such images are referred to as hyperspectral images (HSI), represented as data cubes of 3rd-order.

With the improvements in spatial resolution of acquired images, it has become easier to identify relatively small objects and more relevant to include spatial information as features. In fact, many works such in [START_REF] Bruzzone | A multilevel context-based system for classification of very high spatial resolution images[END_REF][START_REF] Duarte-Carvajalino | Multiscale representation and segmentation of hyperspectral imagery using geometric partial differential equations and algebraic multigrid methods[END_REF][START_REF] Li | Spectralspatial hyperspectral image segmentation using subspace multinomial logistic regression and markov random fields[END_REF][START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF][START_REF] Fauvel | A spatial-spectral kernel-based approach for the classification of remote-sensing images[END_REF][START_REF] Fauvel | Advances in spectral-spatial classification of hyperspectral images[END_REF][START_REF] Ghamisi | A survey on spectral-spatial classification techniques based on attribute profiles[END_REF][START_REF] Ghamisi | Frontiers in spectral-spatial classification of hyperspectral images[END_REF] aimed at incorporating spatial features based on neighborhood information in addition to spectral ones, towards what is referred to as spectral-spatial classification. One way to account to spatial information in classification is by extracting features obtained by filtering an image with operators based on Mathematical Morphology.

Creating a new mode based on spatial MM features normally yields a block of data that is of order four: i.e., two spatial dimensions (image rows and columns), one to spectral and one to spatial features. Furthermore, one might want to stack multiple types of spatial transformations based on different concepts (e.g., size and shape) in a fifth mode, making it a data block of order five. As such, merging the first two modes that account to positioning will still give an array that is at least of order three, so in terms of classification, one way to deal with it is to reshape the dataset into a matrix, but this can still suffer from high dimensionality when the size of the spatial features is large, also this technique ignores the multimodal structure of the data and merges them together as one.

One natural approach to deal with the multi-way structure of the data is by using tensor analysis [START_REF] Comon | Tensors: a brief introduction[END_REF]. Tensors are data structures in multi-linear algebra, they are normally the extension of vectors and matrices into higher orders. Tensors are represented by multi-dimensional (or multi-way) arrays, hence they are a natural selection and a direct way to handle the problem of high-order features and account to the interactions between the different modes.

In short, we aim to model high-order-feature hyperspectral data in the framework of tensors and exploit their tools for the sake of classification. In fact, tensor modeling of HSI for classification has been explored previously in [START_REF] Velasco-Forero | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF] where Tensor-Principal Component Analysis (TPCA), a highorder extension of PCA, is used as a method to maintain the spectral-spatial multi-modal structure and reduce the dimensionality of the data.

In this paper, we build a high-order tensor by stacking multiple types of spectral and spatial features separated in different modes, and we choose to jointly handle these features by taking advantage of the CP decomposition. We show through experiments on real image that CP decomposition is a good way to represent the pixel data directly in a "compact matrix form" and thanks to a low dimensionality of the de- composition, leading to a low dimensionality of the feature space despite the high complexity of the data. Moreover, it is possible to incorporate nonnegative constraints into the CP decomposition, which allows to spectral unmixing in the future. Our approach to model multiple types of spatial features is inspired by the work of [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF] in Extended Multi-Attribute Profiles (EMAP), where multiple types of morphological attribute operators are considered in the classification. In [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF], the different Extended Attribute profiles (EAP) are concatenated along the third mode of the HSI. Instead, we keep the different EAP separated by stacking them as a new mode in the tensor. The proposed method incorporating multiple types of features shows promising results in terms of accuracy of classification compared to other methods in the literature.

From here on the paper is organized as follows. Section 2 introduces the proposed method and the challenges within. Section 3 shows the experiments that were carried on. Finally, a conclusion is drawn in Section 4.

THE PROPOSED METHOD

In this section we briefly introduce the main aspects of the proposed method, which is shown in the flowchart of Figure 1. There are originally two main steps before entering classification, the first one is creating the spatial features based on MM and the second one is the CP decomposition. However, due to challenges induced by high dimensionality of the spectral bands and the high number of pixels, two sub-steps are needed. First, we reduce the number of bands with very minimal loss as a preprocessing step, which is usually done using PCA, but since we care about conserving the nonnegativity aspect of the data (for spectral unmixing which is out of the scope of this paper), we perform the reduction using Nonnegative Matrix Factorization (NMF). Second, before the CP decomposition step, we compress the first mode of the tensor in order to reduce its size, which helps a lot with computing the CP decomposition in an efficient way.

Nonnegative Matrix Factorization Preprocessing

First let us denote a HSI by M ∈ R I1×I2×J , where I = I 1 I 2 is the number of pixels and J is the number of spectral bands. We are interested in reducing the number of bands regardless of the pixel positioning, so M is reshaped in lexicographic order into M ∈ R I×J by merging the first two modes. NMF decomposes a matrix M ∈ R I×J into the product of two nonnegative matrices such that M = W H T , where W ∈ R I×R and H ∈ R J×R , and R is the number of components in the decomposition. Our interest is that W can be thought as a way to represent pixel data with R spectral features while conserving nonnegativity for reasons mentioned earlier, where essentially R < J. W is then reshaped back to 3rd-order W ∈ R I1×I2×J .

Spatial Features Using Mathematical Morphology

As a case example, the focus in this paper will be on Extended Attribute Profiles (EAP) [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF] for their effectiveness in extracting distinctive spatial features. Nevertheless, the proposed approach can also extend to different approaches for incorporating spatial features.

The attribute profile (AP) of a 2D image with scalar values is formed by concatenating the attribute transformations obtained using successive values of the attribute parameter, the result is a 3rd-order data block. The EAP of a multivariate image is formed when the AP of the single 2D images are concatenated, the result is a 3rd-order data block, where the spectral and spatial features are stacked in one mode. Our take on EAP is a bit different, we prefer to separate the spectral and spatial modes such that the result of an EAP is a 4thorder block of data denoted as T ∈ R I1×I2×J×K , and K is the number of attribute transformations.

Furthermore, we compute more than one type of EAP of W. These EAP are then stacked in parallel to form a fifth mode which can be seen as an acquisition of different types of EAP. Finally, the fifth-order data block is reshaped by merg-ing the first two modes to form a fourth-order data block denoted by Y ∈ R I×J×K×L , which is ready to be dealt with using tensor decomposition, and L denotes the number of EAP being computed.

Compressed Tensor Decomposition

Tensor Y can be seen as a high-order-feature hyperspectral data. Seeing as the first mode of Y is the one of pixels, i.e. samples, and seeing as Y is a fourth order tensor, this means that each pixel has a third-order tensor of features which is indeed complex and heavy to deal with for the sake of classification, and here comes the importance of tensor decomposition. Tensor decomposition is a tool used to break the high-order complex structure of tensors into multiple simpler structures of matrices, while accounting to the intrinsic interactions between the elements of the tensor.

One of the most important of such tools is the CP decomposition, which can be seen in two different ways. To begin with, an N th-order tensor D ∈ R I1ו••×I N is referred to as decomposable if it can be written as the outer product of

N vectors v n ∈ R In s.t. n ∈ {1 . . . N }: D def = v 1 ⊗ . . . ⊗ v N .
On one hand CP decomposition allows to decompose a tensor into the sum of rank-1 decomposable tensors:

Y = R r=1 λ r D r , (1)
where the tensor rank, denoted by R, is the least number of terms such that the CP decomposition is exact. On another hand, it allows to decompose an N -th order tensor into N matrices, usually referred to as "factor matrices" which are formed by stacking each set of R vectors as columns, and a diagonal core tensor Λ, formed by filling the diagonal with values of λ r , under the following notation:

Y = (V 1 , . . . , V N ) • Λ. (2)
CP decomposition can adopt the nonnegativity constraints in its computation in order to yield results that can be physically interpretable, and this will be more relevant for future works. Hence, one aspect that appears in computing the CP decomposition is that each factor matrix V n is directly related to the n-th mode of Y. In our case of application, we are interested that V 1 directly describes the first mode of pixels by means of R components which could be seen as a compact form of features, and one could say that the complex structure of pixels in Y boils down to row vectors of R elements, and due to redundancies that are present in Y, R could be relatively of a small value that is acceptable when feeding V 1 to the classifier.

It is worth noting here that since the number of pixels is huge compared to the other dimensions, a compressed nonnegative CP decomposition inspired by [START_REF] Cohen | Fast decomposition of large nonnegative tensors[END_REF] and adopted by the algorithm Alternating Optimization -Alternating Direction Method of Multipliers (AO-ADMM) [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF] is carried out.

EXPERIMENTS AND RESULTS

Description of the dataset

Our dataset is described by a real hyperspectral image, University of Pavia. It is an image acquired by the ROSIS sensor with geometric resolution of 1.3 meters. The image has 610×340 pixels and 103 spectral bands. The dataset contains a groundtruth image that consists of 9 classes: trees, asphalt, bitumen, gravel, metal sheets, shadows, self-blocking bricks, meadows, and bare soil. There are also 42776 labeled pixels available as test set and 3921 pixels available as training set, the latter is fixed at that.

Experimental set-up

NMF. We start with the HSI such that M ∈ R 610×340×103 and the number of pixels is 610 × 340 = 207400. As a preprocessing step in order to reduce the number of 103 bands, we find that the NMF with 30 components yields a reconstruction error of 0.89% where W 30 ∈ R 610×340×30 , and that with 40 components yields a reconstruction error of 0.65% where W 40 ∈ R 610×340×40 . We use both datasets for the following steps.

EAP features. We use three different types of EAP as mentioned in [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF], the area of the regions (EAP a ) with values of λ a = {100, 500, 1000, 5000}, the diagonal of the box bounding the region (EAP d ) with values of λ d = {10, 25, 50, 100}, and the standard deviation of the grey-level values of the pixels in the regions (EAP s ) with values of λ s = {20, 30, 40, 50}, bringing a wide diversity of features to the image. Having 4 values of λ corresponds to K = 9.

First we create the 3 EAP for each of W 30 and W 40 corresponding to L = 3, then we stack them in a 5th order tensor to which the outputs are Y 30 ∈ R 207400×30×9×3 and Y 40 ∈ R 207400×40×9×3 respectively. Denote by Extended Multi-Attribute Profile (EMAP) as the stacking of multiple types of EAP.

Compressed nonnegative CP decomposition. Now that we have the 4th-order tensor, we pass it to the CP decomposition step. As mentioned earlier, this part is adopted by the algorithm AO-ADMM with constraints of nonnegativity on all the modes in addition to constraint of compression only on the first mode. The factor matrices are randomly initialized using the absolute value of i.i.d. standard Gaussian distribution, and only the first factor matrix is compressed after that.

Classification. For all the classifications, we use a Random Forest (RF) classifier with the number of trees set to 100. Table 1 shows some results of various conducted tests, where we compare the proposed method to others in the literature, namely PCA as carried out in [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF] for EMAP and TPCA as carried out in [START_REF] Velasco-Forero | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF]. The comparisons are based on the Overall Accuracy (OA) and the Average Accuracy (AA) of the classification. For all the conducted tests, the same set of attribute parameters was chosen.

Results and discussion

Looking at Table 1, we can see that the best results are achieved by computing the EMAP on the NMFs and applying CPD afterwards. The proposed method shows improved accuracy compared to the competing techniques, PCA and TPCA. On a second note, we see that classifying the components obtained by NMF yields better accuracy than classifying the components obtained by PCA. In addition to good accuracy, obtaining nonnegative features using NMF has an advantage of conserving nonnegativity and extracting physical meaning from the data compared to features obtained using PCA.

CONCLUSION

In this paper, creating multiple spatial diversities thanks to mathematical morphology was explored for the sake of spectral-spatial classification of HSI. The corresponding high-order data was handled using tensor decomposition by jointly accounting to spectral and spatial features. Some challenges like reducing the high dimensionality of the spectral bands while conserving the nonnegativity of the data was dealt with using NMF, which also showed some advantages over choosing PCA. Experimental results using a real HSI showed that the CP decomposition is a promising candidate to represent high-order feature data in low dimensional feature space with good accuracy. In future works, we plan to dig more into the properties of nonnegativity and investigate other potential interests of the CP decomposition, like spectral unmixing of HSI.

Fig. 1 :

 1 Fig. 1: Flow chart of the proposed method.

Table 1 :

 1 Records of various tests in terms of Overall Accuracy (OA) and Average Accuracy (AA). The number of features for each test appears in the second columns. The best results are in bold and the second best are underlined.

	Method	Features	OA	AA
	PCA	4	68.34 78.01
	NMF 30	30	73.79 85.10
	NMF 40	40	77.74 86.33
	PCA + EMAP [4]	108	91.12 90.36
	EMAP + TPCA [10]	100	90.60 91.65
	NMF 30 + EMAP + CPD 30	30	93.96 92.50
	NMF 40 + EMAP + CPD 30	30	94.20 95.85
	NMF 40 + EMAP + CPD 40	40	94.76 94.38


Rodaina Aboul
email: rodaina.aboulhosn@3sr-grenoble.fr

Hosn • Luc

Nadia Benahmed

Bruno Chareyre

R Aboul Hosn

L Sibille

Discrete numerical modeling of loose soil with spherical particles and interparticle rolling friction

Keywords: Discrete element method, Rolling friction, Calibration, Compaction technique

Discrete numerical simulations were carried out to reproduce experimental results obtained on loose cohesionless soil samples subjected to triaxial tests. Periodic boundary conditions were adopted and 3D spherical discrete elements were chosen. However, to overcome excessive rolling of such an oversimplified particle's shape, contact rolling resistance was taken into consideration. The influence of both the elastic and the plastic local parameters is discussed. It is shown that the plastic macroscopic behavior of the granular assembly depends only on the plastic parameters at the microscopic scale, and mainly on the plastic rolling moment reflecting the particle's shape. Moreover, a procedure to obtain an initial density, ranging from loose to dense samples, is proposed by adding adhesion at contacts during the preparation phase. Finally, a calibration procedure is proposed to reproduce experimental results and the limitations of the model are discussed.

Introduction

Most studies based on the discrete (or distinct) element method, pioneered by Cundall and Strack [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], use simplified particle shapes as disks (in 2D) [START_REF] Ting | Discrete numerical model for soil mechanics[END_REF] or spheres (in 3D) [START_REF] Thornton | Numerical simulations of deviatoric shear deformation of granular media[END_REF]. Yet, such shapes lead to excessive rolling as compared to real granular materials [START_REF] Iwashita | Micro-deformation mechanism of shear banding process based on modified distinct element method[END_REF][START_REF] Bardet | Observations on the effects of particle rotations on the failure of idealized granular materials[END_REF][START_REF] Kozicki | Numerical simulations of sand behaviour using DEM with two different descriptions of grain roughness[END_REF]. To overcome this setback, complex shapes as polygons [START_REF] Szarf | Polygons vs. clumps of discs: A numerical study of the influence of grain shape on the mechanical behaviour of granular materials[END_REF] and polyhedrons can be used [START_REF] Lee | Simulation of triaxial compression tests with polyhedral discrete elements[END_REF]. In fact, such shapes can be complicated and computationally expensive in detecting contacts and calculating forces and torques [START_REF] Radjaï | Discrete-element modeling of granular materials[END_REF]. Another alternative solution is to use aggregates or clumps of disks or spheres [START_REF] Salot | Influence of relative density on granular materials behavior: DEM simulations of triaxial tests[END_REF]. However, this requires a larger number of spheres and the contribution of the non-convex surface of clumps becomes difficult to evaluate. Thus, one solution is to fix the particle's rotation [START_REF] Calvetti | A numerical investigation of the incremental behavior of granular soils[END_REF] or to take into account the transfer of a moment between elements in the local constitutive law as pioneered by Iwashita and Oda [START_REF] Oda | Experimental micromechanical evaluation of strength of granular materials: Effects of particle rolling[END_REF][START_REF] Iwashita | Micro-deformation mechanism of shear banding process based on modified distinct element method[END_REF][START_REF] Oda | Study on couple stress and shear band development in granular media based on numerical simulation analyses[END_REF], or more recently by Estrada [START_REF] Estrada | Comparison of the effects of rolling resistance and angularity in sheared granular media[END_REF] to reflect the grain's angularity. In 3D, moment transfer laws can be defined with respect to both twisting (around the contact normal) and rolling (orthogonal to the contact normal) interparticle relative rotation. Only the transfer of rolling moment is addressed in this paper.

Previous studies investigated the influence of the contact rolling resistance properties on the macroscopic behavior [START_REF] Plassiard | A spherical discrete element model: calibration procedure and incremental response[END_REF][START_REF] Widuliński | Numerical simulations of triaxial test with sand using DEM[END_REF][START_REF] Wang | Simulation of triaxial response of granular materials by modified DEM[END_REF][START_REF] Kozicki | Study of some micro-structural phenomena in granular shear zones[END_REF]. All the results showed that increasing the rolling elastic stiffness causes an increase in the macroscopic internal friction angle, but has little influence on the dilatancy [START_REF] Plassiard | A spherical discrete element model: calibration procedure and incremental response[END_REF][START_REF] Widuliński | Numerical simulations of triaxial test with sand using DEM[END_REF]. On the other hand, [START_REF] Wang | Simulation of triaxial response of granular materials by modified DEM[END_REF] showed that when the rolling stiffness is high enough, its influence on both the peak friction angle and the dilatancy angle becomes negligible. However, in all cases, it is suggested to calibrate the contact rolling stiffness (in addition to the plastic rolling moment) to match the macroscopic internal friction angle. This leads to a surprising cross dependency of microscopic elastic parameters and macroscopic plastic ones. Thereby, the first objective of this paper is to show that such a cross dependency can be avoided.

In addition to the particle's shape and the contact law, the mechanical behavior of a granular material is highly dependent on its initial density. This latter is often reached numerically by playing with the contact friction angle during the compaction phase which proves to be powerful in preparing medium dense to dense samples [START_REF] Tong | Pore-scale flow simulations: Model predictions compared with experiments on bi-dispersed granular assemblies[END_REF][START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF][START_REF] Scholtés | Micromechanics of granular materials with capillary effects[END_REF], but not very loose ones. Indeed, discrete numerical models of loose soils are quite rare in the literature. Hence, the second objective of the paper is to address this issue by adding contact adhesion during the compaction phase of the numerical sample. In this way, we mimic the moist tamping technique used in the laboratory to reconstitute loose soil samples.

Before dealing with these objectives, the contact law including rolling resistance is presented, and the dependency of the stability of explicit numerical integration scheme on the rolling stiffness is defined. Then, a calibration procedure is proposed in the last section by comparing numerical simulations to experimental results of triaxial tests performed on Camargue sandy soil [START_REF] Nguyen | Etude expérimentale du comportement instable d'un sable silteux: Application aux digues de protection[END_REF][START_REF] Benahmed | An experimental investigation into the effects of low plastic fines content on the behaviour of sand/silt mixtures[END_REF]. Finally, the limitations of the defined model in reproducing the dependency of the volumetric strains on the mean pressure are discussed.

Note that the effects of some parameters as damping, strain rate or inertial number and tangential stiffness are already known from previous works [START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF][START_REF] Combe | Origines géométrique du comportement quasistatique des assemblages granulaires denses: étude par simulations numériques[END_REF] and are not addressed in this paper.

2 The discrete numerical model

Inter-particle contact law

Let us consider two spheres in contact with radii R 1 and R 2 . In the normal and tangential directions to the contact plane, the interaction law is very classical and characterized by constant normal and tangential (or shear) stiffnesses, K n and K s , respectively, and a contact friction angle, ϕ c , such that:

F n = K n δ n n (1)
∆ F s = -K s ∆ U s with || F s || ≤ || F n || tan ϕ c ( 2
)
where n is the normal to the contact plane, δ n , the overlapping distance between spheres, and U s , the relative tangential (or shear) displacement at the contact point.

Only compressive normal forces are modeled, and the contact is lost as soon as the overlap, δ n , vanishes.

The rolling resistance at the contact is defined by the rolling stiffness, K r , and the coefficient of rolling friction, η r (it is commonly called friction because of the mathematical form it takes although it doesn't directly involve the friction between two surfaces). Thus, the rolling moment, M r , acting against the relative rolling rotation of particles, θ r , is expressed as:

∆ M r = -K r ∆ θ r with || M r || ≤ || F n || η r min(R 1 , R 2 ) (3)
∆ θ r is computed as the tangential component of the incremental relative rotation, ∆ θ, of the contacting spheres.

For spheres with incremental rotations, ∆ ω 1 and ∆ ω 2 respectively, ∆ θ = ∆ ω 2 -∆ ω 1 and:

∆ θ r = ∆ θ -(∆ θ . n) n (4)
Note that Equation (3) defining the contact resistance to rolling is similar to Equation (2) representing the sliding resistance due to dry friction. Therefore, contact resistance to rolling can be interpreted as a rolling friction.

To avoid any dependency of macroscopic elastic properties on particle size, contact stiffnesses are defined from a stiffness modulus, E c , and dimensionless shear and rolling coefficients, α s and α r , respectively:

K n = 2 E c R 1 R 2 R 1 + R 2 ; K s = α s K n ; K r = α r R 1 R 2 K s . (5)
Finally, adhesive normal and tangential forces, F A n and F A s , can be added to the contact law (used here during the sample compaction only) and defined from an adhesive stress σ A such that:

F A n = F A s = σ A [min(R 1 , R 2 )] 2 (6)
Then, the contact presents a resistance to a tensile normal force as long as:

F n . n > -F A n , (7)
and a maximum shear force against sliding is given as:

|| F s || ≤ || F n || tan ϕ c + F A s . (8)

Stability condition of the explicit integration scheme

The stability condition of the explicit centered finite difference scheme, integrating motion equations of particles, can be derived by determining a stiffness matrix K, linking displacement and rotation of an inertial particle to force and torque [START_REF]PFC3D theory and background manual[END_REF][START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF]. K refers to a particular contact network, thus changes during the simulation. The full details of the derivations are given in the Appendix. In order to understand the role of each stiffness in the stability condition, we consider the following simplifying assumptions (not assumed in the numerical model):

1. all spheres are identical (same size and same inertia) 2. all contacts have the same stiffness values 3. the stiffness tensors are isotropic.

The small movements u and θ, in translation and in rotation, respectively, of a particle relative to an initially stable configuration are described by (see Appendix):

ü = - N c K n (1 + 2α s ) M u, (9) θ
= - N c K n (5α s + 5α rn + 2.5α tw ) M θ (10)
where α rn and α tw are respectively the rolling and twisting coefficients (defined in Appendix), K n is the normal stiffness, M is the mass of the particle and N c is the number of contacts formed by the particle. Numerators in the right hand sides of the uncoupled equations of motion 9 and 10 represent the global stiffnesses of the particle in translation and rotation respectively. The critical time step will be fixed by the highest global stiffness. So how the different kinds of contact stiffnesses can affect the order of magnitude of the timestep?

if the normal stiffness is dominant (α k ≪ 1 for all k), the stability condition is controlled by the translational motion (eq. ( 9)); -if α s is not negligible, precisely as soon as α s ≥ 1/3, the rotational motion governs the stability condition (eq. ( 10)); -if α rn ≥ 1/3 or α tw ≥ 2/3 (even if α s remains small), the rotational motion governs the stability condition; -if α k ≈ 1 for all k, the timestep imposed by the rotational motion is less than half of the one corresponding to translational motion if α k > 1 for at least one k, estimating the timestep from k n alone leads to overestimate the maximum allowed timestep significantly.

When the numerical scheme is unstable due to the translational motion, particles are usually ejected away from the granular mass after few timesteps. In most cases it cannot remain unnoticed by the user, so that mistakes on the timestep are immediately detected. On the other hand, preliminary numerical experiments (not reported here) showed that the consequences of numerical instabilities for rotational modes may not be obvious at the macroscale. Only more detailed examinations may reveal spurious relative motion at some contacts and high angular velocities of few particles. It is thus highly required to use a robust and validated procedure for choosing the timestep in the general case. The numerical results may otherwise suffer from (apparently) unexplained bugs. The fact that this problem has not been clearly described previously may explain why the stiff limit for contact moments was not shown by most previous authors.

Macroscopic constitutive behaviour

Numerical triaxial compression tests were performed using YADE software [34] to investigate the effects of some local mechanical parameters on the constitutive behavior of the model. It was made up of 10,000 spherical discrete particles whose radius is equal to 0.014 and 0.026 to avoid crystallization. Periodic boundary conditions were adopted. The periodic cell was formed as a parallelepipedic block filled with a cloud of spheres (i.e. assembly of nonoverlapping particles). The non-overlap constraint necessarily leads to rather loose clouds. Since this is only the starting point of the compaction, it has no consequences on the equilibrium state after compaction. Then, isotropic compaction was applied until reaching the required confining pressure P c . Finally, the triaxial compression was performed. During the compaction phase, the stress is controlled to reach a confined pressure of 80 kPa. In the triaxial compression phase, the inertia number is chosen 1.65 x10 -4 to avoid inertial effects. In this work, the timestep was chosen to be 80% of the maximum value defined by the stability condition (eq. 20 in appendix). Moreover, no damping is considered in the triaxial phase. Typical simulated responses can be seen in Figures 3 and5. Table 1 summarizes the parameters used for these studies.

Elastic local parameters and plastic macroscopic properties

In this section, the influence of elastic local parameters, such as the normal and the rolling stiffnesses, on the macroscopic shear strength will be examined. Shear Table 1: Summary of the parameters used in each case of the parametric studies.

strength at failure is characterized by the peak friction angle, ϕ p , reached at the peak of the deviator stress, q, along a drained compression (constant effective confining stress), while, ϕ 0 , characterizes the shear strength at the residual state (or the critical state).

We investigate first the effect of the normal contact stiffness represented via the dimensionless stiffness κ:

κ = E c 2P c (11)
The response shown in Figure 1, where no rolling resistance was introduced (case 1 in Table 1), demonstrates that the peak friction angle, ϕ p , decreases with the normal contact stiffness particularly for κ < 50. Similar results appear in [START_REF] Roux | Discrete numerical simulation and the mechanical behavior of granular materials[END_REF][START_REF] Roux | IUTAM-ISIMM Symposium on Mathematical Modeling and Physical Instances of Granular Flow[END_REF]. This is a remarkable disagreement with the experimental results on sand, in which the peak friction angle decreases for decreasing κ (i.e. for increasing the contact stiffness). This experimental trend may result simply from the dependency of the contact friction on the normal force [START_REF] Suhr | On the effect of stress dependent interparticle friction in direct shear tests[END_REF]. In numerical models, the contact friction coefficient is usually constant. Thus the decreasing trend of the internal friction is explained by the decreasing number of contacts. Moreover, Figure 1 shows that for a higher dimensionless stiffness, where the overlap becomes negligible, there is still a scattering among the values of ϕ p , but quite low (about ±0.5 • around the mean value). Thereby, for ϕ c = 30 • , plastic failure seems not significantly dependent of the contact normal stiffness for a value of κ greater than 50.

A more extended parametric study was performed concerning the parameters involved in the contact resistance to rolling. The values of the parameters defined in Table 1, case 2, were chosen in order to test different pairs of α r and η r (note that for isotropic shape particles, η r should be ≤ 1 [START_REF] Iwashita | Micro-deformation mechanism of shear banding process based on modified distinct element method[END_REF]. However, η r >1 are investigated in this preliminary section). Macroscopic shear strength in terms of peak, ϕ p , and residual, ϕ 0 , friction angles is shown in Figure 2. Note that the residual Fig. 1: The influence of the reduced contact stiffness κ on the peak friction angle (case 1 in Table 1).

friction angle is calculated from the average of stresses corresponding to ǫ a >25%. The results demonstrate that for a rolling stiffness, α r , sufficiently high (here for α r ≥ 1.25), the macroscopic plastic parameters reach constant values which depend only on η r [START_REF] Wang | Simulation of triaxial response of granular materials by modified DEM[END_REF].

Thus, if the contacts are sufficiently stiff, then the microelastic parameters, including rolling stiffness, have no influence on macro-plastic parameters. Approaching these normal and rolling stiff limits, should always be preferred to avoid cross dependencies between macroscopic plastic properties and elastic local parameters, in contrast with what was suggested in [START_REF] Plassiard | A spherical discrete element model: calibration procedure and incremental response[END_REF][START_REF] Widuliński | Numerical simulations of triaxial test with sand using DEM[END_REF].

Plastic local parameters and plastic macroscopic properties

The influence of the plastic local parameters, η r and ϕ c , on the plastic macroscopic properties can be analyzed independently from the local elastic stiffnesses (as

(a) (b)
Fig. 2: The variation of (a) peak friction angle, ϕ p , and (b) residual friction angle, ϕ 0 , with respect to α r for different values of η r (case 2 in Table 1).

long as the elastic stiffnesses are beyond the stiff limits). Therefore, E c is set in the following to 3 10 5 kPa (κ = 1875) and α r to 1.25. In addition, only cases for η r ≤ 1 are considered. Figure 3 shows the variation of the stress ratio ( q p ) and porosity (n= Vv Vt , where V v is the volume of voids and V t is the total volume of the sample) versus axial deformation during a triaxial compression for different η r values. As can be seen, both peak and residual shear strengths depend on η r , as well as on the dilatancy. The contact resistance to rolling increases with η r , mimicking more and more angular particles, leading to higher peak and residual friction angles and to more dilatant behaviors. This is also shown in Fig-ure 4 where the peak shear strength and the residual shear strength depend only on rolling friction, η r , for sufficiently high values of α r (i.e. without the effect of an elastic parameter).

On the other hand, a series of triaxial compressions were simulated with different values of the contact friction angle (Table 1, case 3). The response of the model is illustrated in Figure 5. It is shown that this parameter has a strong effect on the peak strength and dilatancy. However, for the initial porosity considered, only dilatant behaviors are observed, and purely contractant responses are not reproduced even for low values of ϕ c . Finally, for ϕ c > 15 • , its influence on the residual shear strength becomes insignificant, as observed similarly by [START_REF] Huang | Exploring the influence of interparticle friction on critical state behavior using DEM[END_REF]. Thus, for ϕ c > 15 • , the residual shear strength is governed only by η r , whereas the peak strength can be fine tuned by ϕ c for a fixed value of η r . Volumetric strain cannot be adjusted with ϕ c and/or η r alone since, for instance, it was not possible to reach a contractant behavior during the parametric studies displayed in Figures 3 and5. Nevertheless, the volumetric strain is also strongly dependent on the initial density, and a purely contractant behavior may be observed for loose or very loose granular materials. Consequently, a numerical preparation methodology allowing to reach a wide range of the initial porosity of the model is presented in the next section.

Preparation methodology for a large range of initial density

In the literature, several techniques to generate numerical discrete samples are suggested. They mainly include the isotropic-compression method [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] and the particle expansion method [START_REF] Plassiard | A spherical discrete element model: calibration procedure and incremental response[END_REF][START_REF] Wang | Simulation of triaxial response of granular materials by modified DEM[END_REF]. Both methods can be combined with a progressive lubrication of contacts consisting in decreasing the inter-particle friction until reaching the target porosity (as proposed for instance in [START_REF] Chareyre | Theoretical versus experimental modeling of the anchorage capacity of geotextiles in trenches[END_REF][START_REF] Tong | Pore-scale flow simulations: Model predictions compared with experiments on bi-dispersed granular assemblies[END_REF]). Such methods are efficient means to produce dense samples but they fail in repoducing very loose ones. Alternatively, the multi-layer with under-compaction method proposed in [START_REF] Jiang | An efficient technique for generating homogeneous specimens for DEM studies[END_REF] is capable of generating homogeneous samples with a variety of density conditions ranging from very loose to dense states. It is based on dividing evenly the sample into several layers, each compacted to a state looser that the upper one. However, little studies as the latter have been done on loose samples.

We propose thereafter a new method for that purpose. 1 for α r =1.25).

It imitates numerically the experimental reconstitution technique of soil samples by the moist tamping method

where a low moist content (about 2 to 3 %) introduces adhesion between particles allowing to reach even higher void ratios from the ones corresponding to the loosest state achieved by [START_REF] Jiang | An efficient technique for generating homogeneous specimens for DEM studies[END_REF]. Moreover, as explained in [START_REF] Gilabert | Computer simulation of model cohesive powders: Influence of assembling procedure and contact laws on low consolidation states[END_REF], introducing cohesion in the model, in the presence of rolling resistance, stabilizes looser and less coordinated samples leading to a reduction in the solid fraction and coordination number.

The methodology comprises three steps as follows: Fig. 4: The variation of peak and residual friction angles with respect to η r for different values of α r (case 2 in Table 1 for α r =1.25 and 2.5).

Random generation of a very loose cloud of nonoverlapping particles. 2. Isotropic compaction until a target confining pres-

sure is performed with the presence of an adhesive contact stress, σ A (Eq. 6) and keeping the normal contact friction angle (i.e. the same value used later for compression tests), first at a high strain rate to save computational time, then slower when contacts start to percolate, until equilibrium is achieved (see Figure 6). 3. Removal of the contact adhesion (σ A = 0) and waiting for a new equilibrium to be attained.

Figure 6 shows the evolution of the porosity and the unbalanced forces in each step. It can be seen that during the compaction process that the porosity decreases until an equilibrium state is reached. Equilibrium is assessed from the dimensenionless unbalanced force, U f . It is computed as the ratio of the mean resultant particle force to the mean contact force. U f tends to zero for a perfect static equilibrium and a value of U f = 10 -3 was considered here as representative of a state sufficiently close to this limit (Figure 6). Notice that during the last step where contact adhesion is removed, the sample experiences a slight additional compaction (see Figure 6), as observed experimentally with the moist tamping technique during the saturation stage of the soil sample.

Following the proposed procedure, the accessible densities reached by this method are shown in Figure 7. It can be noticed that as the contact adhesion, σ A , in- 1). creases, the initial porosity after compaction increases too, and the sample shows more and more an obvious contractive behavior as shown in Figure 8. This happens up to a threshold adhesion value (here equal to σ A = 80 kPa or σ A Pc =1), beyond which the numerical sample strongly collapses during the third step of the compaction method leading to lower porosities and more dilatant behavior (see the results for σ A Pc =1.25 and σ A Pc =2.5 in Figures 7 and8). Figure 8 shows the triaxial compressions simulated from the accessible densities reached by this method. It can be observed that removing completely the adhesion during the compaction phase is not enough to reach dense samples presenting a strong dilatant behavior. For that, contacts need to Fig. 6: Evolution of the porosity and the dimensionless unbalanced force during the three steps of compaction.

be lubricated during the compaction phase by fixing a low friction angle, as described for instance in [START_REF] Tong | Pore-scale flow simulations: Model predictions compared with experiments on bi-dispersed granular assemblies[END_REF], and illustrated in Figures 7 and8 for a lower values of the contact friction angle during compaction. Hence, a large range of initial density can be reached by assigning during compaction either an artificial contact adhesion or a reduced contact friction.

Modelling of the soil mechanical behaviour

In this section, a calibration procedure is defined and then applied to reproduce the experimental data obtained from drained triaxial compression tests with P c = 200 kPa (i.e. constant effective confining stress). Moreover, to validate our model, it was tested with the data obtained from: triaxial compression tests with different confining pressures, undrained triaxial test (i.e. no volume change) as well as constant deviator stress test (constant q) [START_REF] Nguyen | Etude expérimentale du comportement instable d'un sable silteux: Application aux digues de protection[END_REF].

Calibration process

After studying the influence of different parameters incorporated in our DEM model, a calibration procedure based on the simulated response to a triaxial compression was defined as follows:

1. As the residual shear strength is independent of the initial porosity and of the contact friction for ϕ c > 15 as a first approximation of the contact friction between two silica particles). For this step, the numerical sample may be compacted without the addition of contact adhesion or contact lubrication in order to start from a medium dense material. 2. Volumetric strain (i.e. dilatancy / contractancy) and peak stress (if any) are approached as close as possible by fine tuning the initial sample density. To get a looser material, adhesion is introduced during the compaction phase, while contacts are lubricated during the latter phase to reach a denser assembly, as described in the previous section. 1).

Contact friction angle ϕ c is calibrated (if necessary)

to improve the reproduction of the volumetric strain and the peak stress. Knowing that a change of ϕ c can affect the initial porosity reached after the compaction phase, it may be necessary to reiterate steps 2 and 3 until a satisfying calibration of both the initial porosity and ϕ c is achieved.

Calibration of elastic parameters, contact elastic modulus, E c , and shear factor, α s , appears secondary to plastic parameters (as long as they are beyond the stiff limit), with respect to their effects on the macroscopic mechanical behavior. They can be adjusted if necessary to better reproduce the macroscopic initial Young modulus and Poisson ratio as discussed in [START_REF] Plassiard | A spherical discrete element model: calibration procedure and incremental response[END_REF]. Rolling stiffness, α r , is just kept close to the rigidity limit, not lower to avoid any impact on plastic macroscopic properties, but not higher to avoid possible reduction of the critical time step as discussed in section 2.2.

It is worth noting that this proposed calibration procedure stresses more on the rolling friction than on the contact friction angle. Indeed, one can expect the contact friction angle not to be very different for different granular soils since grains are generally made with similar constituents (silica and alumine), whereas particle shapes, reflected by the rolling friction, may change a lot and directly impact the soil mechanical behavior.

Validation on laboratory triaxial compression paths

The calibration procedure is applied to calibrate the numerical model on experimental data obtained from drained triaxial compression tests on a Camargue sandy soil. Laboratory triaxial compression results are presented in [START_REF] Nguyen | Etude expérimentale du comportement instable d'un sable silteux: Application aux digues de protection[END_REF][START_REF] Benahmed | An experimental investigation into the effects of low plastic fines content on the behaviour of sand/silt mixtures[END_REF]. The triaxial tests were performed on very loose Camargue sand samples prepared with an initial (i.e. before isotropic compression) relative density of 20%. The particle size distribution of the numerical model partially follows the one of the sand as presented in Figure 9. It has been simplified by removing 3% of the largest particles and 3% of the smallest ones to limit the total number of discrete elements and so to reduce the computational cost. When applied to compressions with different confining pressures than the one used for calibration (Figure 10), the model succeeded to describe the shear strength but partially failed in reproducing quantitatively the dependence of the volumetric strain on the confining pressure. The important reduction of contractancy of the sand obtained experimentally for P c = 100 kPa, with respect to P c = 200 kPa and 400 kPa, is underestimated by the model.

Responses of the calibrated model and of the sand subjected to an undrained compression, after an initial isotropic compression, P c = 200 kPa, are compared in Figure 11. The model exhibits sudden and total liquefaction at very low axial strain which was not the case experimentally for the sand, as if the ability of the numerical model to contract was more important than for the sandy sample. Responses were also compared to a constant deviator stress loading path performed after an initial drained compression. Stress and strain paths are displayed in Figure 12. The numerical model correctly describes the stress state at failure occurrence but predicts a much more contractant response than the experimental data. The constant shear drained path, such as the q-constant path, lead to a reduction of the mean pressure p. Consequently, the overestimation of contractancy along these loading paths is in agreement with the results obtained from the drained compressions (Figure 10). In the latter case, we showed that the model underestimates the reduction of contractancy with the decrease of the mean pressure.

The observed difference of the dependency of volumetric strain on mean pressure can be related to different factors that differ between experimental and numerical tests:

1. Difference of the compressibility behavior induced by the isotropic compression: for the numerical sample, the relative variations of porosity after compaction when confining pressure is increased from Pc = 100 kPa to 200 kPa, and up to Pc = 400 kPa, are -0.09 % and -0.81 % respectively, whereas experimentally the sand sample presents a more important compressibility with relative porosity variations equal to -0.82 % and -2.22 % (always with respect to the state for Pc=100 kPa). This difference is particularly important in the range of confining pressure 100 kPa-200 kPa where the relative reduction of porosity for the numerical model is about 9 times smaller than what is found experimentally. This may contribute to the discrepancy between numerical and experimental volumetric strains ob-served, especially for the compression at Pc = 100 kPa. 2. Different particle shapes: this is related to the ability for angular or elongated particles to untangle. A reduction of the mean pressure, allowing more flexibility for angular or elongated particles to rotate, can promote their disentanglement, and thus a lower contractive behavior. In addition, even without any translation motion of particles, new contacts may be formed with elongated particles in rotation, limiting further compaction of the granular assembly. Such mechanisms are not reproduced with the numerical model made of spheres, even with the introduction of a contact resistance to rolling, and apparently, constitutes a limitation of the model.

Conclusion

This paper presented the DEM modeling of the soil mechanical behaviors by using spherical particles and contact rolling resistance. It was demonstrated that the macroscopic behavior of the granular assembly mainly depends, on one hand, on the initial packing density and, on the other hand, on the model parameter fixing the plastic rolling moment at contacts. Other parameters as normal contact stiffness or contact friction angle prove to play a secondary role (bearing in mind that the contacts should be sufficiently stiff and based on a given particle gradation). This is consistent with the fact that the shear strength in cohesionless soil is more strongly related to particle shape (presented here by the plastic rolling moment) than to the inter-particle friction angle, as shown previously by [START_REF] Estrada | Comparison of the effects of rolling resistance and angularity in sheared granular media[END_REF]. Moreover, it is greatly affected, together with dilatancy, by the initial density. Consequently, a new numerical method of compaction was proposed, involving either the introduction of interparticle adhesion forces or the lubrication of contacts, to reach a wide range of initial density of granular assemblies, from very loose to very dense. However, it was shown that even if our model succeeds in reproducing the experimental data, there are some limitations to be taken into account related to the description of the dependency of volumetric strains on the mean pressure. This is due to the difference in the relative variation of new force and torque induced by u and θ can be expressed as follows, in which the translational and the rotational degrees of freedom are uncoupled (this uncoupling assumption is also found in [START_REF]PFC3D theory and background manual[END_REF]): (Note that since F = T = 0, then we can write ∆ F = F and ∆ T = T ):

F ≃ - Nc c=1 (Kn -Ks)( n ⊗ n) + KsI u (12)
T ≃ -r 2

Nc c=1

Ks(I -n ⊗ n) θ

I being the inertia tensor of the particle and r its radius. These forces are given by neglecting the rotational term in equation 12 and the translational term in equation 13.

If contact moments ( Mr and Mtw) are present, the moments (initially null) induced at each contact (remaining elastic) by u and θ are:

Mtw = -Ktw( θ • n) n = -Ktw( n ⊗ n) θ, (14)
Mr = -Kr( θ -( θ • n) n) = -Kr(I -n ⊗ n) θ. (15)
where Kr and Ktw are the rolling and twisting stiffnesses, respectively. Note that, in our study, only the rolling moment was taken into consideration. However, for generalization, we added the twisting moment here.

Adding these contributions (eq. 14 and 15) to the global rotational stiffness (eq. 13) yields The uncoupled equations of motion for the particle can be expressed as:

ü = - Ku M u, (17) θ
= - K θ J θ, (18)
where Ku and K θ are the tensors defined by equations 12 and 16 respectively, M is the mass of the particle, and J is the moment of inertia (equal to 2 5 r 2 M for a sphere).

Ku

M and K θ J are symmetric positive definite tensors from which six eigen values can be calculated, corresponding to six natural frequencies of the particle. The explicit centered finite difference scheme is stable for the above equations -for a single particle b -as long as ∆t is less than 2/ λ b max , λ b max being the maximum eigenvalue of the particle.

A multi-body system as a whole may have natural oscillation modes combining the motion of multiple particles, with frequencies higher than the maximum frequency found for individual particles. An accurate determination of the critical timestep would thus need to consider the stiffness matrix of the system as a whole where the motion of two particles in contact would be included. For N b particles it would lead to a matrix of size (6N b ) 2 from which eigenvalues would have to be extracted. It can be avoided by remarking that this large assembled matrix would have every contact stiffness appearing twice in each row (once in the diagonal term and once in the column whose index corresponds to the other particle forming the contact). Since the absolute row-sum of the components can be used as an upper bound of the maximum eigenvalue (Perron-Froebenius theorem [START_REF] Brauer | The greatest distance between two characteristic roots of a matrix[END_REF]) an upper bound is obtained simply by multiplying every stiffness by two when calculating the frequencies of one single particle. Equivalently, we can keep the stiffness unchanged and define the upper bound using the eigenvalues associated to individual particles as

λmax < 2 max b (λ b max ). ( 19
)
[13] suggested a similar inequality without proof. Finally, a sufficient condition of stability is:

∆t < 2 max b (λ b max )

Fig. 3 :

 3 Fig.3: The variation of (a) the stress ratio and (b) porosity during triaxial drained compressions under different values of η r (case 2 in Table1for α r =1.25).

Fig. 5 :

 5 Fig. 5: The variation of (a) the stress ratio and (b) porosity during triaxial compressions under different values of ϕ c (case 3 in Table1).

Fig. 7 :

 7 Fig. 7: Range of porosity reached after compaction for: (a) introduction of contact adhesion σ a during compaction (with (φ c ) comp =30 • ), and (b) reduction of contact friction angle, (φ c ) comp , during compaction.

Fig. 8 :

 8 Fig. 8: Effect of adhesion, σ A , during compaction, on (a) the stress-strain response to a triaxial compression and (b) porosity (case 4 in Table1).

Fig. 9 :Fig. 10 :

 910 Fig. 9: The gradation curve of the Camargue sand.

Fig. 11 :

 11 Fig. 11: Comparison between experimental and numerical responses under undrained triaxial compression.

Fig. 12 :

 12 Fig. 12: Comparison between experimental and numerical responses under constant deviator stress loading path just up to failure point.

{r 2

 2 Ks(I -n ⊗ n) + Ktw n ⊗ n + Kr(I -n ⊗ n)} θ[START_REF] Jiang | An efficient technique for generating homogeneous specimens for DEM studies[END_REF]

Table 2 :

 2 Summary of the parameters calibrated on the loose Camargue sand.

	Pc	Ec	αs	ϕc	αr	ηr	σ A paction only) (for com-
	(kPa)	(kPa)		(deg)			(kPa)
	200	2 10 5	0.2	25	7.5	0.22	200

Acknowledgements

The environment and support provided by the French research group GDR MeGe 3176 is gratefully acknowledged. Laboratory 3SR is part of the LabEx Tec 21 (Investissements d'Avenir -grant agreement n • ANR-11-LABX-0030)

Appendix

Let u, θ define small movements of a particle relative to an initially stable configuration in which the total force and torque on the particle are F = T = 0. If the particle is involved in Nc contacts and every contact remains elastic (worst case), the