Luis F Contreras-Samamé

Salvador Domínguez-Quijada
email: salvador.dominguezquijada@ls2n.fr

Olivier Kermorgant
email: olivier.kermorgant@ls2n.fr

Philippe Martinet
email: philippe.martinet@inria.fr

CoMapping: Efficient 3D-Map Sharing Methodology for Decentralized cases

Keywords:

CoMapping is a framework to efficient manage, share, and merge 3D map data between mobile robots. The main objective of this framework is to implement a Collaborative Mapping for outdoor environments. The framework structure is based on two stages. During the first one, the Pre-Local Mapping stage, each robot constructs a real time pre-local map of its environment using Laser Rangefinder data and low cost GPS information only in certain situations. Afterwards, the second one is the Local Mapping stage where the robots share their pre-local maps and merge them in a decentralized way in order to improve their new maps, renamed now as local maps. An experimental study for the case of decentralized cooperative 3D mapping is presented, where tests were conducted using three intelligent cars equipped with LiDAR and GPS receiver devices in urban outdoor scenarios. We also discuss the performance of the cooperative system in terms of map alignments.

I. INTRODUCTION

Mapping the environment can be complex since in certain situations, e.g. in large regions, it may require a group of robots that build the maps in a reasonable amount of time with regards to the expected accuracy [START_REF] Dinnissen | Map merging of multi-robot SLAM using reinforcement learning[END_REF]. A set of robots extends the capability of a single robot by merging measurements from group members, providing each robot with information beyond their individual sensors range. This leads to a better usage of resources and execution of tasks which are not feasible by a single robot. Multi-robot mapping is considered as a centralized approach when it requires all the data to be analysed and merged at a single computation unit. Otherwise, in a decentralized approach, each robot builds their local maps independent of one another and merge their maps upon rendezvous.

Fig. 1. Scheme of our CoMapping System considering a decentralized case Figure 1 depicts the scheme of work proposed in this article for a group of robots where it was assumed that ZOE robot have direct exchange of data (as pose, size and limits of maps) with FLUENCE and GOLFCAR. And by contrast, FLUENCE and GOLFCAR are in a scenario of non-direct communication, with limited access conditions to a same environment, and without any meeting point for map sharing between these mobile units.

Following this scenario, this paper presents the development and validation of a new Cooperative Mapping framework (CoMapping) where:

• In the first stage named "Pre-Local Mapping", each individual robot builds its map by processing range measurements from a 3D LiDAR moving in six degrees of freedom (6-DOF) and using low cost GPS data (GPS/GGA). • For the second stage named "Local Mapping", the robots send a certain part of their pre-local maps to the other robots based on our proposed Sharing algorithm. The registration process includes an intersecting technique of maps to accelerate processing. This work is addressed to outdoor environments applications denied of a continuous GPS service by using a decentralized approach. Our proposal has been tested and validated in an outdoor environment, with data acquired on the surroundings of the ECN ( École Centrale Nantes) campus.

II. RELATED WORKS

In a scenario of cooperative mapping, robots first operate independently to generate individual maps. Here the registration method plays a fundamental role. Many registration applications use LiDAR as a Rangefinder sensor for map building [START_REF] Kohlbrecher | A flexible and scalable SLAM system with full 3D motion estimation[END_REF]. However, a high Lidar scan rate can be harmful for this task, since it may create distortion in the map construction. For those cases, ICP approaches [START_REF] Pomerleau | Comparing ICP variants on real-world data sets -open-source library and experimental protocol[END_REF] can be applied to match different scans. Implementations, for 2 or 3-axis and geometric structures matches of a generated local point set, were presented in [START_REF] Zlot | Efficient Large-Scale 3D Mobile Mapping and Surface Reconstruction of an Underground Mine[END_REF] [START_REF] Bosse | Zebedee: Design of a Spring-Mounted 3-D Range Sensor with Application to Mobile Mapping[END_REF]. Those methods use batch processing to build maps offline. In the first stage of our implementation we reconstruct maps as 3D point clouds in real-time using 3axis LiDAR by extraction and matching of geometric features in Cartesian space based in [START_REF] Zhang | LOAM: Lidar odometry and mapping in realtime[END_REF]. Then our system uses GPS position data to coarsely re-localize that cloud in a global frame.

Once all the maps have been projected on a global frame, they have to be merged together to form a global map. In this context, [START_REF] Jessup | Merging of octree based 3D occupancy grid maps[END_REF] proposed a method for 3D merging of occupancy grid maps based on octrees [START_REF] Hornung | OctoMap: An Efficient Probabilistic 3D Mapping Framework Based on Octrees[END_REF] for multi-robots.

The method was validated in a simulation environment. For the merging step, an accurate transformation between maps was assumed as known, nevertheless in real applications, that information is not accurate, since it is obtained by means of uncertain sensor observations. In our case, real experiments were performed for a multi-robot application without perfect knowledge of the transformations between the maps. In [START_REF] Jessup | Robust and Efficient Multirobot 3-D Mapping Merging With Octree-Based Occupancy Grids[END_REF] a pre-merging technique is proposed, which consists in selecting the subset of points included in the common region between maps bounding. Then, a centralized merging process refines the transformation estimate between maps by ICP registration [START_REF] Pomerleau | Comparing ICP variants on real-world data sets -open-source library and experimental protocol[END_REF]. We use a variation of [START_REF] Jessup | Robust and Efficient Multirobot 3-D Mapping Merging With Octree-Based Occupancy Grids[END_REF] but previously we include an efficient technique to exchange maps between robots in order to optimize bandwidth.

Other solutions may be used in order to merge maps for a group of robots with a centralized approach [START_REF] Birk | Merging occupancy grid maps from multiple robots[END_REF], [START_REF] Jessup | Robust and Efficient Multirobot 3-D Mapping Merging With Octree-Based Occupancy Grids[END_REF], [START_REF] Dubé | An online multi-robot SLAM system for 3D LiDARs[END_REF], which are generally found in the literature. However this kind of solutions can compromise the team performance because merging and map construction depend exclusively on a processing Another approach less explored and analysed is the decentralized one [START_REF] Özkucur | Supervised feature type selection for topological mapping in indoor environments[END_REF], [START_REF] Dinnissen | Map merging of multi-robot SLAM using reinforcement learning[END_REF], [START_REF] Zhang | Aerial and Ground-based Collaborative Mapping: An Experimental Study[END_REF]. The advantage in this case lies in the independence and robustness of the system because map construction is not affected even if one of the robots has failures in communication or processing, since map merging can be executed in different units after traversing the environment. This kind of approach can consider a meeting point for the vehicles in order to exchange their maps and other data. This approach is also investigated in our final experiments.

III. METHODOLOGY

A. Pre-Local Mapping Stage

Each mobile robot executes a Pre-Local Mapping system using data provided by a LidarSLAM process. We just use coarse GPS position to project the generated map on a global frame. In order to reduce implementation costs, a beneficial cheap GPS service was used, specifically GPS/GGA (Global Positioning System Fix Data) at an accuracy of about 2 to 7 meters. Another advantage of our Pre-Local Mapping Stage is its versatile configuration, since it does not depend on a specific LidarSLAM method.

A modified version of the LOAM technique 1 [START_REF] Zhang | LOAM: Lidar odometry and mapping in realtime[END_REF] was thus chosen as the LidarSLAM method for this work. is denoted P k and processed by a Lidar Odometry algorithm, which runs at a frequency around 10Hz and computes the lidar motion (transform T k ) between two consecutive sweeps. The distortion in P k is corrected using the estimated lidar motion. The resulting undistorted P k is processed at a frequency of 1Hz by an algorithm knows as Lidar Mapping, which performs the matching and registration of the undistorted cloud onto a map. Using the GPS information of the vehicle pose, it is possible to coarsely project the map of each robot into common coordinate frame for all the robots. This projected cloud is denoted as the Pre-Local Map.

1) Lidar Odometry step:

The step begins with feature points extraction from the cloud P k . The feature points are selected from sharp edges and planar surface patches. Let us define S as the set of consecutive points i returned by the laser scanner in the same scan, where i ∈ P k . An indicator proposed in [START_REF] Zhang | LOAM: Lidar odometry and mapping in realtime[END_REF] evaluates the smoothness of the local surface as following:

c = 1 | S | . X L (k,i) j∈S,j =i (X L (k,i) -X L (k,j) ) , (1)
where X L (k,i) and X L (k,j) are the coordinate of two points from the set S.

Moreover, a scan is split into four subregions to uniformly distribute the selected feature points within the environment. In each subregion is determined maximally two edge points and four planar points. The criteria to select the feature points as edge points is related to maximum c values, and by contrast the planar points selection to minimum c values. When a point is selected, it is thus mandatory that none of its surrounding points are already selected. Besides, selected points on a surface patch cannot be approximately parallel to the laser beam, or on boundary of an occluded region.

When the correspondences of the feature points are found, then the distances from a feature point to its correspondence are calculated. Those distances are named as d E and d H for edge points and planar points respectively. The minimization of the overall distances of the feature points leads to the Lidar odometry. That motion estimation is modeled with constant angular and linear velocities during a sweep.

Let us define E k+1 and H k+1 as the sets of edge points and planar points extracted from P k+1 , for a sweep k+1. The lidar motion relies on establishing a geometric relationship between an edge point in E k+1 and the corresponding edge line:

f E (X L (k+1,i) , T L k+1 ) = d E , i ∈ E k+1 , (2)
where T L k+1 is the lidar pose transform between the starting time of sweep k + 1 and the current time t i . Analogously, the relationship between a planar point in H k+1 and the corresponding planar patch is:

f H (X L (k+1,i) , T L k+1 ) = d H , i ∈ H k+1 , (3)
Equations ( 2) and ( 3) can be reduced to a general case for each feature point in E k+1 and H k+1 , leading to a nonlinear function:

f (T L k+1 ) = d, (4)
in which each row of f is related to a feature point, and d possesses the corresponding distances. Levenberg-Marquardt method [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] is used to solve the Equation (4). Jacobian matrix (J) of f with respect to T L k+1 is computed. Then, the minimization of d through nonlinear iterations allows solving the sensor motion estimation:

T L k+1 ←-T L k+1 -(J T J + λdiag(J T J)) -1 J T d, ( 5
)
where λ is the Levenberg-Marquardt gain.

Finally, the Lidar Odometry algorithm produces a pose transform T L k+1 that contains the lidar tracking during the sweep between [t k+1 , t k+2 ] and simultaneously an undistorted point cloud Pk+1 . Both outputs will be used by the Lidar Mapping step, detailed in the next section.

2) Lidar Mapping step: This algorithm is used only once per sweep and runs at a lower frequency (1 Hz) than the Lidar Odometry step (10 Hz). The technique matches, registers and projects the cloud Pk+1 provided by the Lidar Odometry as a map into the coordinate system of a vehicle, noted as {V }. To understand the technique behaviour, let us defined Q k as point cloud accumulated until sweep k, and T V k as the sensor pose on the map at the end of sweep k, t k+1 . The algorithm extends T V k for one sweep from t k+1 to t k+2 , to get T V k+1 , and projects Pk+1 on the robot coordinate system, denoted as Qk+1 . Then, by optimizing the lidar pose T V k+1 , the matching of Qk+1 with Q k is obtained.

In this step the feature points extraction and their correspondences are calculated in the same way as in Lidar Odometry, the difference just lies in that all points in Qk+1 share the time stamp, t k+2 .

In that context, nonlinear optimization is solved also by the Levenberg-Marquardt method [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], registering Qk+1 on a new accumulated cloud map. To get a uniform points distribution, down-sampling is performed to the cloud using a voxel grid filter [START_REF] Rusu | 3D is here: Point Cloud Library (PCL)[END_REF] with a voxel size of 5 cm cubes.

Finally, since we have to work with multiple robots, we use a common coordinate system for their maps, {W }, coming from rough GPS position estimation of the 1st accumulated cloud frame Q k .

B. Local Mapping Stage

In this section the Local Mapping is detailed, considering that the process is executed on the robot "i" with a shared map by robot "n" (see Figure 3).

1) Map Sharing

Step: When the generation of Pre-Local Maps is done, the robots would have to exchange their maps to start the map alignment process. In several cases the sharing and processing of large maps can affect negatively the performance of the system with respect to runtime and

Map Intersecting Algorithm

Map Sharing Step

Robot "n" memory usage. A sharing technique is presented in order to overcome this problem, in which each vehicle only sends a certain subset of points of its map to the other robots. When the maps are ready for transferring, they are compressed in octree format using OctoMap library [START_REF] Hornung | OctoMap: An Efficient Probabilistic 3D Mapping Framework Based on Octrees[END_REF] in order to optimize the robot-communication.

Local Mapping

Map limits A

The proposed sharing technique is based on the method developed in [START_REF] Contreras | Efficient Decentralized Collaborative Mapping for Outdoor Environments[END_REF]. Figure 4 depicts the behaviour, wherein point clouds A and B represent the Pre-Local Maps from two robots "i" and "n" respectively. In each robot the algorithm first receives only information about the 3D limits of the maps (i.e. bounding cubic lattice of the point clouds) and then decides what part of its map will be shared to the other robot. These limits were determined previously using the function GetBounds() that returns two vectors: in the first one, Amin, their components represent the lowest displacement from the origin along each axis in the point cloud; and the other vector, Amax, is related to the point of the highest displacement. Pseudo-code of the map sharing step is described in Algorithm 1. Inside the code, the function GetV alues() sorts in ascending order the array of components along each axis of the vectors Amin, Amax, Bmin, Bmax and returns the 2nd and 3rd values from this sorted array, denoted V 2 and V 3 respectively. Then for each axis, the average of the two values obtained by the function GetV alues() is used in order to determine the Cartesian coordinate (C x ,C y ,C z ) of the geometric center of the sharing region S. This map sharing region is a cube whose edge length 2L is determined iteratively. The Points from A contained in this cube region are extracted to generate a new point cloud A sel . In each iteration the cube region is reduced until the number of points from A sel is smaller than the manual parameter N p max , which represents the number of points maximum that the user wants to exchange between robots. Once the loop ends, A sel is sent to the other robot. Analogously on the other mobile robot "n", the points from B included in this region are also extracted to obtain and share B sel with the another robot "i". Finally, the clouds A sel and B sel are encoded and sent in octree format to reduce the usage of bandwidth resources of the multi-robot network. Then maps are decoded and reconverted in 3D point cloud format to be used in the Registration step.

2) Registration

Step: The intersecting volumes of the two maps A sel and B sel are computed and denoted as A int and B int , obtained from the exchanged map bounds [START_REF] Jessup | Robust and Efficient Multirobot 3-D Mapping Merging With Octree-Based Occupancy Grids[END_REF]. In order to improve the computation speed, point clouds A int to B int first go through a down-sampling process to reduce the number of points in the cloud alignment. Feature descriptors as surface normals and curvature are used to improve the matching, which is the most expensive stage of the registration algorithm [START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF]. These generated normal-point clouds A intN and B intN are then used by Iterative Closest Point (ICP) algorithm [START_REF] Besl | A Method for Registration of 3-D Shapes[END_REF]. This method refines an initial alignment between clouds, which basically consists in estimating the best transformation to align a source cloud B intN to a target cloud A intN by iterative minimization of an error metric function. At each iteration, the algorithm determines the corresponding pairs (b', a'), which are the points from A intN and B intN respectively, with the least Euclidean distance.

Then, least squares registration is computed and the mean squared distance E is minimized with regards to estimated translation t and rotation R:

E(R, t) = 1 N p b' N p b' i=1 a' i -(R b' i + t) 2 , (6)
where N p b is the number of points b'.

The resulting rotation matrix and translation vector can be express in a homogeneous coordinate representation (4×4 transformation matrix T j ) and are applied to B intN . The algorithm then re-computes matches between points from A intN and B intN , until the variation of mean square error between iterations is less than a defined threshold. The final ICP refinement for n iterations can be obtained by multiplying the individual transformations: T ICP = n j=1 T j . Finally the transformation T ICP is applied to the point cloud B sel to align and merge with the original point cloud A, generating the Local Map A L then. Each robot thus performed its own merging according to data from other agents within communication range. We now present the corresponding experimental results. In this section we show results validating the presented concepts and the functionality of our system. As we consider ground vehicles, the ENU (East-North-Up) coordinate system is used as external reference of the world frame {W }, where y-axis corresponds to North and x-axis corresponds to East, The proposed framework is validated considering three vehicles for experiments, a ZOE Renault, a FLUENCE Renault and a GOLFCAR (see Figure 5) customized and equipped with a Velodyne VLP-16 3D LiDAR, with 360 • horizontal and a 30 • vertical field of view. All data come from the campus outdoor environment in an area of approximately 1000m x 700m. The vehicles traversed that environment following different paths and collected sensor observations about the world, running pre-local mapping process in real-time.

IV. RESULTS

For the validation, the vehicles build clouds from different paths (see Figure 6). Results of the Pre-Local Mapping of this experiment are shown in Figure 7. As we assume a decentralized scenario, each robot performs a relative registration process considering its Pre-Local map as target cloud for alignment reference. Each vehicle also executes the intersecting algorithm and then an ICP refinement to obtain an improved transform between each map. Figures 8 and9 depict in yellow color the intersection between the shared point clouds during the alignment process. Once the refined transformation is obtained, it is then applied to the shared map. I corresponds to the merging process in ZOE, when this robot received the map shared by FLUENCE and it aligned that map to its own pre-local map. The decentralized system demonstrated alignments in opposite directions for both robots, since each robot performs the merging process considering its Pre-Local map as target cloud for alignment reference. Table II also reveals this symmetrical behavior, where the algorithm on ZOE converged to the value of displacement of -0.1782 m and -3.2605 m along the x-axis and y-axis respectively. On the other hand on the GOLFCAR robot, the algorithm converged to a value of displacement of 0.2213 m and 3.3857 m along the x-axis and y-axis respectively, reconfirming relative alignments in opposite directions.

Figure 10 shows one of the merging results corresponding to the ZOE robot, in which the cloud represents the final 3D local map projected on a 2D map in order to make qualitative comparisons. Experiments showed the impact of working with intersecting regions, since it can accelerate the alignment process by decreasing the number of points to compute. In the same way, tests demonstrated that our proposed map sharing technique developed a transcendental position in the performance of the entire mapping collaborative system by reducing the map size to transmit. Finally, the sharing algorithm proves to be a suitable candidate to exchange efficiently maps between robots considering the use of clouds of large dimensions.

V. CONCLUSION AND FUTURE WORK

A framework was presented for decentralized 3D mapping system for multiple robots. The work has showed that maps from different robots can be successfully merged, from a coarse initial registration and a suitable exchange of data volume. The system uses initially range measurements from a 3D LiDAR, generating a pre-local maps for each robot. The complete system solves the mapping problem in an efficient and versatile way that can run in computers dedicated to three vehicles for experiments, leading to merge maps independently on each vehicle for partially GPS-denied environments. Future work will focus on the analysis of the consistency of the final maps estimated on each robots

Fig. 2 .

 2 Fig. 2. Architecture of Pre-Local Mapping Stage

Fig. 4 .

 4 Fig. 4. Graphical representation of the Map Sharing technique (Top view of plane XY). Aminx, Amaxx, Bminx and Bmaxx represent the point cloud limits along the x-axis.

Algorithm 1 :

 1 Selection of Point Cloud to share with another robot.

Fig. 5 .

 5 Fig. 5. Vehicles used in the tests: ZOE, FLUENCE and GOLFCAR.

Fig. 6 .

 6 Fig. 6. Paths followed by ZOE (green one), FLUENCE (red one) and GOLFCAR robot (blue one) during experiments. Image source: Google Earth.

Fig. 7 .

 7 Fig. 7. Top view of unaligned Pre-Local Maps generated by ZOE (green one), FLUENCE (red one) and GOLFCAR robot (blue one) projected on common coordinate system

Figure 7

 7 Figure 7 also depicts the "sharing region" determined during the map exchange process in each robot. It was assumed that all the vehicles have the constraint of exchanging a maximum number of points N p max of 410000 to simulate restrictions in resources of bandwidth network or memory usage in robots. The tests were divided in two. In the first one, test A, ZOE and FLUENCE car define a meeting point to transfer their maps. Around this position, ZOE car exchanges and updates its local map and then a new point of rendezvous for map sharing is determined by ZOE and GOLFCAR in the following test B.

Fig. 8 .

 8 Fig. 8. Test A: Alignment of the intersecting regions with ICP refinement performed in ZOE robot, when it received the FLUENCE map (a) Green and red maps represent the target and source clouds pre ICP, top view (b) Green and blue maps represent the target and aligned source clouds post ICP, top view. Alignment can be better appreciated in yellow box

Fig. 10 .

 10 Fig. 10. Top view of final 3D Local Map of ZOE robot (color). Uniform green point clouds come from the aligned maps of FLUENCE (left) and GOLFCAR (right)

TABLE I TEST

 I A: RELATIVE ICP TRANSFORMATIONS IN EULER FORMAT BETWEEN ZOE AND FLUENCE ROBOT

	Robot	x	y	z	roll	pitch	yaw
	ZOE	-1.6517	3.0966	-9.9729	0.0132	0.0730	0.0022
	FLU.	4.5748	-4.4556	6.6061	-0.0054	-0.0624	-0.0084

TABLE II TEST

 II B: RELATIVE ICP TRANSFORMATIONS IN EULER FORMAT BETWEEN ZOE AND GOLFCAR ROBOT

	Robot	x	y	z	roll	pitch	yaw
	ZOE	-0.1782	-3.2605	1.7771	-0.0516	0.0115	0.0356
	GOL.	0.2213	3.3857	-2.6070	0.0411	-0.0256	-0.0380

ACKNOWLEDGMENT This article is based upon work supported by the Erasmus Mundus Action 2 programme through the Sustain-T Project. The authors also gratefully acknowledge all the members of ECN, LS2N and INRIA institution on this work. Parts of the equipments used here were funded by the project ROBOTEX, reference ANR-10-EQPX-44-01.


Mariano Tepper
email: mariano.tepper@duke.edu

Pablo Musé

Andrés Almansa
email: andres.almansa@telecom-paristech.fr

Marta Mejail

Finding Contrasted and Regular Edges by a Contrario Detection of Periodic Subsequences

Keywords: topographic maps, level lines, periodic binary sequences, edge detection, Helmholtz principle

A new method to detect salient pieces of boundaries in an image is presented. After detecting perceptually meaningful level lines, periodic binary sequences are built by labeling each point in close curves as salient or non-salient. We propose a general and automatic method to detect meaningful subsequences within these binary sequences. Experimental results show its good performance, when tested with different saliency criteria, such as contrast, regularity, and the combination of both.

Introduction

Shape plays a key role in our cognitive system: in the perception of shape lies the beginning of concept formation. Formally, shapes in an image can be defined by extracting contours from solid objects. Shapes can be represented and analyzed as the locus of an infinite number of points, which leads to levelsets methods [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF].

We define an image as a lower (or upper) semi-continuous function u : R 2 → R. Level sets [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF], or level lines, provide a complete and contrast-invariant image description. We define the boundaries of the connected components of a level set as a level line. These level lines have the following properties: (1) level lines are closed Jordan curves; (2) level lines at different levels are disjoint; (3) by topological inclusion, level lines form a partially ordered set.

We call the collection of level lines (along with their level) a topographic map. The inclusion relation allows to embed the topographic map in a treelike representation. For extracting the level lines of a digital image we use the Fast Level Set Transform (FLST) [START_REF] Monasse | Fast Computation of a Contrast Invariant Image Representation[END_REF] which computes level lines by bilinear interpolation. In general, the topographic map is an infinite set and so only quantized grey levels are considered, ensuring that the set is finite.

Edge detectors, from which the most renowned is Canny's [START_REF] Canny | A Computational Approach to Edge Detection[END_REF], rely on the fact that information is concentrated along contours (regions where contrast changes abruptly). From one side, only a subset of the topographic map is necessary to obtain a perceptually complete description. Going to a deeper level, perceptually important level lines, in general, are so because they contain contrasted pieces. In summary, we have to prune the topographic map and then prune inside the level lines themselves.

The search for perceptually important contours will focus on unexpected configurations, rising from the perceptual laws of Gestalt Theory [START_REF] Kanizsa | Organization in Vision: Essays on Gestalt Perception[END_REF]. From an algorithmic point of view, the main problem with the Gestalt rules is their qualitative nature. Desolneux et al. [START_REF] Desolneux | From Gestalt Theory to Image Analysis[END_REF] developed the Computational Gestalt detection theory which seeks to provide a quantitative assessment of gestalts. It is primarily based on the Helmholtz principle which states that conspicuous structures may be viewed as exceptions to randomness. In this approach, there is no need to characterize the elements one wishes to detect but contrarily, the elements one wishes to avoid detecting, i.e., the background model. When an element sufficiently deviates from the background model, it is considered meaningful and thus, detected.

Within this framework, Desolneux et al. [START_REF] Desolneux | From Gestalt Theory to Image Analysis[END_REF] proposed an algorithm to detect contrasted level lines in grey level images, called meaningful boundaries. Further improvements to this algorithm were proposed by Cao et al. [START_REF] Cao | A Theory of Shape Identification[END_REF] and by Tepper et al. [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF], which include for example, the use of regularity as the saliency measure instead of the contrast.

In this work we address the dissection of meaningful boundaries, developing an algorithm to select salient pieces contained in them. Each level line is considered as a periodic binary sequence where, following a partial saliency model, each point is labeled as salient or non-salient. Then, the goal is to extract meaningful subsequences of salient points. To do so, in this work we extend to the periodic case an algorithm for detecting binary subsequences.

Grompone et al. [START_REF] Grompone | On Straight Line Segment Detection[END_REF] proposed a method for accurately detecting straight line segments in a digital image. It is based on the Helmholtz principle and hence it is parameterless. In the authors' words, "at the core of the work lies a new way to interpret binary sequences in terms of unions of segments." In the same spirit, pieces of level lines correspond to object contours and can be recovered extending this approach to the periodic case (level lines are closed curves).

The remainder of this paper is organized as follows. We begin, in Section 2, by explaining the proposed algorithm for detecting periodic subsequences in binary sequences. Then, in Section 3 we recall several variants of the meaningful boundaries algorithm [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF]. In Section 4 we show how to apply the proposed algorithm to the problem of selecting salient pieces of level lines. In Section 5 we discuss the pertinence of the approach and provide some final remarks.

Detecting Periodic Subsequences

A sequence S = (s i ) 1≤i≤L of length L is binary if ∀i, s i ∈ {0, 1}. A subsequence a ⊆ S is defined by a pair of indices a (1) , a (2) with 1 ≤ a (1) < a (2) ≤ L, such that ∀s i , a (1) ≤ i ≤ a (2) s i ∈ a. Given a binary sequence S of length L, an n-subsequence is an n-tuple (a 1 , . . . , a n ) of n disjoints subsequences a i ⊆ S.

The set of all n-subsequences in S will be denoted by M(n, S). We define k(a) = #{s i | i ∈ [a (1) , a (2) ] ∧ s i = 1} and l(a) = a (2) -a (1) + 1 (i.e. the length of a). Notice that #M(n, S) = L 2n [START_REF] Grompone | On Straight Line Segment Detection[END_REF]. Definition 1. (Grompone et al. [START_REF] Grompone | On Straight Line Segment Detection[END_REF]) Given a binary sequence S of length L, an n-subsequence

(a 1 , . . . , a n ) in M(n, S) is said ε-meaningful if NFA(a 1 , . . . , a n ) = L 2n n i=1 (l(a i ) + 1)B(l(a i ), k(a i ), p) < ε (1)
where p = Pr(s i = 1), 1 ≤ i ≤ L. This number is called number of false alarms (NFA) of (a 1 , . . . , a n ).

Proposition 1. The expected number of ε-meaningful n-subsequences in a random binary sequence is smaller than ε.

We refer to the work by Grompone et al. [START_REF] Grompone | On Straight Line Segment Detection[END_REF] for a complete proof.

A run in S is a maximal subsequence only containing ones, i.e.

∀i ∈ [a (1) , a (2) ],

s i = 1 ∧ a (1) = 1 ∨ s a (1) -1 = 0 ∧ a (2) = L ∨ s a (2) +1 = 0 .
One can restrict the search for n-subsequences to the ones where each of the n subsequences starts at a run start and ends at a run end [START_REF] Grompone | On Straight Line Segment Detection[END_REF]. We denote by R the number of runs in S.

Definition 2. Given a binary sequence S, its maximal ε-meaningful subsequence (a 1 , . . . , a n ) * is defined as

(a 1 , . . . , a n ) * = arg min 1≤n≤R (a1,...,an)∈M(n,S) NFA(a 1 , . . . , a n ).
We propose now to extend the above definitions to support periodic binary sequences. A binary sequence S = (s i ) 1≤i≤L is made periodic by considering L its period. Periodic sequences are different in nature from their non-periodic counterparts, see Fig. 1. A definition suitable for the periodic case is needed.

In the periodic case, a subsequence must be defined more carefully. Now a subsequence a ⊆ S, defined by a pair of indices a (1) , a (2) , can belong to one of two different types: intra-subsequences: if a (1) < a (2) then the non-periodic definition holds, i.e., 1 ≤ a (1) < a (2) ≤ L, and ∀s i , a (1) ≤ s i ≤ a (2) s i ∈ a.

inter-subsequences: if a (1) > a (2)

∀s i , 1 ≤ s i ≤ a (2) ∨ a (1) ≤ s i ≤ L s i ∈ a.
Runs are modified accordingly to also cover inter-subsequences. Given a periodic binary sequence S of period L, a periodic n-subsequence is an n-tuple (a 1 , . . . , a n ) of n disjoints subsequences a i ⊆ S. The set of all n-subsequences in S will be denoted by M(n, S).

We define ) , a (2) ] ∧ s i = 1} and the length of a as

k(a) = #{s i | i ∈ [a (1
l(a) = a (2) -a (1) + 1, if a is an intra-subsequence; a (2) + L -a (1) + 1, if a is an inter-subsequence.
Notice that #M(n, S) = 2 L 2n since from each pair of points in S two subsequences can be constructed. Definition 3. Given a periodic binary sequence S of period L, an n-subsequence

(a 1 , . . . , a n ) in M(n, S) is said ε-meaningful if NFA(a 1 , . . . , a n ) = 2 L 2n n i=1 (l(a i ) + 1) B(l(a i ), k(a i ), p) < ε
where p = Pr(s i = 1), 1 ≤ i ≤ L. This number is called number of false alarms (NFA) of (a 1 , . . . , a n ).

Proposition 2. The expected number of ε-meaningful n-subsequences in a random periodic binary sequence is smaller than ε.

Proof. This proof follows closely the one by Grompone et al.

P (NFA(a 1 , . . . , a n ) < ε) . NFA(a 1 , . . . , a n ) < ε implies that n i=1 B(l(a i ), k(a i ), p) < ε 2 L 2n n i=1 (l(a i ) + 1) . Let U i = B(l(a i ), k(a i ), p) be a random variable, let α ∈ R + , and let P α U = P ( n i=1 U i < α). Then, P α U = u2,...,un P n i=1 U i < α U 2 = u 2 , . . . , U n = u n P (U 2 = u 2 , . . . , U n = u n ) .
Since the a i are disjoint, the U i are independent. Then

P α U = u2,...,un P n i=1 U i < α u 2 . . . u n • P (U 2 = u 2 , . . . , U n = u n ) .
Using the classical lemma P(

U i < α) < α, that P (U 2 = u 2 , . . . , U n = u n ) ≤ P (U 2 ≤ u 2 , . . . , U n ≤ u n ),
and that there are l(a i ) + 1 possible values for U i ,

P n i=1 U i < α < n i=2 (l(a i ) + 1) α < n i=1 (l(a i ) + 1) α.
Let us recall that #M(n, S) = 2 L 2n , then setting α = ε 2 ( L 2n ) n i=1 (l(ai)+1) gives the wanted result.

The maximality rule from Def. 2 holds unchanged in the periodic case. On the implementation side, Grompone et al. [START_REF] Grompone | On Straight Line Segment Detection[END_REF] describe a dynamic programming scheme for the non-periodic case that eases the heavy computational burden. We show now that implementing the algorithm for detecting periodic subsequences is indeed straightforward.

We begin by shifting the periodic sequence S (with R runs), to transform inter-subsequences into intra-subsequences. A circular shift to the left is used. We then form a non-periodic sequence S (2) of length 2L from two periods of the periodic sequence S of period L. Let R (2) be the number of runs in S (2) . Two key tricks allow us to solve the problem:

1. restrict the number of tested subsequences. In the non-periodic case, we test for n-subsequences for S (2) where 1 ≤ n ≤ R (2) . In the periodic case, we only test for n-subsequences where 1 ≤ n ≤ R; 2. subsequences longer than L are not tested.

With these two restrictions, one can simply detect non-periodic subsequences in non-periodic sequence S (2) and the result will be optimal.

Meaningful Boundaries

We briefly reproduce the formal definitions that lead to the meaningful boundaries algorithm, as defined by Tepper et al. [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF]. Several criteria are used to determine the saliency of a boundary: contrast, regularity and the combination of both. All these variants will be used to test the proposed algorithm.

Contrasted Meaningful Boundaries

Let C be a level line of the image u and let us denote by {x i } i=0...n-1 the n regularly sampled points of C, with arc-length two pixels, which in the a contrario noise model are assumed to be independent. In particular the gradients at these points are independent random variables (the image gradient norm |Du| can be computed by standard finite differences on a 2 × 2 neighborhood). We note by µ k (0 ≤ k < n) the k-th value of the values |Du|(x i ) sorted in ascending order.

The detection algorithm consists in rejecting the null hypothesis H 0 : the line C with contrasts {µ k } k=0...n-1 is observed only by chance. For this we assume that the values of |Du| are i.i.d., extracted from a noise image with the same gradient histogram as the image u itself.

Desolneux et al. [START_REF] Desolneux | From Gestalt Theory to Image Analysis[END_REF] present a thorough study of the binomial tail B(n, k; p) and its use in the detection of geometric structures. The regularized incomplete beta function, defined by I(x; a, b) is an interpolation B of the binomial tail to the continuous domain B(n, k; p) = I(p; k, n -k + 1) where n, k ∈ R [START_REF] Desolneux | From Gestalt Theory to Image Analysis[END_REF]. Additionally the regularized incomplete beta function can be computed very efficiently.

Let H c (µ) = P(|Du| > µ). For a given line of length l, the probability under H 0 that, some parts with total length greater or equal than l (s,n) (n -k) have a contrast greater than µ can be modeled by B(n • l (s,n) , k • l (s,n) ; H c (µ)), where l (s,n) = l s•n acts as a normalization factor [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF].

Definition 4. Let C be a finite set of N ll level lines of u. A level line C ∈ C is a ε-meaningful boundary if NFA K (C) = N ll • K • min k<K B(n • l (2,n) , k • l (2,n) ; H c (µ k )) < ε, (2)
where K is a parameter of the algorithm. We also note

k min = arg min k<K B(n • l (2,n) , k • l (2,n) ; H c (µ k )). (3)
The parameter K controls the number of points that we allow to be likely generated by noise, that is a line must have no more than K points with a "high" probability of belonging to the background model. It is simply chosen as a percentile of the total number of points in the line.

Def. 4 is motivated by the following proposition (we refer to the work by Tepper [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF] for a complete proof).

Proposition 3. The expected number of ε-meaningful boundaries, in a finite set of random level lines is smaller than ε.

x = C(0) C(-s) C(s) s × R s (x)
Figure 2: Reproduced from the work by Cao et al. [START_REF] Cao | Extracting Meaningful Curves from Images[END_REF]. The regularity at x is obtained by comparing the radius of the circle with s. The radius is equal to s if and only if the curve is a straight line. If the curve has a large curvature, the radius will be small compared to s.

Combining contrast and smoothness

As already stated, in natural images contrasted boundaries often locally coincide with object edges. Thus, they are also incidentally smooth. Active contours [START_REF] Kass | Snakes: Active contour models[END_REF] rely on this combination of good contrast and smoothness to provide well localized contours. In this section, we reproduce results that use smoothness as the saliency measure in the a contrario detection process [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF].

Let C be a rectifiable planar curve, parameterized by its length. Let l be the length of C and x = C(τ ) ∈ C. Without loss of generality, we assume that τ = 0. Definition 5. (Cao et al. [START_REF] Cao | Extracting Meaningful Curves from Images[END_REF]) Let s > 0 be a fixed positive value such that 2s < l. We call regularity of C at x (at scale s) the quantity

R s (x) = max(|x -C(-s)|, |x -C(s)|) s , (4)
where |x i -x j | represents the Euclidean distance between x i and x j .

Figure 2 visually explains the pertinence of this definition. Only when one of the subcurves C((-s, 0)) or C((0, s)) is a line segment, R s (x) = 1; in all other cases R s (x) < 1. When s is small enough, regularity is inversely proportional to the curve's curvature around x [START_REF] Cao | Extracting Meaningful Curves from Images[END_REF].

The question about the choice of s arises naturally and was studied in detail by Cao et al. [START_REF] Cao | Extracting Meaningful Curves from Images[END_REF] and Musé [START_REF] Musé | On the definition and recognition of planar shapes in digital images[END_REF]. We will limit ourselves to state that a larger value of s (thus at less local scale of analysis) is more robust to noise. On the other side, s should not be too large either. In practice, and following Cao et al. [START_REF] Cao | Extracting Meaningful Curves from Images[END_REF] one may safely set s = 5, which is the value we use in our experiments.

Let us denote by H s (r) the distribution of the regularity in white noise level lines, i.e., H s (r) = P R s (x) > r, x ∈ C, C is a white noise level line , [START_REF] Desolneux | From Gestalt Theory to Image Analysis[END_REF] which depends only on s and can be empirically estimated. Again, the curve detection algorithm consists in adequately rejecting the null hypothesis H 0 : the values of |R s | are i.i.d., extracted from a noise image. We assume that, in the background model, contrast and regularity are independent. Definition 6. Let C be a finite set of N ll level lines of u.

A level line C ∈ C is a ε-meaningful regular boundary (MRB) if NFA R K (C) = N ll K s min k<Ks B(n • l (2s,n) , k • l (2s,n) ; H s (ρ k )) < ε, (6)
and K s is a parameter of the algorithm. We also note

k min = arg min k<Ks B(n • l (2s,n) , k • l (2s,n) ; H s (ρ k )). ( 7
)
Proposition 4. The expected number of ε-meaningful contrasted regular boundaries in a finite set E of random curves is smaller than ε.

A proof is given in [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF]. We can also combine the two criteria, as done by Tepper et al. [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF], obtaining the following formulation.

Definition 7. Let C be a finite set of N ll level lines of u. A level line C ∈ C is a ε-meaningful contrasted regular boundary (MCRB) if NFA CR K (C) = N ll K c K s max   min k<Kc I c (C, k) 2 min k<Ks I s (C, k) 2   < ε, (8)
where

I c (C, k) = B(n • l (2,n) , k • l (2,n) ; H c (µ k )), I s (C, k) = B(n • l (2s,n) , k • l (2s,n) ; H s (ρ k )),
and K c and K s are parameters of the algorithm. In this case, k min is defined as in Eq. 3 or as in Eq. 7, depending on whether contrast or regularity attain the maximum.

Here K c and K s have the same meaning as K in Definition 4 and they are also set as a percentile of the total number of points in the curve.

Proposition 5. The expected number of ε-meaningful contrasted regular boundaries in a finite set E of random curves is smaller than ε.

A proof is given in [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF]. Definition 7 exhibits some interesting properties:

• A contrasted but irregular curve will not be detected;

• A regular but non-contrasted curve will not be detected;

• An irregular and non-contrasted curve will not be detected;

• A regular and contrasted curve will be detected.

Both gestalts, i.e., contrast and good continuation, interact in a novel way: they compete for the "control" of the curve.

Meaningful boundaries are further pruned to eliminate the packet effect (i.e., a bundle of level lines passing through a very salient edge) using the maximality procedure detailed in [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF].

Boundary clean-up

Prop. 3 asserts that if a level line is a meaningful boundary, then it cannot be entirely generated in white noise (up to ε false detections on the average) but it can have parts that are likely to be contained in noise.

Cao et al. [START_REF] Cao | A Theory of Shape Identification[END_REF] propose to give an upper bound to the size of those parts. Assume that C is a piece of level line with L independent points, contained in a non-edge part, described by the noise model. The probability that L is larger than l > 0 needs to be estimated, knowing that |Du| ≥ µ. This is exactly the a posteriori length distribution p(µ; l) = P (L ≥ l | |Du| ≥ µ). The estimation of this distribution was studied by Cao et al. [START_REF] Cao | A Theory of Shape Identification[END_REF].

Let us now consider an image u with N ll (quantized) level lines. Let us also denote by N l the number of all possible sampled subcurves of these level lines.

(N l = N ll i=1 n i (n i -1)/2
, where n i is the number of independent points in line i). As in Prop. 3, it can be proved that N l • p(µ; l) is an upper bound of the expected number of pieces of lines of length larger than l with gradient larger than µ. For a fixed µ, let be l such that N l •p(µ; l) < ε. Then, we know that on a white noise image, on the average, we cannot observe more than ε pieces of level line with a length larger than l and a gradient everywhere larger than µ. Then one can define L(µ) = inf{l, N l • p(µ; l) < ε} and keep every subcurve of any meaningful boundary with length equal or greater than L(µ), where |Du| ≥ µ.

The value of µ can be seen as a new parameter of the method. Its value can be fixed arbitrarily using a conservative approach [START_REF] Cao | A Theory of Shape Identification[END_REF]: letting |Du| be less than 1, means that edges with an accuracy less than one pixel may be detected. Thus, taking µ = 1 is the least restrictive choice. For µ about 1, values of L(µ) less than a few hundreds are obtained.

Since L(µ) is a decreasing function of µ, fixing it at a small value produces large lengths. We are imposing that the contrasted pieces have to be very large and this is not always the case, as argued before. Furthermore the probability distribution p(µ; l) has to be estimated. We propose to take a different path to remove non-contrasted boundary parts.

In Def. 4, pieces of a meaningful boundary are explicitly allowed to be generated in white noise. We are certainly not interested in these pieces and this relaxation responds to the fact that we want to retrieve the remaining pieces of that boundary (i.e. edge region). The desired detection of contrasted parts in a boundary is very close in spirit to periodic subsequence detection.

Algorithm 1: Computation of ε-meaningful boundaries in image u and application of the proposed clean-up procedure.

Boundary clean-up by detecting meaningful periodic subsequences

We now explain how to adapt the detector of meaningful periodic subsequences (MPS) for cleaning-up boundaries.

Before applying the detector to any boundary, we need to binarize it since its contrast (or its regularity) takes on real values. This former problem is solved by thresholding on the contrast (or on the regularity). In this direction, we claim that µ kmin and ρ kmin are natural choices (see Defs. 4, 6, and 7). A maximal ε-meaningful boundary is thus converted into a periodic binary sequence. We want to apply the periodic subsequence detection algorithm from Defs. 3 and 2 to that sequence. The only parameter left is p = Pr(s i = 1), 1 ≤ i ≤ L and it is straightforward defined, for each boundary, as p = H c (µ kmin ) or p = H s (ρ kmin ), depending on the case.

We finally define the following clean-up rule: For any meaningful boundary, keep every subcurve belonging to its maximal 1-meaningful subsequence.

Algorithm 1 shows a possible procedure to obtain all ε-meaningful boundaries and then apply the proposed clean-up method to them.

Discussion

The proposed clean-up mechanism does not impose a minimal length to contrasted parts. The length is adjusted automatically, by choosing the most meaningful subsequence in the level line. As an additional advantage, there is no need to estimate any probability distribution. Fig. 3 shows an example of the benefits of the proposed clean-up method over the one by Cao et al. [START_REF] Cao | A Theory of Shape Identification[END_REF]. Their version clearly produces underdetection: visually important structures are missed (notice the face in the third image). The proposed algorithm produces a very mild overdetection: some small noisy parts are not eliminated but no important structure is lost. Figs. 4, 6 and7 show examples on images from the BSD database [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. In all cases MPS eliminates the vast majority of spurious pieces of level lines.

Notice that, on the last row of Fig. 4, MPS does not remove a few pieces of lines that should be removed (e.g., the lower wall and the roof). This does not occur because of a failure in MPS, but because of a faulty binarization, that is, the µ kmin was not optimal in those cases. Fig. 5 depicts the same example, with a more aggressive binarization threshold. These spurious structures are now discarded. In any case, there is no universal threshold for all cases and it should be tuned for the application at hand. Figure 5: Altering the binarization threshold might help, in some cases, to eliminate spurious structures. We modified it in a simple way that allows to visualize its effect, setting for the result on the right the threshold p 2 = min(6 * p 1 , 1), where p 1 is the threshold for the result on the center. Notice that the window is better recovered with the original threshold.

Before concluding, we need to address an important objection. Instead of applying MPS to the meaningful boundaries, one could directly apply MPS to the raw collection of level lines. For each level lines, we could test several thresholds, as we do in MB, and select the periodic subsequence for which the NFA is minimized. This procedure would avoid using MB as an intermediary step. There are two main reasons for not following this more straightforward approach. First, we have a computational reason: applying MPS with different thresholds would be prohibitive; MB is much faster than MPS, allowing to render a practical algorithm. Second, since the maximality constraint [START_REF] Tepper | On the role of contrast and regularity in perceptual boundary saliency[END_REF] can be applied thanks to the inclusion properties between level lines. If we dissecting the level lines before applying the maximality rule, we cannot apply it afterwards. A different, and much more complex, mechanism would have to be used to solve the packet effect.

In summary, we presented a general and fully automatic algorithm to detect meaningful subsequences within periodic binary sequences. As a useful application, we used it to select salient pieces of level lines in an image. We have also shown that different saliency criteria can be used, such as contrast, regularity or a combination of both. In all cases the proposed algorithm exhibits good results on natural images. We recommend to zoom-in the results in the electronic version (they are vectorial graphics) for better visualization.

Figure 1 :

 1 Figure 1: A periodic sequence where runs are represented in green. If treated as a non-periodic sequence, any subsequence detector would detect four subsequences at best, when in fact the desired result is to detect three subsequences.

  [START_REF] Grompone | On Straight Line Segment Detection[END_REF] but adapted to periodic sequences. The expected number of ε-meaningful n-subsequences is given by E   (a1,...,an)∈M(n,S) 1 NFA(a1,...,an)<ε   = (a1,...,an)∈M(n,S)

Figure 3 :

 3 Figure 3: Comparison of the results obtained with both clean-up algorithms. The one by Cao et al. (CU)[START_REF] Cao | A Theory of Shape Identification[END_REF] produces underdetection; this is corrected by using MPS.

Figure 4 :

 4 Figure 4: Results of the presented clean-up algorithm on MCB. MPS eliminates the vast majority of the unwanted pieces of level line.

Figure 8

 8 Figure 8 presents a comparison of the results obtained with the different meaningful boundaries algorithms.

Figure 6 :

 6 Figure 6: Results of the presented clean-up algorithm on MRB. MPS eliminates the vast majority of the unwanted pieces of level line.

Figure 7 :

 7 Figure 7: Results of the presented clean-up algorithm on MCRB. MPS eliminates the vast majority of the unwanted pieces of level line.

Figure 8 :

 8 Figure 8: Comparison of the results obtained with the presented clean-up algorithm on the different versions of MB. MPS eliminates the vast majority of the unwanted pieces of line. We recommend to zoom-in the results in the electronic version (they are vectorial graphics) for better visualization.


Laouafa Farid
email: farid.laouafa@ineris.fr

Thermomechanical impact of Underground Coal Gasification exploitation

Keywords: Underground, coal, combustion, mining, thermomechanic, modelling

Underground mining by coal combustion (Underground Coal Gasification -UCG) raises questions regarding the mechanical behavior of the site and the stability of the overburdened rock layers. By studying the underground reactor, its inlet and outlet, the key roles played by mechanical damage and thermo-mechanical phenomena were confirmed. Deformation or collapse above the cavity may cause a collapse in the overlying layers or a subsidence at the surface level. These phenomena highly depend on the thermo-poro-mechanical behaviour of the surrounding rock (host rocks). The numerical results presented in this paper were derived from models based on different assumptions describing a raw geological background. Both 3D and 2D nonlinear finite element modelling was conducted based on two different approaches. Based on the analysis of the numerical results, we were able to highlight the main factors influencing the behaviour and mechanical stability of the overburden, and consequently the UCG process evolution.

INTRODUCTION

Underground Coal Gasification (UCG) is a process in which coal is burnt in situ without an explicit extraction of the matter [START_REF] Bell | Coal Gasification And Its Applications[END_REF][START_REF] Bhutto | Underground Coal Gasification: From fundamentals to applications[END_REF][START_REF] Shafirovich | Underground Coal Gasification: A Brief Review of Current Status[END_REF]. The efficiency of the process depends on the combustion properties of coal rock and other thermodynamic features of the physical phenomenon (temperature, pressure, mass transport related variable) [START_REF] Williams | Combustion Theory The Fundamental Theory of Chemically Reacting Flow Systems Second Edition[END_REF][START_REF] Speight | The Chemistry and Technology of Coal[END_REF].

In terms of cavity shape and its time evolution, the UCG resembles the longwall mining method [START_REF] Peng | Roof Stability in longwall coal face[END_REF][START_REF] Peng | Longwall Mining[END_REF][START_REF] Peng | Longwall Mining[END_REF]Shabanimashcool, 2012).

The UCG process involves high temperatures of around 1000°C, which will affect the host rock, changing its mechanical properties. Based on different tests and research of mudstone specimens from coal seams, [START_REF] Blackwell Luo | High-Temperature Mechanical Properties of Mudstone in the Process of Underground Coal Gasification[END_REF] shows that high temperatures cause thermal expansion of mudstone, and that the expansion increases with temperature, producing expansive stress in the surrounding rock of underground gasifiers. Brittle cracks are the main type of deformation. The temperature dependency of rock behavior has been the subject of intense research (Bauer 1983;Deng'gina et al. 1994;[START_REF] Dwivedia | Thermo-mechanical properties of Indian and other granites[END_REF][START_REF] Fjaer | Fissure formation and subsurface subsidence in a coalbed fire[END_REF][START_REF] Fjaer | Fissure formation and subsurface subsidence in a coalbed fire[END_REF][START_REF] Somerton | Thermal properties and temperature-related behavior of rock/fluid systems[END_REF][START_REF] Wolf | The Interaction between Coal Fires and Their Roof Rocks[END_REF].

The UCG process may also have other environmental impacts in terms of pollution of groundwater and surface subsidence. The groundwater may be affected by gas or others chemical species transported through the rock porous matrix, and the initial or induced cracks and microcracks. At the ground surface, soil subsidence is expected. The extent of subsidence is dependent on the location of the coal seam, its thickness, and the mechanical properties of the overburden.

The purpose of our article is related to the thermomechanical effects induced by the UCG process.

In this framework we performed finite element analyses of a case study coal site that could be exploited by this method. The numerical results presented in this paper were derived from models based on different assumptions describing a raw geological background. Both 3D and 2D nonlinear finite element modelling was conducted based on two approaches. The first one assumed the rock media as a perfect thermo-elastoplastic continuum. In the second one, and in order to simulate explicitly large space scale crack propagation, a method based on finite element deactivation was developed, which is built on a finite element mesh refinement and uses two failure criteria.

On the basis of the analysis of the numerical results, we were able to highlight two main factors that influence the behavior and mechanical stability of overburden, and consequently the sustainability of the UCG process.

UCG PROCESS

In the case of UCG, the shape of the cavity (first factor) and its evolution with time can be strongly modified by the thermo-mechanical behaviour of the host rocks. The second factor is the presence of a heat source whose location and intensity evolve with time. The thermal properties of rocks depend strongly on their mineralogical composition, porosity (φ), saturation, etc [START_REF] Fjaer | Fissure formation and subsurface subsidence in a coalbed fire[END_REF]. Table 1 below gives values of porosity (φ), thermal conductivity (λ) and specific heat capacity (ρ c p ), where ρ is the material density. [START_REF] Somerton | Thermal properties and temperature-related behavior of rock/fluid systems[END_REF] has shown that for different rock materials, the thermal properties are not constant but temperature dependent.

Table 1: Thermal properties of some materials [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF].

In Figure 1, we describe the Underground Coal Gasification (UCG) process. The injection of air is performed by a vertical well, while syngas (synthetic gas) resulting from coal combustion is extracted by the other well. The coal seams are made to react underground with insufficient oxygen for complete combustion to create syngas. Air is usually used as an oxidant. The oxygen and water within the coal seam react with the coal to produce syngas, which is withdrawn through a production well.

Figure 1 depicts the so called CRIP (Controlled, Retractable Injection Point) method, in which the injection point is gradually withdrawn as coal is consumed. This specific UCG industrial process uses two parallel wells also called Knife-edge CRIP. In this case the distance between the two wells axes is taken to be 20 m. The thickness of the coal seam is equal to 17.3 m. It is assumed that the total width of the burning front is equal to 30 m. The temperature at the combustion front is about 1000 °C (after Total, a French multinational integrated oil and gas company). The burning front is assumed to be a straight plane at 1000 °C and its normal velocity is about 0.3 mm/min (after Total Company). The temperature within the cavity can reach several hundred degrees Celsius.

Figure 2: Stratigraphy and location of coal seam (on the left) and 3D numerical model (on the right). A 2D plane passing through the middle of exploitation used in plane strain numerical analyses. V stand for the burning front velocity.

In the case of a cavity (or gallery) having a geometric configuration such as that shown above, the deformation mechanism is mainly plane strain (plane perpendicular to the axis of the cavity). The plane strain analyses are based on a cavity width of 30 m. These dimensions assume that combustion extends 5 m transversely from two casings installed 20 m away from each other.

The reader can refer to [START_REF] Fjaer | Fissure formation and subsurface subsidence in a coalbed fire[END_REF] for the mechanical and thermal properties of some of the rocks. Among all these properties that have been identified, some representative averaged values were extracted. The properties used in these models are as follows): Thermal conductivity: 2 W/m/K, Heat capacity: 800 J/kg/ K, density: 2500 kg/m3. The Young's Modulus of the coal is taken as 2.5 GPa and Poisson's ratio as 0.3. The Young's modulus of overburden is assumed to vary from 5 GPa to 25 GPa and the Poisson's ratio is equal to 0.3. Several computations have been performed by varying these parameters within realistic value ranges, and they do not change our conclusions based on our modeling work. Low values for the Young's Modulus are assigned to the layers close to the surface in order to consider them as dead load. The thermal expansion coefficient of rocks is equal to 10 -5 °C-1 .

2D AND 3D THERMOMECHNICAL MODELLING

The objective is to analyze the spatial extent of plasticity (damage) in the vicinity of the cavity. The model size is deliberately reduced to the dimensions defined in Figure 3 dt with

ε ε ε ε ε = = ∫    
We observe that plasticity affects a part of the cover and the coal wall. However, the extension is limited but the coal layer is also affected.

We can assume that the hypotheses of material continuity and classical thermo-elasticoplastic framework are too conservative. In other words, if the analysis was performed in the framework of fracture mechanics, linear or nonlinear [START_REF] Broek | Elementary Engineering Fracture Mechanics[END_REF], we would certainly have found a lower overall stability rather than local rockfall from the roof and major collapse in coal (Figure 3, right). Likewise, it is reasonable to assume that overburden massive coal generally consists of calcareous clay layers, sometimes sandy, is layered, heterogeneous, and marked by discontinuities at different scales. The overall strength properties are quite low. Therefore, the choice of pair of values (cohesion, friction angle) representative of the whole overburden is critical for the prediction of the extent of the damaged area. We can a priori assume that the uniaxial compression strength of the overburden seldom exceeds 3 to 5 MPa overall. This leads to an extension of the damage larger than calculated above. These failure types can affect the combustion process as the combustion front may break and, due to the failures at the edge of the cavity, change the kinetics, the size, and shape of the combustion front, and therefore the shape of the cavity. They can also affect the direction of the cavity evolution. It is observed that the greatest displacement naturally lies in the cavity roof. The effect on the surface (subsidence) will be presented in the section dedicated to 3D models. If it is not explicitly specified, the criterion for the adaptive mesh refinement is based upon the value of the equivalent plastic strain. The model on which we performed our analysis is illustrated by its finite element mesh in Figure 4. The width of the cavity is 30 m. This approach required the use of another finite element code that allows adaptive mesh refinement and the possibility of deactivating a number of elements (MARC TM finite elements code, MSC Software, 2014). Figure 4 shows also spatial distribution of plasticity and Figure 5 presents the mesh density when using finite element refinement and deactivation in order to simulate crack propagation.

It was observed that the extension of the plasticity zone for a cavity of 30 m width is limited in space but still significant.

Figure 5: Example of result with adaptive mesh refinement and element deactivation (L = 30 m). The effective plastic deformation range from 0 to 0.17.

The results in Figure 4 are obtained through a classical approach and the elastoplastic continuum medium (i.e., overburden and coal) is assumed to be continuous and remains so regardless of the degree of "damage". Figures 5,6, and 7 show the same computation, but using a deactivation (roughly simulating fracturing) of finite elements based on adaptive mesh refinement. We observe that with the most intense element concentration (small elements and therefore finer mesh), the damage zone becomes more detailed. Note also that the computation highlight the "failure" in the walls of the cavity and the collapse of part of the roof of the cavity.

On the basis of the spatial distribution of the mesh, it is easy to foresee the shape and damaged areas. The extent of mesh refinement is limited only by the power of the computer. The refinement is an h-refinement (based on subdivision of the finite element), which is different than p-refinement (based on the degree of the finite element).

Initial Damaged

The effect of the fineness of the mesh is shown in Figures 6 and7.

The rock properties (reduction of cohesion and friction angle by 40%) were decreased and the cavity was widened to 30 m. The initial mesh was regular and relatively coarse (in size). Figure 6 shows the distribution of fracturing in the overburden, in the coal walls, and the mesh density in the plasticity zones. Figure 7 (right) shows a focus of quality of mesh refinement required in order to carry out the calculation. The "triangular" shape of the collapse area in the roof and high degradation of coal wall are to be noted. This failure mode, described here in a very simplified way, is commonly observed in situ in the case of underground mines or quarries. The figures represent results at an instantaneous state (observation at a given time t) because the elements are allowed to move in the dynamic analysis. The inertial terms induced by gravity and the fractures play an important role and this analysis therefore differs fundamentally from the quasi-static approach. The cavity extends laterally until the coal zone is sufficiently undamaged (by failure). The combustion process which is already a chaotic phenomenon would be difficult to formalize in such condition. The combination of burning and fracturing/block falling should influence the rate of combustion in the cavity. The modelling results depicted in Figures 6 and7 present the main failure mechanisms as described by Whittles et al. (2007). The coal wall, the roof, and part of the floor are damaged by the cavity creation. This numerical computation using adaptive mesh refinement and element deactivation seems relevant for this study as they contain much information performed even on an idealized geometry such as failure and damage of the roof and wall being important outcomes. It should be noted that such mechanisms are those generally observed during underground mining. In the case of UCG, these failures can disrupt the function of the reactor, alter the path of the combustion front, and so on.

SUBSIDENCE AND ROOF DEFORMATION

The 3D computations are performed on the geometric model shown in Figure 8. The cavity and its spatio-temporal evolution are not explicitly created. More precisely, the initial cavity, a purely geometrical shape, has the dimensions of the final "cavity". However, it applies to the boundaries of the cavity pressure and temperature whose laws of spatial and temporal evolution describe the expansion of the cavity. The first remark concerns the vertical displacement of the roof of the cavity. Asymmetry (with respect to the reference solution T=0) occurs. Note in the modelling the "front" of combustion moves from left to right. The higher the temperature, the greater the intensity of movement (left increases), and the asymmetry is more pronounced. This is explained by the fact that by conduction and due to the intensity of the temperature difference, the part to the left is longer subjected to a high temperature. Consequently the extension of the temperature will be greater in the rock matrix. There is also a shift to the right for the maximum subsidence. The surface subsidence is induced not only by the excavation of the cavity, but also by thermal effects. That the subsidence is negative indicates that the contribution of the rock mass thermal expansion in the uprising is lower than the displacement induced by stress modification due to underground excavation. The following 3D results are given as an illustration and show the burning front evolution. Elastoplastic models with adaptive refinement and 3D finite element deactivation were performed in order to assess the consequences of a possible "fracturing".

Figure 10: Example of induced "fracturing" during time and the space evolution of the cavity.

STRATIFIED ROOF

In this short subsection we analyze briefly the behavior of a stratified roof composed of marl and sandstone (Figure 11). The elliptical cavity is located in the coal layer. he remaining is sandstones.

Figure 11: Mesh of the elliptical burning cavity and stratified roof.

After a given extend of the burning cavity (assumed elliptical in this sub-section), we applied a given temperature to its boundary (200 °C). After a given time we obtained the following displacement field as illustrated in Figure 12. One can observe that the intensity of the vertical displacement and the final shape is strongly dependent on the assumption of a perfect continuum. The effect of the presence of layers with different mechanical parameters is show in Figure 13. The spatial distribution of the effective plastic strain is totally different when we take into account a layered roof. In this configuration that exhibits a stratified roof we highlight the transient effect of temperature. Indeed, thermal effects are more important when the medium is stratified and thus heterogeneous. The differential thermal and/or mechanical properties induce incremental strains and shear at interfaces layers. The case presented here is an academic and theoretical case. In real environments, we encounter heterogeneities at the porous matrix level, where it is well known and demonstrated that they induce an increase of the Von Mises stress. In addition, the results for a stratified medium show that solving such problems by using continuous media assumptions is far too conservative.

CONCLUSIONS

Geomechanical evaluations of the rock mass properties remain very important processes during Underground Coal Gasification. An accurate stability analysis of the consequence of UCG requires a global characterization of the rock constituting the overburden. All of the different factors contribute to instability of the host rock during and after the gasification process. Unlike the stresses induced by void creation and their redistribution (as encountered in conventional mining), UCG introduced additional thermal stresses in the vicinity of the reactor. Thermal expansion may induce buckling or elastoplastic instability due to the increase of compressive stresses. The problem, which can be defined as transient bifurcation, is quite complex. The shape of the cavity and its evolution in time remain an important parameter which influences the stability of the host rock. Nevertheless, the evolution of the cavity in space and time has some random character due to the non-control of the parameters driving the combustion process and also due to the effects of host rock behavior. Caution must be taken when using a specific technique for UCG simulation.

Figure 1 :

 1 Figure 1: Schematic illustration of the adopted UCG industrial process with two wells.

  . Although the model has a vertical plane of symmetry (through the axis of the cavity), the calculations were carried out on the whole 2D model to analyze a component of the robustness of the numerical model (maintaining the symmetry of the mechanical response).

Figure 3 :

 3 Figure 3: Spatial distribution of the effective plastic strain (in the left) and schematic of potential movements (failure, rockfalls, collapse,..) on the right.

Figure 4 :

 4 Figure 4: 2D finite elements mesh (left) and effective plastic strain (right) (L=30 m).

Figure 6 :

 6 Figure 6: Initial mesh (left) and after dynamic computation, using adaptive mesh refinement and finite element deactivation (right). (Cavity width=30 m).

Figure 7 :

 7 Figure 7: Focus of the coal wall and roof meshes (several thousand elements).

Figure 8 :

 8 Figure 8: 3D mesh and location of the coal seam.

Figure 9 :

 9 Figure 9: Vertical displacement of the roof of the cavity (segment CD) when it is created (left) and surface subsidence (segment AB) when the whole cavity is created (right) for 4 temperatures fields inside the cavity.

Figure 12 :

 12 Figure 12: Vertical displacements field using two methods: Pseudo fracturing (top) and continuum model (bottom).

Figure 13 :

 13 Figure 13: Effective plastic strain distribution with crack propagation.

ACKNOWLEDGEMENT

The work performed in this paper is supported by the EU project COGAR, Grant Agreement RFCR-CT-2013-00002.


Mountaka Souley
email: mountaka.souley@ineris.fr

Marwan Al Heib

Vincent Renaud

Numerical investigation of EDZ development around a deep polymetallic ore mine

Keywords: Brittle-ductile behaviour, constitutive model, verification, Garpenberg mine, numerical modelling

This paper deals with the development of a non-linear constitutive model of rock mass and its verification to predict a damaged zone. Simulations of triaxial compressions provide a verification of the implementation with a good agreement between predictions and theoretical values of peak and residual strengths as well as the transition between brittle failure and ductile response. The applicability of the model to predict potential failure around stopes of a deep polymetallic ore mine is checked. The present study highlights the interest to consider more realistic rheology of hard rock masses compared with the elastic perfectly plastic models of underground deep mines.

INTRODUCTION

Damage by microcracking is the main dissipation process associated with inelastic behaviour and failure in most brittle materials such as rocks, concrete and ceramic composites. Under high in situ and induced stresses and high anisotropic stress ratios, an excavation damaged zone (EDZ) may be formed around underground openings excavated into brittle rocks. The failure mechanism in this damaged zone is the initiation, growth and coalescence of cracks and fractures and is directly related to the constitutive behaviour of rock mass. Experimental studies on brittle rocks have shown that there are many different mechanisms by which cracks can be initiated and grown under compressive stresses [START_REF] Steif | Crack extension under compressive loading[END_REF][START_REF] Horii | Brittle failure in compression: splitting, faulting and brittleductile transition[END_REF]. The involved mechanisms include sliding along pre-existing cracks and grain boundaries, pore crushing, elastic mismatch between mineral grains and dislocation movement, etc. Indeed, irreversible deformations and failure of rocks subjected to stresses occur through progressive damage as microcracks initiate and grow at small scale and coalesce to form large-scale fractures and faults. Over the past decade, several constitutive models have been developed providing a more realistic description of damage for brittle rocks in relation to experimental observations [START_REF] Zhu | Homogenization-based analysis of anisotropic damage in brittle materials with unilateral effect and interactions between microcracks[END_REF]).

In the framework of the European program I 2 Mine (http://www.i2mine.eu/), INERIS has developed an innovative tool (stress monitoring tool) aimed to better understand the behaviour of the rock mass resulting from the mining process of Garpenberg deep mine (superposed rooms and pillars and vertical stopes between depths 1100 and 1300 m, Figure 1) in Sweden. Several seismic sensors and strain gauges were installed in a deep area of the mine (approximately at 1150 m deep) to monitor the evolution of induced stresses during the exploitation and the related induced seismic activities.

The objective of this paper is to study the mechanical behaviour of the main galleries and the pillars by incorporating the stages of the underground work through advancing numerical modelling. At first approach, we are limited to 2D models that take into account the initial state of stress (anisotropic) and the mine operations.

More precisely, the purpose of this paper is to present: (a) a numerical implementation of an elastoplastic model obeying to the Hoek-Brown criterion and taking into account the brittle/ductile transition between the brittle failure and the ductile behaviour depending on the mean stress as generally observed on most rock samples: that is to say, the initiation and growth of cracks for brittle rocks are modeled by a softening behaviour in the post-peak in the framework of plasticity theory, (b) the corresponding verification based on simulation of triaxial compression tests. Finally, the numerical modelling provides the localization of EDZ and the potential instability area in the vicinity of stopes located around the pilot area of the mine works at the Garpenberg mine (New Boliden company).

PROPOSED RHEOLOGICAL MODEL AND NUMERICAL IMPLEMENTATION

The mechanical behaviour of geomaterials is widely varied and depends mainly on the confining stress (or mean stress) and the loading paths. At low stress confining levels, rocks break by the creation of one or more shear planes or bands accompanied by a strain-softening behaviour characterized by microcracks dilatancy and grains rotation at the microscopic scale. Under high mean stresses, rocks undergo a hardening behaviour which is microscopically associated with volumetric strains. More precisely, for most geomaterials subjected to triaxial compressions under low to moderate confining pressures, it is generally observed the following typical characteristics of stress-strain curves: (a) a linear isotropic and elastic behaviour after a short non-linear phase corresponding to the closure of the initial pores and the progressive contact walls of the original micropores and microcracks (b) a strain-hardening in the pre-peak region corresponding to the initiation and the growth of microcracks modeled generally by the plasticity theory or the damage mechanics through the concept of effective stress and the hypothesis of strain equivalence (c) a strain-softening after reaching the peak strength (failure) associated with a progressive loss in material cohesion and then a decrease in strength (d) a phase where the rock strength remains practically constant. Relative to the elasto-brittle materials where the initial microcracks closure and the growth of microcracks in the pre-peak domain can be neglected in a first approximation (as it is the typical case of hard rock masses studied herein at the Garpenberg mine), the present study considers an isotropic elastic linear behaviour before the peak.

The Hoek-Brown criterion widely used in the field of rock mechanics and rock engineering, whose limitations were documented by detailed discussions of the simplifying assumptions [START_REF] Hoek | Underground excavations in rock[END_REF], was generalized in terms of the three stress invariants in this study. One of the most important limitations is the non-dependency of the criterion on the intermediate principal stress, σ 2 . Subsequent experimental studies have since suggested that the intermediate principal stress has a substantial influence on the rock strength (e.g. [START_REF] Haimson | True triaxial stresses and the brittle fracture of rock[END_REF]. This has led to the development of several 3D versions of the Hoek-Brown failure criterion [START_REF] Priest | Determination of shear strength and three-dimensional yield strength for the Hoek-Brown criterion[END_REF][START_REF] Zhang | A generalized threedimensional Hoek-Brown strength criterion[END_REF].

It is well known that the three principal stresses σ 1 , σ 2 and σ 3 can alternatively be expressed in terms of mean stress (p), generalized deviatoric stress (q) and Lode's angle (θ) and vice versa, which allows to express the classical Hoek-Brown criterion in 3Dstresses space (p, q, θ). After some manipulations, the generalized yield function is proposed on the basis of [START_REF] Souley | Modelling of the viscoplastic behaviour including damage for deep argillaceous rocks[END_REF]:

𝐹 𝑠 = 4 3 cos 2 𝜃 𝜎 𝑐 𝑞 2 -𝑚 � cos 𝜃 √3 - sin 𝜃 3 � 𝑞 + 𝑚 𝑝 -𝑠𝜎 𝑐 (1)
where σ c is the uniaxial compressive strength of the intact rock, m, s represents the rheological behaviour parameters whose significance and quantification differ from the dimensionless empirical parameters of Hoek-Brown reflecting the strength of the intact rock, the natural fractured rock mass and therefore consideration of a possible scale effect and fracturing at the large scale. The parameter m quantifies the frictional strength and s (measuring the induced cracks) is related to the rock cohesion. Let (m p , s p ) and (m r , s r ) be respectively the values of parameters (m, s) at the peak and the beginning residual phase. Assuming a transition between brittle and ductile behaviour, equation (1) makes it possible to express a relationship between the peak and residual strength parameters:

𝑚 𝑟 = 𝑚 𝑝 + �𝑠 𝑝 -𝑠 𝑟 � 𝜎 𝑐 𝜎 3 𝑏-𝑑 (2)
In the absence of experimental data on the residual strength, almost null value of s r can be considered as a good approximation in relation to sliding along the macroscopic failure planes or shear bands. For softening phase, the same yield function (equation 1) is used by introducing the softening internal variable (ξ). In accordance with the yield function used, (a) the generalized plastic strain is chosen to represent the internal variable, (b) m linearly varies with ξ between the peak and the beginning of residual phase (0 < ξ ≤ ξ r representing the shear plastic strain) whereas s varies non linearly (c) m and s remain constant in the residual domain:

𝑑𝜉 = � 2 3 𝑑𝑒 𝑝 : 𝑑𝑒 𝑝 𝑚 = 𝑚 𝑝 + �𝑚 𝑟 -𝑚 𝑝 � 𝜉 𝜉 𝑟 (3) 𝑠 = 𝑠 𝑝 + �𝑠 𝑟 -𝑠 𝑝 � �2 - 𝜉 𝜉 𝑟 � 𝜉 𝜉 𝑟 (4)
where de p is increment of deviatoric plastic strain tensor in post-peak or residual regions.

Assuming for this type of rock materials that the plastic dilatancy occurs only near the peak strength and cannot exceed a threshold rate (β m ) for a certain value of deformation and by adapting a non associated flow rule (generally accepted for rock materials), we propose the following plastic potential (G s ) which allows to express the constitutive relationship between the increments of stresses (dσ) and total strains as follows (dε):

𝐺 𝑠 = 𝑞 -𝛽 𝑚 �1 -𝑒 -𝑏 1 𝜉 �𝑝 (5) 𝑑𝜎 = �𝐶 - �𝐶: 𝜕𝐹 𝑠 𝜕𝜎 � ⨂ �𝐶: 𝜕𝐹 𝑠 𝜕𝜎 � 𝜕𝐹 𝑠 𝜕𝜎 : 𝐶: 𝜕𝐺 𝑠 𝜕𝜎 - 𝜕𝐹 𝑠 𝜕𝜉 𝜕𝐺 𝑠 𝜕𝑞 � 𝑑𝜀 (6)
In addition to the two elastic constants: Young The most input data used in this paper are from laboratory tests, field studies, reports and papers. More precisely, the Young's modulus (E), Poisson's ratio (ν) and the unconfined compressive strength (σ c ) have been obtained from the point and biaxial load tests or uniaxial compression loading carried out on samples (complex ores containing zinc, lead, silver, copper and gold) taken from the measuring holes at -883 m and -967 m levels [START_REF] References Edelbro | Strength, fallouts and numerical modeling of hard rock masses[END_REF][START_REF] Van Koppen | Estimation of the Risk for Mine Induced Seismicity in Large Scale Mining in the Garpenberg Mine[END_REF]. These parameters and the empirical constants m p and s p , and dilatancy at failure were confirmed by many numerical studies aimed to the calibration of the mechanical characteristics, the determination of the in situ stresses (at -880 m level and from one access gallery at -1155 m level, respectively carried out by Boliden and INERIS) and the understanding of the failure mechanism (for instance spalling, shear failure, slabbing, buckling or other types of compressive failure) as well as the back-evaluation of fallout extent observed in the field.

It is known that it is often very difficult to get the full stress-strain curves of the post-peak behaviour particularly for hard rock masses exhibiting brittle behaviour. As first approach, we retained the values of additional parameters (s r , σ 3 b-d , ξ r , β m and b 1 )

summarized in Table 1: a sensitivity analysis on these parameters will be conducted in parallel to the interpretation of the in situ stress measurements at 1155 m coming in the future and whose deformation measurement systems have already been installed from an access gallery at -1155 m level (Figure 1c).

NUMERICAL IMPLEMENTATION AND MODEL VERIFICATIONS

In order to verify the implemented model, seven triaxial compression tests with confining pressures ranged from 0.01 to 75 MPa have been simulated. The input parameters used for this exercise are summarized in Table 1.

Figure 2 presents the deviatoric stress versus axial and lateral strains curves for different confining pressure (σ 3 ). We note that the post-peak behaviour is confining pressure dependent: the transition stress between brittle failure and ductile behaviour is clearly marked and the numerical transition stress, σ 3 b-d , occurs for a confining pressure of 63 MPa in accordance with the input numerical value.

Figure 3 shows a comparison in terms of the peak and residual strengths between the predictions and theoretical expression in equation (1) (with θ = 30° for triaxial compression loading path). The match is very good as may be seen in this figure, where numerical and analytical solutions coincide.

The relative error for peak and residual strengths is less than 0.3 % and 0.7 % respectively for peak and residual strengths. This validates the numerical implementation of the proposed elastoplastic model in FLAC 3D .

APPLICATION TO GARPENBERG MINE

The aim of this section is to provide a verification of the implementation for non-triaxial stress paths and to show numerically the ability of the model to evaluate the extent of potential damaged zones in relation to the excavation of one or more stopes.

In the Garpenberg mine, polymetallic ore is extracted which contains zinc, silver, lead, copper and gold. The ore has fold axis and the depth is 500-1300 m below the surface. During the mining operations of the upper level (-883 m and -967 m levels), several field studies and numerical modellings were performed with the aim to calibrate certain mechanical characteristics, back-determine the initial and induced in situ stresses and to understand the potential failure mechanisms. The results of these observations and measurements were also used as input data for the preliminary 3D modelling of the initial stress field prevailing at the lower level around 1155 m deep before the future excavation of stopes and for which stress measurements are available using CSIRO cells (Figure 4).

To retrieve the initial stress state (at 1155 m), a 3D numerical model has been carried out with several assumptions. A specific methodology for verifying the correlation of measured stresses (where locations are illustrated in Figure 4) with computed stresses (after mining operations of the upper level) is developed: it consists of creating a stress data bank with the numerical model according to the range of measured stress tensors. Because the dependence between initial and induced stresses is strongly nonlinear, it is very difficult to find what is the initial stress tensor that produces the correct measured stresses after the excavation of many galleries and stopes and after a partial backfilling of the stopes (chambers). Nevertheless, it was possible to estimate the initial stress tensor corresponding to level -1155 m. The horizontal and vertical directions are principal with the following principal values 47.5 MPa (x-horizontal direction), 55 MPa (yhorizontal direction) and 20.5 MPa (vertically). This stress state was used to uniformly load the numerical 2D-model shown in Figure 4.

As illustrated in Figure 1c, to model the mine geometry in 2D as first stage, we can consider a double symmetry in the vertical and horizontal directions. The geometrical model for 2D-modelling with FLAC 3D assuming plane strain conditions is shown in Figure 4. Null normal displacements are imposed at the boundary limits (symmetry). The access galleries are 10 m wide and 6 m high. Finally, typical primary stope sizes are: 80 m long (along Ydirection), 10 m wide, and 19 m high.

The modelling sequences were performed as follows: (a) model without excavations was consolidated under the previous initial in situ stresses, using roller boundaries to the model sides respectively parallel to x-and z-axis for seeking symmetry (step 0); (b) four galleries are simultaneously excavated, numerical model is stepped until equilibrium (step 1) and (c) stopes 1 and 2 are simultaneously deleted and the numerical model is stepped until equilibrium.

The plastic zones corresponding to the galleries excavation (phase 1) (or failed with behaviour in residual or transition between peak and residual) are limited to the galleries skin. This is consistent with the observed fallouts (whose extension is around 0.05 m) or predicted by numerical analysis conducted at -880 m level. This is also consistent with the strength of the rock masses (σ c √s≈63 MPa) and the magnitude of induced stresses. In considering for example the case of the gallery sides and assuming the linear elasticity behaviour and neglecting variations of axial stress, the expected orthoradial, radial, and axial After excavating the two stopes in a 2D model, the extent of plastic zones is concentrated in the pillar (Figure 5), in the compressive region where the maximum deviatoric stress is located (at the heart of pillars between stopes along the symmetry axis). Figure 5a shows the distribution of plastic zones in the pillars (and therefore strain localization). In fact, the most compressive stresses are located at the singularity levels (corners of galleries), but as the mean stress also increases in magnitude, it is the same for the failure strength. Excavation of stopes results in complete release of the horizontal stress, causing predominantly tensile stresses over the whole width and substantially all the height of pillars. Therefore, the failures can be initiated at the center of pillars and propagate along an inclined plane at first and then sub-vertically following the localization of the plastic zones. As a consequence of stress projections on the current yield envelop, we observed decreases in the vertical and axial stresses. This is well illustrated in Figure 6, which shows the distribution of horizontal and vertical stresses at the end of step 2.

For comparison purpose and in order to capture the contribution of the proposed rheological model (#1), another simulation is performed with a basic perfectly plastic model with Hoek-Brown's criterion (model #2). In this case, the plastic zones are also developed in the heart of the pillars but in a very limited area compared to the results of the model #1 (previously presented) without plastic strain localization (Figure 5b). Specifically, the examination of the principal stresses (σ 1 , σ 2 , σ 3 ) at the pillar center delimited by the two stopes for the two models 1 and 2 provides the following stress values: (-20.5, -17.5, 0.75 MPa) and (-48, -35.5, 0.85 MPa) respectively for #1 and #2 (compressive stresses < 0). Finally, Figure 6b clearly shows a significant increase in vertical stress (σ zz ) near the changes of orientations of failure plane locations, indicating some eventual seats of potential risks of instability. The quantitative comparison between the two models was conducted on the stress values at the end of step 2 (points 1 to 4 illustrated in Figure 4) and reported in Table 1 in terms of stress ratios with respect to the initial stress state before the excavations of galleries and stopes.

One can note that in both cases, the four points are outside the plastic domain induced by galleries and stopes. The points P 1 , P 2 , and P 3 are located in the pillar between the two lower galleries and point P 4 is located at the floor of a stope. Table 2 reveals a slight difference between the models on the relaxation at P 1 , P 2 , and P 3 and compression at P4 of the horizontal stresses (along x and perpendicular to the model cross section). Conversely, over-stresses develop in the vertical direction at P 1 , P 2 , and P 3 positions with more magnitude (15 to 35%) for #1 compared to #2, resulting from the failure occurred a little earlier in the pillar of stopes, and then a new stresses redistribution.

CONCLUSION

This paper presents the development of a constitutive model for brittle hard rock masses, its numerical implementation in the three-dimensional code FLAC 3D , as well as its verification. Simulation of triaxial compression tests at different levels of confining pressure provides a verification of the implemented model. The resulting curves display three regions (elastic, softening in post-peak, and residual phase) as well as a transition between brittle failure and ductile behaviour as generally observed on hard and brittle rocks. The ability of the proposed model to predict the potential failure regions (damaged zones) around deep underground excavations in hard rocks is successfully tested. In particular, the importance of considering the postpeak behaviour to predict strains and failure zones localization, and consequently the advantage over the approaches without strain softening is shown.

It is often difficult to characterize the post-peak behaviour parameters from laboratory tests on such hard rocks in relation to the brittle nature of their behaviour. These parameters can be viewed as calibration parameters, which can be backdetermined on the base of the in situ observations. The prospect of this work after a full 3D modelling back-analysis would complete this numerical development by implementing the energy balance in view of assessing the proneness to rockbursts or strain bursts.

Figure 1 :

 1 Figure 1: Experimental set-up (stress and microseismic measurements) in the Swedish mine Garpenberg: (a-left) the mine location; (b-centre) operating plan and location of the pilot area; (c-right) 3D mesh of galleries and stopes.

  modulus and Poisson ratio (E, ν), the needed parameters for the proposed model are: m p , m r , σ 3 b-d , s p , s r , σ c , ξ r , β m and b 1 . All of these parameters can be identifiable from triaxial compression tests (3 minimum) carried out from low confining pressure (or without) to high confining pressure (characterizing σ 3 b-d ). Finally the proposed model was implemented in the three-dimensional explicit finite-difference code, FLAC 3D , where a perfect elastoplastic model with the Hoek-Brown criterion in the (σ 1 , σ 3 ) plane already exists in the commercial version.

Figure 2 :

 2 Figure 2: Numerical results of triaxial compression tests.

Figure 3 :

 3 Figure 3: Peak and residual strengths: numerical and analytical solutions.

  16, 0, and 50 MPa. This confirms the failure at the walls. Moreover, away from the wall, the presence of a non null radial stress significantly increases the rock mass strength by the term (m p σ r σ c ), which explains the limited character of failed area around the galleries.

Figure 4 :

 4 Figure 4: Geometrical model with position of history points.

Figure 5 :

 5 Figure 5: Extent of plastic zones around the stopes: (a) model #1 (b) model #2.

Figure 6 :

 6 Figure 6: Horizontal (top) and vertical (down) stresses: model #1.

Table 1 :

 1 Numerical values of input model parameters.

	Parameter	Value	Parameter	Value
	E (MPa)	80000	s r (-)	10 -5 s p

Table 2 :

 2 Ratio of horizontal (σ xx ), vertical (σ zz ) and axial (σ yy ) stresses.Model

		σ xx * P 1	σ zz * P 1	σ yy * P 1
	#1	0.38	1.45	0.93
	#2	0.41	1.59	0.94
	Model	σ xx * P 2	σ zz * P 2	σ yy * P 2
	#1	0.39	1.31	0.92
	#2	0.42	1.57	0.94
	Model	σ xx * P 3	σ zz * P 3	σ yy * P 3
	#1	0.41	1.22	0.91
	#2	0.41	1.58	0.94
	Model	σ xx * P 4	σ zz * P 4	σ yy * P 4
	#1	1.04	0.96	1
	#2	1.05	0.87	1

(a) (MPa) (b) (MPa)