- Email: [email protected]
3D volumetric modelling of Cadomian terranes (Northern Brittany, France): an automatic method using Voronoı̈ diagrams
Tectonophysics 331 (2001) 181±196
www.elsevier.com/locate/tecto
3D volumetric modelling of Cadomian terranes (Northern Brittany, France): an automatic method using VoronoõÈ diagrams Gabriel Courrioux a,*, SteÂphane Nullans b, Antonio Guillen a, Jean Daniel Boissonnat b, Philippe Repusseau a, Xavier Renaud a, Muriel Thibaut a b
a BRGM B.P. 6009, 45060 OrleÂans Cedex 02, France INRIA 2004, Route des Lucioles, B.P. 93, 06902 Sophia-Antipolis Cedex, France
Received 1 December 1997; accepted 1 February 2000
Abstract An automatic method based on the use of VoronoõÈ diagrams is proposed for the 3D volumetric reconstruction of geological objects. It enables volumes to be constructed starting from the combined use of geological maps and cross-sections de®ned in multiple directions, and can also take into account incomplete information on geological interfaces and faults. The method is suitable for global modelling of a set of geological objects; it gives a consistent 3D volumic partition of space according to geological information contained in the maps and the cross-sections. The method is described and applied to modelling the main domains of the Cadomian collisional orogenic belt of Panafrican age in Northern Brittany (France). Thanks to the contrasting densities of the main units involved in the collisional process, the gravimetric effect of the volumic model can be calculated. Its comparison with the observed anomaly allows us to discuss the validity of the geological model. Although some improvements appear to be necessary, it is concluded that the method can make geological modelling more ef®cient through providing the possibility of rapidly testing many hypotheses. q 2001 Elsevier Science B.V. All rights reserved. Keywords: 3D modelling; VoronoõÈ diagrams; Cadomian belt; Brittany
1. Introduction The Cadomian orogenic belt of of Northern Brittany (France) is interpreted classically in terms of the evolution of an active margin ending with the collision of a continental margin, a marginal basin and a volcanic arc (Chantraine et al., 2001) from 650 to 450 Ma (Panafrican orogeny). A three-dimensional representation of the main units and discontinuities of the North Armorican basement is a major objective of the Armor GeoFrance 3D programme, since one of the major problems for understanding the Cadomian * Corresponding author. Fax: 133-2-38-64-33-61. E-mail address: [email protected] (G. Courrioux).
orogeny concerns the geometry of these elements at depth. As any 3D modelling must enable different geological hypotheses to be illustrated and the validity of cross-sections to be checked, ef®cient methods of volume reconstruction must be developed that take into account geological data diversity. The general goal, therefore, is to achieve models that are (1) geometrically correct Ð i.e. ®t known geometrical features, (2) topologically consistent Ð i.e. respect the relationships between different components of a geological object, (3) geologically realistic, and (4) reproducible. In addition, the underlying data structure must be suitable for the end applications, such as geophysical modelling. We propose a method based on VoronoõÈ diagrams, which enables volumes to be
0040-1951/01/$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0040-195 1(00)00242-0
182
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
constructed starting from the combined use of geological maps and cross-sections de®ned in multiple directions, and which can also take into account incomplete information on geological interfaces and faults drawn on the cross-sections. The aim of this paper is to describe the method and to discuss its application to the geological data of the Cadomian belt in Northern Brittany (France). The resulting 3D model corresponds to an a priori model based on geological maps and cross-sections. 2. Review of current methods With regard to sur®cial interpolations, classical geostatistical methods (Blanchin and Chiles, 1993) have proved their ef®ciency when suf®cient data points are known on relatively simple geological interfaces. Interpolation methods such as BeÂzier parametric surfaces (Farin, 1988; Auerbach and Schaeben, 1990) and triangulated surfaces (Mallet, 1992, 1997) also have been applied to geometric surface modelling (Siehl et al., 1992; Renard and Courrioux, 1994). These methods all start from geological interface datasets and are aimed at modelling single surfaces. Meshed surfaces can handle uncertainty as a property assigned to vertices; however, to date, only geostatistical methods can give a good quanti®cation of uncertainty on the result. Different approaches are currently being explored for 3D solid modelling (Foley et al., 1995): (1) Surface boundary representations. Volumes are built by sewing meshed or parametric surfaces derived from surface reconstruction methods (Bertrand et al., 1992; Halbwachs et al., 1996). This approach requires fairly large interface datasets with homogeneous distributions. It is dif®cult to treat complex geometries with various topologies when few data are known. The major problem here is that surfaces are constructed independently of each other and so it is dif®cult to guarantee their correct relationships, and to arrive at a consistent ®nal topology. This approach therefore requires many intersection calculations between all the surfaces. (2) Volumetric meshing. Volumes are described by sets of small elementary cubes (voxels) each of them being assigned a property. The voxel values are computed from interpolating discrete values
(Mayoraz, 1993). This results in interpolating a Boolean value that describes a point as being inside or outside a given formation. A volume is then described as the set of voxels having the same property. This method can be used provided that appropriate geostatistical tools are available. It requires a large amount of memory space. (3) Volume reconstruction starting from crosssections (Boissonnat, 1988; Meyers et al., 1991). Current algorithms, mainly developed for medical applications, require complete, close and more or less parallel sections. For such algorithms to be transposed to geology requires that the interfaces be completely described, that they bound entire formations on each of the sections, and that these sections do not intersect one another. Such conditions are far from respected in natural case studies; thus these methods are unable to combine the geological map and partial cross-section data provided by geologists. In this paper we propose an automatic construction based on the VoronoõÈ cell method where volumes are described as sets of connected polygonal shapes (polytopes). In that sense, it can be viewed as a combination of boundary representations and pure voxel methods since a voxel is a particular case of polytope. The following constraints on data must be satis®ed: ² formations are sampled from the interior at known data points; ² some data on the geological interfaces are known; ² data are indistinctly located on any support, i.e. maps, cross-sections, or drill-holes; ² interfaces and formations can be incompletely described on any support; ² uncertainty on the location of data points must be taken into account; ² faults can be taken into account. After describing the principle of the method, we illustrate different aspects through its application to Cadomian geological data. 3. Principle of the method 3.1. The VoronoõÈ diagram and the data discretization Taking P {p1 ; p2 ; ::::::::::; pn } as a ®nite set of points, called sites, of the Euclidean space E 3, then the VoronoõÈ diagram of these sites is the partition of
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
the space assigning each point to its nearest site (see Fig. 1). The VoronoõÈ diagram consists of n cells, one per site, so that the VoronoõÈ cell Ci consists of all the points at least as close to pi ; as to any other site: Ci Vpi {x [ E3 ; ;q [ P=d x; pi # d x; q} where d is the Euclidean distance in E 3. An ambiguous case may occur if four points are co-cyclic in 2D (®ve co-spheric points in 3D), then
Fig. 1. Principle of VoronoõÈ reconstruction. (a) VoronoõÈ diagram and Delaunay triangulation in two dimensions. (b) Construction of boundary between two regions A and B. The boundary is a subset of the VoronoõÈ diagram. The union of circum disks to Delaunay triangles gives the space where this boundary is allowed to move while smoothing.
183
two possible triangulations are possible that satisfy the de®nition. This indetermination can be raised by slightly perturbing the position of sites which makes the triangulation non-ambiguous (Seidel, 1994). The dual graph of the VoronoõÈ diagram is obtained by using segments to connect the pairs of site belonging to adjacent VoronoõÈ cells (Fig. 1a). The dual graph of the VoronoõÈ diagram is a triangulation called the Delaunay triangulation; it has the property that the circumcircle of every triangle contains no sites in its interior (Fig. 1b). In practice, a VoronoõÈ model is described as a set of cells. Each cell surrounds one of the sites and is made up of connected polygonal faces. Each face of a cell is the mediator of one edge of the Delaunay triangulation. Obviously two adjacent cells share a common face. Implementation principles and data structures are reviewed by Aurenhammer (1991). The VoronoõÈ diagram and the Delaunay triangulation of a set of n points can be constructed optimally in O n log n time. Based on these principles, an incremental algorithm using a tree data structure for Delaunay tetrahedrons is used (Boissonnat and Teillaud, 1993). It allows us to insert new sites ef®ciently. Although incremental algorithms may be quadratic in the worst case, this algorithm requires expected O n log n time for inserting a new site, which is optimal. Boissonnat and Nullans (1996) provide a detailed description of the method and its application to geological objects. Initial data may sample either the interior of geological formations or their boundaries. The geometric elements that bear these two kinds of information are: ² Isolated points Ð outcrop observations provide an indication either on the interior of a geological formation or, in some cases, on a contact between two formations. ² Segment lines or curves Ð drilling data enable the interior of a formation to be sampled along a curved line, whereas geological maps and crosssections give an indication on the trace of interfaces between two formations. ² Surfaces representing already modelled interfaces Ð a formation from the geological map which is projected onto a digital elevation model (DEM) represents the 3D interface between the given formation and the atmosphere.
184
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
Fig. 2. Data discretization and generation of points: (a) on each side of a given interface along the drilling line, (b) on each side of an interface separating two geological formations, (c) on each side of the surface separating two formations according to the slope, (d) following an observed contact according to structural data, and (e) according to uncertainty on data.
All these data types have to be transformed into a discrete set of points that all belong to the interior of a given formation, i.e. any data point must belong to a given formation. These will be used as VoronoõÈ sites. This is done in such a way that successive points on a discretized curve or surface belong to adjacent cells in the VoronoõÈ diagram (or equivalently, are joined by an edge in the Delaunay triangulation). This is always possible provided a suf®cient number of points are
taken (Boissonnat, 1988; Edelsbrunner and Tan, 1992). Data are discretized as follows (Fig. 2): ² Drilling lines are discretized according to a given step (Fig. 2a). ² Points on interfaces at formation boundaries are discretized according to a given step along the boundary. At each step two points are generated
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
on each side of the interface (Fig. 2b), at a distance from the contact relative to the precision at which the contact is known. The vector de®ned by these two points is normal to the known azimuth and dip of the interface at this place (Fig. 2d). ² Surfaces representing an interface are discretized in their interior and the points are then duplicated on each side of the interface according to the local normal at the surface (Fig. 2c).
185
² When a surface is at the boundary of the model (geological map), a particular formation which represents the exterior of the concerned zone needs to be introduced. ² Each point is then assigned to its geological formation. ² The uncertainty on data can be taken into account by de®ning the displacement value of the duplicated points on each side of an interface in relation
Fig. 3. Construction of homogeneous regions in 2D. (a) Initial heterogeneous data set: some drilling data and pieces of interfaces are known at some places. (b) VoronoõÈ diagram built from the points data set. (c) Extraction of boundaries according to the assigned formations at VoronoõÈ sites. (d) Smoothing boundaries and taking faults into account.
186
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
to this uncertainty (Fig. 2e). If an interface point pi is known with an uncertainty d, then discretization must be such that the distance between pia and pib is at least equal to d, ipia 2 pib i ! 2 £ d; and such that the direction of vector piapib is the normal direction at this place. If the normal direction is unknown, points are generated within the support they belong to (maps or sections or drill-holes). If an interface point pi is precisely known d is assumed to be very close to 0.
3.2. Taking faults into account Two methods are being explored for taking faults into account: (1) Faults are inserted into the model after VoronoõÈ reconstruction, and formations are remeshed according to the intersections between the VoronoõÈ cells and the faults. All intersection vertices between faults and formations are split into two nodes, one at each side of the fault. All split nodes are constrained to move only along the fault surface during the smoothing process. If the throw function of a fault is known, then it can be introduced as an additional constraint. (2) A fault is modelled as a surface de®ned by a set of discrete points (Fig. 6) so that it becomes an inherent part of the VoronoõÈ diagram. The main difference is that for a geological interface one knows which formation is on either side, whereas in the case of a fault this is generally unknown. Therefore the VoronoõÈ cells around the faults are assigned an indeterminate colour (Fig. 6a). The method then consists of correctly reassigning colours according to the neighbouring cells. The whole diagram is constructed in a single phase with fault and formation datasets. The colours around faults are then assigned according to the maximum common surfaces with neighbouring formations (Fig. 6b). Though the ®rst method has been successfully applied in some two-dimensional cases (Fig. 3), it raises some problems in three dimensions for the following reasons: ² intersection problems are much more dif®cult to solve in three dimensions especially handling with tangency; ² when large displacements occur along faults, the initial solution (i.e. without faults) may be very far
Fig. 4. Splitting VoronoõÈ cells. Node 1 with neighbours 2±5 is split into two nodes with neighbours 2±4 and 3±5.
from what should be the actual solution, and there is then little chance to readjust properly the geological formations along the faults; ² the initial VoronoõÈ partition is lost due to intersections, consequently the model is no more suitable for insertion of new sites, and queries on point location cannot be performed on the model. The second method is therefore favoured for its algorithmic simplicity and for its more general application.
Fig. 5. Curvature of the mesh on a node Pi. An approximation of the curvature of the mesh is 1/R, R being the radius of the circumcircle of Pi, Pi1, Pi2 and Pi3. PiNi is the normal to the plane given by Pi1, Pi2 and Pi3. Position of Pi is moved along the normal PiNi according to smoothing criteria.
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
187
² geometric regularity constraint Ð smoothing must minimize some ®xed bending; ² topology constraint Ð smoothing must respect the initial topology of the reconstruction.
Vi 1; ¼; 4 in the graph of the mesh. Accounting for the geometry regularity constraint is performed locally on the current vertex and three neighbours by using the radius of the circumsphere as minimizing criteria. It is therefore necessary to split the nodes of degree 4 into two nodes of degree 3 (Fig. 4). All the vertices of the initial VoronoõÈ mesh are moved until an equilibrium is reached (the equilibrium criteria is that the total displacement of the vertices is less than a given value (e). After smoothing, the vertices of an initial polygonal face are no longer coplanar. They are therefore triangulated.
Since smoothing results in moving the vertices of the initial VoronoõÈ cells, care must be taken that data points de®ned inside geological formations remain inside after the smoothing process. From the VoronoõÈ diagram, the boundaries of geological formations are extracted as polygonal shapes. These can be seen as 3D meshes by considering their edges and vertices where each vertex V has at least three, and at most four, neighbouring vertices.
3.3.1. Geometric regularity constraint The curvature of a curve (surface) is a local property which quanti®es how much it bends. A geometrical de®nition can be given: the largest circle (sphere) that is tangential to a curve (surface) on its concave side at a point has the same curvature than the curve at that point. On each vertex we can de®ne a discrete approximation
3.3. Smoothing process Because of the roughness of the VoronoõÈ reconstruction (Fig. 3c), the shapes of the interfaces do not look quite ªnaturalº. It is therefore necessary to smooth the boundaries according to the following constraints:
Fig. 6. Taking faults into account. (a) VoronoõÈ sites with an indeterminate colour are generated on each side of the fault. (b) VoronoõÈ diagram and boundary extraction. The cells around the fault are assigned a colour according to the maximum common surfaces with neighbouring formations.
188
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
189
of the normal and the curvature of the mesh M. The neighbours Pi1, Pi2 and Pi3 de®ne a plane PLi and the normal vector Ni of PLi will be the normal of the mesh at Pi (Fig. 5). An approximation of the curvature of the mesh can be given by 1/R, R being the radius of the circumsphere of Pi, Pi1, Pi2 and Pi3. Indeed, the larger the radius, the closer the points to a planar distribution. It is possible to adapt the position of Pi with respect to a given curvature, and the position of the three neighbours. The smoothing method consists of obtaining a more regular distribution of curvatures along all points of the mesh. This can be achieved by interpolating a new curvature in Pi from the curvatures of the three neighbours and then to recalculate Pi the interpolated curvature and three neighbours. This is done on all vertices of the mesh one after the other in random order. This process can be iterated several times.
regions can be obtained (Fig. 1b). The surfaces bounding two volumes are described as sets of connected polygonal faces. A two-dimensional example of the reconstruction process is illustrated in Fig. 3. The regions are a ®rst approximation of the geological formations obtained when applying the following rule: a point p x; y; z belongs to a formation F if the site closest to p x; y; z belongs to F. This rule does not take into account the heterogeneous spatial distribution of the data. However, it is used to provide a solution that is, in most cases, topologically correct, i.e. homotopic to the actual geological model. In order to improve the quality of modelling, several constraints may be added in this reconstruction algorithm; for example, minimizing or maximizing the volume of a formation, preserving the thickness of some formations, or taking into account the anisotropy of the medium.
3.3.2. Topology constraint: Delaunay ball constraint This is the most particular constraint of the reconstruction method. It allows the topology of the different shapes of the model to be respected during the smoothing step. The Delaunay triangulation has the property that the circumcircle of every triangle contains no sites in its interior; such a circle (or disk) will be simply called an empty disk for short. Maintaining the displacement of the mesh edges and vertices in the union of these empty disks will guarantee that the initial topology is respected and that data points are inside their proper volumes.
4. 3D Armor model
3.3.3. Construction of homogeneous regions After construction of the VoronoõÈ diagram, each site pi [ P has one corresponding VoronoõÈ cell, so that each cell can be assigned to the geological formation (indexed by a colour) of its corresponding site. By merging adjacent VoronoõÈ cells that have the same colour, a partition of the space into colour-connected
4.1. Geological context The Cadomian orogenic belt of northern Brittany (France) is classically interpreted in terms of the evolution of an active margin ending with the collision of a continental margin, a marginal basin and a volcanic arc (Chantraine et al., 1998). Metamorphic, magmatic and structural evolution resulting from this collision is recorded from 650 to 540 Ma (Panafrican orogeny). Shortening resulted in the accretion of different units and was accommodated mainly by the combination of a southwest-directed wrench and strike-slip tectonism (Bale and Brun, 1989). The convergence kinematics and the way in which the different units are involved is still being debated (Egal et al., 1996; Strachan et al., 1996). Especially with regard to the geometry some questions arise: Do the structures ¯atten at depth? How are they branched? What is the relative importance of the
Fig. 7. Simpli®ed geological map of the main Cadomian units and geological cross-sections. Section 1: x1 217:643; y1 2395:084; x2 220:000; y2 2430:004: Section 2: x1 210:626; y1 2398:404; x2 226:505; y2 2423:857: Section 3: x1 203:488; y1 2408:585; x2 219:493; y2 2427:790: Section 4: x1 200:261; y1 2411:856; x2 216:931; y2 2437:662: Section 5: x1 201:749; y1 2410:827; x2 231:645; y2 2413:321: Section 6: x1 191:249; y1 2417:121; x2 231:168; y2 2414:590: Section 7: x1 230:706; y1 2383:129; x2 222:016; y2 2422:174: Section 8: x1 242:037; y1 2386:942; x2 219:850; y2 2420:225: Section 9: x1 258:209; y1 2386:664; x2 236:678; y2 2414:258: Coordinates are expressed in kilometres (kilometric Lambert2 coordinates system).
190
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
different units? The realization of a 3D geometric model serving as a reference model, allowing quantitative analysis, is a necessary starting point for answering these questions. Moreover, the contrast between physical and mechanical properties of the blocks makes an appropriate area for a testing methodology based on combined volumic and gravity modelling. 4.2. Lithology Units that have been modelled here are (Chantraine et al., 2001): (1) The Saint-Brieuc Unit composed of (1) the high metamorphic complex of the Yf®niac Formation, (2) Pentevrian orthogneisses including the Port-Morvan Formation, (3) dominantly basic metavolcanics of the Lanvollon Formation, (4) clastic metasediments of the Binic Formation overlying the metavolcanics, (5) basic plutonic intrusions (Squif®ec, St-Quay and Fort-Lalatte plutons). (2) The Saint-Malo±Guingamp Unit, originally composed of clastic brioverian sediments, characterized by the widespread development of migmatites. (3) Paleozoic granites intruding the Brioverian sediments in the southern part of the map. (4) Ordovician overstep sequence, composed of clastic sediments. Cadomian Units (Saint-Brieuc and Saint-Malo± Guingamp) are bounded by major faults along which they were accreted to the Central Armorican continental block.
4.3. Data Geometric modelling is based on a geological map and nine cross-sections oriented radially with respect to the arcuate shape of the belt (Fig. 7). The different units are discretized in terms of points belonging to each unit so that they can be suitable for VoronoõÈ reconstruction. In the particular case of the studied area, the major faults can be considered as geological interfaces since they form major boundaries between distinct units. Traces of geological contacts are transformed into points on each side of the contact. Assuming that sections are nearly perpendicular to the main structures, these points are generated within the section normal to the contact trace. They are then assigned their respective geological formation. The way contacts are discretized is illustrated in Fig. 8. Point density is globally equal to one point per 1 km along contacts but the distribution may vary according to the size and complexity of structures. Special care must be taken that point distribution re¯ects the characteristic features of the sections as, for example, triple points or thin narrow elongated bodies. Where a contact has an uncertain position or is poorly constrained, as is the case with most of the contacts at depth, the points are generated farther away from the drawn contact. 4.4. Results Starting from a given data set, the VoronoõÈ method allows one to automatically model the volume geometry of different geological units, while ensuring
Fig. 8. Example of discretization on Section 9. Points density is globally equal to one point per 1 km along contacts but their distribution may vary according to the size and complexity of structures. Care must be taken that points distribution re¯ects the characteristic features of the sections as for example triple points or thin narrow elongated bodies. Where a contact has an uncertain position or is poorly constrained, as is the case with most of the contacts at depth, the points are generated farther away from the drawn contact.
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
191
Fig. 9. 3D VoronoõÈ diagram of the Cadomian belt. (a) Volumic construction of VoronoõÈ diagram with no smoothing; (b) smoothing boundaries; and (c) extraction and visualization of different volumes from the model. The model is made up of different volumes that share common boundaries. 1 Ð Ordovician overstep sequence; 2 Ð Paleozoic granite intrusion; 3 Ð the Saint-Malo±Guingamp units; 4 Ð Pentevrian orthogneiss including the Port-Morvan Formation; 5 Ð undifferentiated Lanvollon and Binic formations.
192
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
topologic consistency between data and volumes. About 5 min on a current work-station are required to provide a model with 2000 data points. The principle of the method is illustrated through modelling the bulk architecture of the belt. In a ®rst step, ®ve volumes were constructed, corresponding to the Palaeozoic granites, the Saint-Malo±Guingamp Unit, Saint-Brieuc Unit where the volcano-sedimentary basin (Lanvollon±Binic) has been distinguished, and Ordovician sediments. Fig. 9a shows the resulting VoronoõÈ diagram without any smoothing. Though this representation ®ts the data, it remains geologically unrealistic. Fig. 9b shows the same volumes after smoothing, resulting in a more realistic model. Smoothing is controlled in order to respect the topological relations between data points and volumes. The arcuate shape of the belt is well illustrated. Major contact units can be extracted from the volumes, showing that the model is not only a set of
single volumes (Fig. 9c) but can be considered as a whole entity containing volumes sharing common boundaries. Fig. 10 shows a more detailed representation of the Saint-Brieuc Unit. The volumic imagery of these different formations makes the model comprehensive and communicable, helping to pose problems by highlighting inconsistencies. As example, southward dips of Port-Morvan Orthogneiss in Eastern part of the model brings some torsion in the shape of the model that suggests that Eastern and Western parts of the model are not continuously and simply correatable. Adding or correcting data in the model is easily done by simply inserting new sites in the diagram. As an example, Fig. 11 shows two models of the Binic Formation and the Ordovician overstep sequence. The ®rst model was poorly constrained in the north-eastern part while some fuzzy information on the estimated depth of the interface between the
Fig. 10. Modelling geological formations of the Cadomian belt. 1 Ð Ordovician overstep sequence; 2 Ð Paleozoic granite intrusion; 3a Ð the Saint-Malo±Guingamp units including 3b Ð the Quessoy orthogneiss; the following formations of the Saint-Brieuc Unit are represented: 4a Ð Pentevrian orthogneiss including the Port-Morvan Formation; 4b Ð basic intrusions including the Saint-Quay gabbrodiorite; 4c Ð the Yf®niac Formation; 4d Ð clastic sediments of the Binic Formation.
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
Fig. 11. Binic Formation and Ordovician overstep sequence. (a) The model for the north-eastern part was poorly constrained. Fuzzy information on the estimated depth of the interface has been added to the second model. 1 Ð Ordovician overstep sequence; 2 Ð clastic sediments of the Binic Formation.
two formations has been added in the second model, leading to a more geologically reasonable model. It is clear that the method only re¯ects the geometry of the underlying geological model contained in the cross-sections. However, some quanti®cation on the three-dimensional consequences of a hypothesis can be performed. For example, the simulation of geophysical effect can be computed based on the assumed physical properties of the medium. Since modelled volumes are closed and de®ne a whole partition of the space, the classical method of Okabe (1979) has been used to compute the gravity effect of the model (Fig. 12). When compared to the observed anomaly, discrepancies are analyzed in order to discuss the geological uncertainties and eventually to converge
193
to the best satisfying model. It can be seen that, compared to the observed one (Truffert et al., 2001), the ®rst-order variations are in agreement but that negative anomalies are overestimated. This is probably due to an insuf®cient precision in Hercynian granite de®nition. The depth and thickness of Port-Morvan orthogneiss are probably under-estimated. Important discrepancies exist at the north-east of the model owing to the lack of data. Great uncertainty exists with regard to the extension of Binic Serie and Lanvollon Serie in the North East under the Ordovician overstep sequence. This area clearly needs to be better constrained. Another reason is that the chosen formation regroupments may not be the best ones with regard to gravity properties. So, new sections calibrated both on geological constraints and on gravity anomaly pro®les have been redrawn. These were used to reconstruct the model and led to a more satisfying model. This latter model and the comparison between the computed gravity effect and observed anomaly are discussed in Truffert et al. (2001). However, a remaining problem concerns how and where to interact on the initial sections or on the model to reduce discrepancies. Some solutions are explored using the Monte-Carlo inverse method to randomly walk around the model space solution (Mosegaard and Tarantola, 1995). 5. Discussion The ®rst interesting aspect of the method is mainly its global character: because all volumes are built simultaneously, it is not necessary to construct each geological interface separately and independently to achieve a complete model, as with most surface methods. Taking into account data on the interior of the formations and not only on their contacts is a great advantage. The simplicity of the method allows the geologist to focus on data speci®cation and validity rather than on fastidious 3D modelling reconstruction methods. Only 3D points and their geological assignment are necessary to build a model. The second interesting aspect is the automatic character of the method. The current methodology for geological modelling consists of testing geological hypotheses and comparing their gravimetric or
194
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
Fig. 12. Gravity anomaly computed from the geometric model. Coordinates are expressed in kilometres (kilometric Lambert2 coordinates system).
magnetic effect against the recorded anomalies. It is clear that an automatic model greatly facilitates testing many hypotheses. The general idea that 3D modelling allows as many hypotheses as required to be tested is therefore given real meaning. However, several points still need to be improved for this idea to become a true reality. ² The ®rst point concerns the data management facilities used for building the models. Testing many hypotheses requires tools that allow easy modi®cation and editing of geological cross-sections and maps. Work is currently in progress on this aspect. ² The second point concerns the method itself. The present model is much too sensitive to the way data have been discretized; for units that are too ªthinº in places, as for example the Yf®niac Formation, the volume is discontinuous where it is expected to be
continuous. This is due to the algorithm that assigns the formation to the nearest neighbour. In order to improve the quality of the results, two approaches are currently being investigated: ± incorporating an anisotropic rule to allow the direction and strain of geological formations to be taken into account; ± improving the discretization methods in order to respect geological constraints (hypotheses on layer continuity, thickness, strain data, anisotropy). Another improvement of the approach will be an automatic modi®cation of the geometric model using the inverse problem to ®t geophysical measurements. This can be done easily because modifying cells or adding new sites in the VoronoõÈ diagram is rapid and so there is no need to rebuild the whole VoronoõÈ.
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
Consequently, there is no need to recompute the whole anomaly for each incremental modi®cation. 6. Conclusions The principal result of this work is the demonstration that volume reconstruction using VoronoõÈ diagrams is suitable for geological objects. The volume reconstructions are automatic, reproducible and give a consistent partition of space according to the data speci®cations. This can make the geological modelling process more ef®cient by providing the possibility of testing many hypotheses rapidly. In the case of the Cadomian belt of northern Brittany, the model that was built illustrates one of the possible hypotheses on the architecture of the belt. Other models based on different interpretations also have been tested using VoronoõÈ reconstruction in this issue (Truffert et al., 2001). Moreover, in the light of this work, 3D inversion algorithms for ®tting geophysical data onto the data structure of the VoronoõÈ diagrams can be engaged positively. Acknowledgements This work was funded by the GeÂofrance3D programme and by a BRGM research project. We are grateful to Jean Chantraine and Emmanuel Egal for their involvement in this work. We acknowledge Delphine Rouby and Jan-Diederik van Wees for their constructive review of the manuscript. GeoFrance 3D publication 19. References Auerbach, S., Schaeben, H., 1990. Computer-aided geometric design of geologic surfaces and bodies. Math. Geol. 22 (8), 957±987. Aurenhammer, F., 1991. VoronoõÈ diagrams: a survey of a fundamental geometric data structure. ACM Comput. Surv. 23, 345±405. BaleÂ, P., Brun, J.P., 1989. Late Precambrian thrust and wrench zones in Northern Brittany (France). J. Struct. Geol. 11 (4), 391±405. Bertrand, P., Dufourd, J.F., FrancËon, J., Lienhardt, P., 1992. ModeÂlisation volumique aÁ base topologique: actes MICAD 92 Paris 1, 59±74. Blanchin, R., Chiles, J.P., 1993. The Channel tunnel: geostatistical
195
prediction of the geological conditions and its validation by the reality. Math. Geol. 25 (7), 963±974. Boissonnat, J.D., 1988. Shape reconstruction from planar crosssections. Comput. Vis. Graph, Image Process. 44 (1), 1±29. Boissonnat, J.D., Nullans, S., 1996. Reconstruction of geological structures from heterogeneous and sparse data. Research report no. 3069, INRIA, France, 24 pp. Boissonnat, J.D., Teillaud, M., 1993. On the randomized construction of the Delaunay tree. Theor. Comput. Sci. 112, 339±354. Chantraine, J., Chauvel, J.J., BaleÂ, P., Denis, E., Rabu, D., 1998. Le BrioveÂrien (ProteÂrozoõÈque supeÂrieur aÁ terminal) et l'orogeÁne cadomien en Bretagne (France). Bull. Soc. GeÂol. Fr. IV, 815±829. Chantraine, J., Egal, E., ThieÂblemont, D., Le Goff, E., Guerrot, C., BalleÁvre, M., Guennoc, P., 2001. The Cadomian active margin (North Armorican Massif, France): a segment of the North Atlantic Panafrican Belt. Tectonophysics 331 (1±2), 1±18. Edelsbrunner, H., Tan, T.S., 1992. An upper bound for conforming Delaunay triangulations. Proc. 8th Annu. ACM Symp. Comput. Geom., pp. 53±62. Egal, E., Guerrot, C., Le Goff, E., ThieÂblemont, D., Chantraine, J., 1996. The Cadomian orogeny revisited in northern Brittany. Avalonian and Related Peri-Gondwana Terranes of the Circum-North Atlantic, Geol. Soc. Am. Spec. Publ. 304, 281±318. Farin, G., 1988. Curves and Surfaces for Computer Aided Design, A Practical Guide. Academic Press, San Diego 401 pp. Foley, J., van Dam, A., Feiner, S., Hughes, J., 1995. C Edition Ð Interactive Computer Graphics: Principles and Practice. Addison-Wesley, Reading, MA 1174 pp. Halbwachs, Y., Courrioux, G., Renaud, X., Repusseau, P., 1996. Topological and geometric characterisation of fault networks using 3-dimensional generalized maps. Math. Geol. 28, 625±656. Mallet, J.L., 1992. Discrete smooth interpolation in geometric modelling. Comput. Aided Des. 24, 178±191. Mallet, J.L., 1997. Discrete modeling for natural objects. Math. Geol. 29 (2), 199±219. Mayoraz, R., 1993. ModeÂlisation et visualisation infographiques tridimensionnelles de structures et proprieÂteÂs geÂologiques. TheÁse, Sciences, EPF Lausanne, No. 1127, XX, 218 pp. Meyers, D., Skinner, S., Sloan, K., 1991. Surfaces form contours: the correspondence and branching problems. Graphics Interface 91, 246±254. Mosegaard, K., Tarantola, A., 1995. Montecarlo sampling of solutions to inverse problems. J. Geophys. Res. 100, 12,431±12,447. Okabe, M., 1979. Analytical expressions for gravity anomalies due to homogeneous polyhedral bodies and translations into magnetic anomalies. Geophysics 44, 730±741. Renard, P., Courrioux, G., 1994. 3-D geometric modeling of a faulted domain: the Soultz horst example (Alsace, France). Comput. Geosci. 20, 1379±1390. Seidel, R., 1994. The nature and meaning of perturbations in geometric computing. Proc. 11th Symp. Theor. Aspects Comput. Sci. (STACS), Lecture Notes in Computer Science, vol. 775. Springer, Berlin, pp. 3±17.
196
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
Siehl, A., Ruber, O., Valdivia-Machego, M., Klaff, J., 1992. Geological maps derived from interactive spatial modeling. Geol. Jahrb. 122, 273±279. Strachan, R.A., D'Lemos, R.S., Dallmeyer, R.D., 1996. Neoproterozoic evolution of an active plate margin: North Armorican Massif, France. Avalonian and Related Peri-Gondwana
Terranes of the Circum-North Atlantic, Geol. Soc. Am. Spec. Publ. 304, 319±332. Truffert, C., Egal, E., Le Goff, E., Courrioux, G., Guennoc, P., 2001. Gravity modelling of the Cadomian active margin of northern Brittany. Tectonophysics 331 (1±2), 81±97.
www.elsevier.com/locate/tecto
3D volumetric modelling of Cadomian terranes (Northern Brittany, France): an automatic method using VoronoõÈ diagrams Gabriel Courrioux a,*, SteÂphane Nullans b, Antonio Guillen a, Jean Daniel Boissonnat b, Philippe Repusseau a, Xavier Renaud a, Muriel Thibaut a b
a BRGM B.P. 6009, 45060 OrleÂans Cedex 02, France INRIA 2004, Route des Lucioles, B.P. 93, 06902 Sophia-Antipolis Cedex, France
Received 1 December 1997; accepted 1 February 2000
Abstract An automatic method based on the use of VoronoõÈ diagrams is proposed for the 3D volumetric reconstruction of geological objects. It enables volumes to be constructed starting from the combined use of geological maps and cross-sections de®ned in multiple directions, and can also take into account incomplete information on geological interfaces and faults. The method is suitable for global modelling of a set of geological objects; it gives a consistent 3D volumic partition of space according to geological information contained in the maps and the cross-sections. The method is described and applied to modelling the main domains of the Cadomian collisional orogenic belt of Panafrican age in Northern Brittany (France). Thanks to the contrasting densities of the main units involved in the collisional process, the gravimetric effect of the volumic model can be calculated. Its comparison with the observed anomaly allows us to discuss the validity of the geological model. Although some improvements appear to be necessary, it is concluded that the method can make geological modelling more ef®cient through providing the possibility of rapidly testing many hypotheses. q 2001 Elsevier Science B.V. All rights reserved. Keywords: 3D modelling; VoronoõÈ diagrams; Cadomian belt; Brittany
1. Introduction The Cadomian orogenic belt of of Northern Brittany (France) is interpreted classically in terms of the evolution of an active margin ending with the collision of a continental margin, a marginal basin and a volcanic arc (Chantraine et al., 2001) from 650 to 450 Ma (Panafrican orogeny). A three-dimensional representation of the main units and discontinuities of the North Armorican basement is a major objective of the Armor GeoFrance 3D programme, since one of the major problems for understanding the Cadomian * Corresponding author. Fax: 133-2-38-64-33-61. E-mail address: [email protected] (G. Courrioux).
orogeny concerns the geometry of these elements at depth. As any 3D modelling must enable different geological hypotheses to be illustrated and the validity of cross-sections to be checked, ef®cient methods of volume reconstruction must be developed that take into account geological data diversity. The general goal, therefore, is to achieve models that are (1) geometrically correct Ð i.e. ®t known geometrical features, (2) topologically consistent Ð i.e. respect the relationships between different components of a geological object, (3) geologically realistic, and (4) reproducible. In addition, the underlying data structure must be suitable for the end applications, such as geophysical modelling. We propose a method based on VoronoõÈ diagrams, which enables volumes to be
0040-1951/01/$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0040-195 1(00)00242-0
182
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
constructed starting from the combined use of geological maps and cross-sections de®ned in multiple directions, and which can also take into account incomplete information on geological interfaces and faults drawn on the cross-sections. The aim of this paper is to describe the method and to discuss its application to the geological data of the Cadomian belt in Northern Brittany (France). The resulting 3D model corresponds to an a priori model based on geological maps and cross-sections. 2. Review of current methods With regard to sur®cial interpolations, classical geostatistical methods (Blanchin and Chiles, 1993) have proved their ef®ciency when suf®cient data points are known on relatively simple geological interfaces. Interpolation methods such as BeÂzier parametric surfaces (Farin, 1988; Auerbach and Schaeben, 1990) and triangulated surfaces (Mallet, 1992, 1997) also have been applied to geometric surface modelling (Siehl et al., 1992; Renard and Courrioux, 1994). These methods all start from geological interface datasets and are aimed at modelling single surfaces. Meshed surfaces can handle uncertainty as a property assigned to vertices; however, to date, only geostatistical methods can give a good quanti®cation of uncertainty on the result. Different approaches are currently being explored for 3D solid modelling (Foley et al., 1995): (1) Surface boundary representations. Volumes are built by sewing meshed or parametric surfaces derived from surface reconstruction methods (Bertrand et al., 1992; Halbwachs et al., 1996). This approach requires fairly large interface datasets with homogeneous distributions. It is dif®cult to treat complex geometries with various topologies when few data are known. The major problem here is that surfaces are constructed independently of each other and so it is dif®cult to guarantee their correct relationships, and to arrive at a consistent ®nal topology. This approach therefore requires many intersection calculations between all the surfaces. (2) Volumetric meshing. Volumes are described by sets of small elementary cubes (voxels) each of them being assigned a property. The voxel values are computed from interpolating discrete values
(Mayoraz, 1993). This results in interpolating a Boolean value that describes a point as being inside or outside a given formation. A volume is then described as the set of voxels having the same property. This method can be used provided that appropriate geostatistical tools are available. It requires a large amount of memory space. (3) Volume reconstruction starting from crosssections (Boissonnat, 1988; Meyers et al., 1991). Current algorithms, mainly developed for medical applications, require complete, close and more or less parallel sections. For such algorithms to be transposed to geology requires that the interfaces be completely described, that they bound entire formations on each of the sections, and that these sections do not intersect one another. Such conditions are far from respected in natural case studies; thus these methods are unable to combine the geological map and partial cross-section data provided by geologists. In this paper we propose an automatic construction based on the VoronoõÈ cell method where volumes are described as sets of connected polygonal shapes (polytopes). In that sense, it can be viewed as a combination of boundary representations and pure voxel methods since a voxel is a particular case of polytope. The following constraints on data must be satis®ed: ² formations are sampled from the interior at known data points; ² some data on the geological interfaces are known; ² data are indistinctly located on any support, i.e. maps, cross-sections, or drill-holes; ² interfaces and formations can be incompletely described on any support; ² uncertainty on the location of data points must be taken into account; ² faults can be taken into account. After describing the principle of the method, we illustrate different aspects through its application to Cadomian geological data. 3. Principle of the method 3.1. The VoronoõÈ diagram and the data discretization Taking P {p1 ; p2 ; ::::::::::; pn } as a ®nite set of points, called sites, of the Euclidean space E 3, then the VoronoõÈ diagram of these sites is the partition of
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
the space assigning each point to its nearest site (see Fig. 1). The VoronoõÈ diagram consists of n cells, one per site, so that the VoronoõÈ cell Ci consists of all the points at least as close to pi ; as to any other site: Ci Vpi {x [ E3 ; ;q [ P=d x; pi # d x; q} where d is the Euclidean distance in E 3. An ambiguous case may occur if four points are co-cyclic in 2D (®ve co-spheric points in 3D), then
Fig. 1. Principle of VoronoõÈ reconstruction. (a) VoronoõÈ diagram and Delaunay triangulation in two dimensions. (b) Construction of boundary between two regions A and B. The boundary is a subset of the VoronoõÈ diagram. The union of circum disks to Delaunay triangles gives the space where this boundary is allowed to move while smoothing.
183
two possible triangulations are possible that satisfy the de®nition. This indetermination can be raised by slightly perturbing the position of sites which makes the triangulation non-ambiguous (Seidel, 1994). The dual graph of the VoronoõÈ diagram is obtained by using segments to connect the pairs of site belonging to adjacent VoronoõÈ cells (Fig. 1a). The dual graph of the VoronoõÈ diagram is a triangulation called the Delaunay triangulation; it has the property that the circumcircle of every triangle contains no sites in its interior (Fig. 1b). In practice, a VoronoõÈ model is described as a set of cells. Each cell surrounds one of the sites and is made up of connected polygonal faces. Each face of a cell is the mediator of one edge of the Delaunay triangulation. Obviously two adjacent cells share a common face. Implementation principles and data structures are reviewed by Aurenhammer (1991). The VoronoõÈ diagram and the Delaunay triangulation of a set of n points can be constructed optimally in O n log n time. Based on these principles, an incremental algorithm using a tree data structure for Delaunay tetrahedrons is used (Boissonnat and Teillaud, 1993). It allows us to insert new sites ef®ciently. Although incremental algorithms may be quadratic in the worst case, this algorithm requires expected O n log n time for inserting a new site, which is optimal. Boissonnat and Nullans (1996) provide a detailed description of the method and its application to geological objects. Initial data may sample either the interior of geological formations or their boundaries. The geometric elements that bear these two kinds of information are: ² Isolated points Ð outcrop observations provide an indication either on the interior of a geological formation or, in some cases, on a contact between two formations. ² Segment lines or curves Ð drilling data enable the interior of a formation to be sampled along a curved line, whereas geological maps and crosssections give an indication on the trace of interfaces between two formations. ² Surfaces representing already modelled interfaces Ð a formation from the geological map which is projected onto a digital elevation model (DEM) represents the 3D interface between the given formation and the atmosphere.
184
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
Fig. 2. Data discretization and generation of points: (a) on each side of a given interface along the drilling line, (b) on each side of an interface separating two geological formations, (c) on each side of the surface separating two formations according to the slope, (d) following an observed contact according to structural data, and (e) according to uncertainty on data.
All these data types have to be transformed into a discrete set of points that all belong to the interior of a given formation, i.e. any data point must belong to a given formation. These will be used as VoronoõÈ sites. This is done in such a way that successive points on a discretized curve or surface belong to adjacent cells in the VoronoõÈ diagram (or equivalently, are joined by an edge in the Delaunay triangulation). This is always possible provided a suf®cient number of points are
taken (Boissonnat, 1988; Edelsbrunner and Tan, 1992). Data are discretized as follows (Fig. 2): ² Drilling lines are discretized according to a given step (Fig. 2a). ² Points on interfaces at formation boundaries are discretized according to a given step along the boundary. At each step two points are generated
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
on each side of the interface (Fig. 2b), at a distance from the contact relative to the precision at which the contact is known. The vector de®ned by these two points is normal to the known azimuth and dip of the interface at this place (Fig. 2d). ² Surfaces representing an interface are discretized in their interior and the points are then duplicated on each side of the interface according to the local normal at the surface (Fig. 2c).
185
² When a surface is at the boundary of the model (geological map), a particular formation which represents the exterior of the concerned zone needs to be introduced. ² Each point is then assigned to its geological formation. ² The uncertainty on data can be taken into account by de®ning the displacement value of the duplicated points on each side of an interface in relation
Fig. 3. Construction of homogeneous regions in 2D. (a) Initial heterogeneous data set: some drilling data and pieces of interfaces are known at some places. (b) VoronoõÈ diagram built from the points data set. (c) Extraction of boundaries according to the assigned formations at VoronoõÈ sites. (d) Smoothing boundaries and taking faults into account.
186
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
to this uncertainty (Fig. 2e). If an interface point pi is known with an uncertainty d, then discretization must be such that the distance between pia and pib is at least equal to d, ipia 2 pib i ! 2 £ d; and such that the direction of vector piapib is the normal direction at this place. If the normal direction is unknown, points are generated within the support they belong to (maps or sections or drill-holes). If an interface point pi is precisely known d is assumed to be very close to 0.
3.2. Taking faults into account Two methods are being explored for taking faults into account: (1) Faults are inserted into the model after VoronoõÈ reconstruction, and formations are remeshed according to the intersections between the VoronoõÈ cells and the faults. All intersection vertices between faults and formations are split into two nodes, one at each side of the fault. All split nodes are constrained to move only along the fault surface during the smoothing process. If the throw function of a fault is known, then it can be introduced as an additional constraint. (2) A fault is modelled as a surface de®ned by a set of discrete points (Fig. 6) so that it becomes an inherent part of the VoronoõÈ diagram. The main difference is that for a geological interface one knows which formation is on either side, whereas in the case of a fault this is generally unknown. Therefore the VoronoõÈ cells around the faults are assigned an indeterminate colour (Fig. 6a). The method then consists of correctly reassigning colours according to the neighbouring cells. The whole diagram is constructed in a single phase with fault and formation datasets. The colours around faults are then assigned according to the maximum common surfaces with neighbouring formations (Fig. 6b). Though the ®rst method has been successfully applied in some two-dimensional cases (Fig. 3), it raises some problems in three dimensions for the following reasons: ² intersection problems are much more dif®cult to solve in three dimensions especially handling with tangency; ² when large displacements occur along faults, the initial solution (i.e. without faults) may be very far
Fig. 4. Splitting VoronoõÈ cells. Node 1 with neighbours 2±5 is split into two nodes with neighbours 2±4 and 3±5.
from what should be the actual solution, and there is then little chance to readjust properly the geological formations along the faults; ² the initial VoronoõÈ partition is lost due to intersections, consequently the model is no more suitable for insertion of new sites, and queries on point location cannot be performed on the model. The second method is therefore favoured for its algorithmic simplicity and for its more general application.
Fig. 5. Curvature of the mesh on a node Pi. An approximation of the curvature of the mesh is 1/R, R being the radius of the circumcircle of Pi, Pi1, Pi2 and Pi3. PiNi is the normal to the plane given by Pi1, Pi2 and Pi3. Position of Pi is moved along the normal PiNi according to smoothing criteria.
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
187
² geometric regularity constraint Ð smoothing must minimize some ®xed bending; ² topology constraint Ð smoothing must respect the initial topology of the reconstruction.
Vi 1; ¼; 4 in the graph of the mesh. Accounting for the geometry regularity constraint is performed locally on the current vertex and three neighbours by using the radius of the circumsphere as minimizing criteria. It is therefore necessary to split the nodes of degree 4 into two nodes of degree 3 (Fig. 4). All the vertices of the initial VoronoõÈ mesh are moved until an equilibrium is reached (the equilibrium criteria is that the total displacement of the vertices is less than a given value (e). After smoothing, the vertices of an initial polygonal face are no longer coplanar. They are therefore triangulated.
Since smoothing results in moving the vertices of the initial VoronoõÈ cells, care must be taken that data points de®ned inside geological formations remain inside after the smoothing process. From the VoronoõÈ diagram, the boundaries of geological formations are extracted as polygonal shapes. These can be seen as 3D meshes by considering their edges and vertices where each vertex V has at least three, and at most four, neighbouring vertices.
3.3.1. Geometric regularity constraint The curvature of a curve (surface) is a local property which quanti®es how much it bends. A geometrical de®nition can be given: the largest circle (sphere) that is tangential to a curve (surface) on its concave side at a point has the same curvature than the curve at that point. On each vertex we can de®ne a discrete approximation
3.3. Smoothing process Because of the roughness of the VoronoõÈ reconstruction (Fig. 3c), the shapes of the interfaces do not look quite ªnaturalº. It is therefore necessary to smooth the boundaries according to the following constraints:
Fig. 6. Taking faults into account. (a) VoronoõÈ sites with an indeterminate colour are generated on each side of the fault. (b) VoronoõÈ diagram and boundary extraction. The cells around the fault are assigned a colour according to the maximum common surfaces with neighbouring formations.
188
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
189
of the normal and the curvature of the mesh M. The neighbours Pi1, Pi2 and Pi3 de®ne a plane PLi and the normal vector Ni of PLi will be the normal of the mesh at Pi (Fig. 5). An approximation of the curvature of the mesh can be given by 1/R, R being the radius of the circumsphere of Pi, Pi1, Pi2 and Pi3. Indeed, the larger the radius, the closer the points to a planar distribution. It is possible to adapt the position of Pi with respect to a given curvature, and the position of the three neighbours. The smoothing method consists of obtaining a more regular distribution of curvatures along all points of the mesh. This can be achieved by interpolating a new curvature in Pi from the curvatures of the three neighbours and then to recalculate Pi the interpolated curvature and three neighbours. This is done on all vertices of the mesh one after the other in random order. This process can be iterated several times.
regions can be obtained (Fig. 1b). The surfaces bounding two volumes are described as sets of connected polygonal faces. A two-dimensional example of the reconstruction process is illustrated in Fig. 3. The regions are a ®rst approximation of the geological formations obtained when applying the following rule: a point p x; y; z belongs to a formation F if the site closest to p x; y; z belongs to F. This rule does not take into account the heterogeneous spatial distribution of the data. However, it is used to provide a solution that is, in most cases, topologically correct, i.e. homotopic to the actual geological model. In order to improve the quality of modelling, several constraints may be added in this reconstruction algorithm; for example, minimizing or maximizing the volume of a formation, preserving the thickness of some formations, or taking into account the anisotropy of the medium.
3.3.2. Topology constraint: Delaunay ball constraint This is the most particular constraint of the reconstruction method. It allows the topology of the different shapes of the model to be respected during the smoothing step. The Delaunay triangulation has the property that the circumcircle of every triangle contains no sites in its interior; such a circle (or disk) will be simply called an empty disk for short. Maintaining the displacement of the mesh edges and vertices in the union of these empty disks will guarantee that the initial topology is respected and that data points are inside their proper volumes.
4. 3D Armor model
3.3.3. Construction of homogeneous regions After construction of the VoronoõÈ diagram, each site pi [ P has one corresponding VoronoõÈ cell, so that each cell can be assigned to the geological formation (indexed by a colour) of its corresponding site. By merging adjacent VoronoõÈ cells that have the same colour, a partition of the space into colour-connected
4.1. Geological context The Cadomian orogenic belt of northern Brittany (France) is classically interpreted in terms of the evolution of an active margin ending with the collision of a continental margin, a marginal basin and a volcanic arc (Chantraine et al., 1998). Metamorphic, magmatic and structural evolution resulting from this collision is recorded from 650 to 540 Ma (Panafrican orogeny). Shortening resulted in the accretion of different units and was accommodated mainly by the combination of a southwest-directed wrench and strike-slip tectonism (Bale and Brun, 1989). The convergence kinematics and the way in which the different units are involved is still being debated (Egal et al., 1996; Strachan et al., 1996). Especially with regard to the geometry some questions arise: Do the structures ¯atten at depth? How are they branched? What is the relative importance of the
Fig. 7. Simpli®ed geological map of the main Cadomian units and geological cross-sections. Section 1: x1 217:643; y1 2395:084; x2 220:000; y2 2430:004: Section 2: x1 210:626; y1 2398:404; x2 226:505; y2 2423:857: Section 3: x1 203:488; y1 2408:585; x2 219:493; y2 2427:790: Section 4: x1 200:261; y1 2411:856; x2 216:931; y2 2437:662: Section 5: x1 201:749; y1 2410:827; x2 231:645; y2 2413:321: Section 6: x1 191:249; y1 2417:121; x2 231:168; y2 2414:590: Section 7: x1 230:706; y1 2383:129; x2 222:016; y2 2422:174: Section 8: x1 242:037; y1 2386:942; x2 219:850; y2 2420:225: Section 9: x1 258:209; y1 2386:664; x2 236:678; y2 2414:258: Coordinates are expressed in kilometres (kilometric Lambert2 coordinates system).
190
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
different units? The realization of a 3D geometric model serving as a reference model, allowing quantitative analysis, is a necessary starting point for answering these questions. Moreover, the contrast between physical and mechanical properties of the blocks makes an appropriate area for a testing methodology based on combined volumic and gravity modelling. 4.2. Lithology Units that have been modelled here are (Chantraine et al., 2001): (1) The Saint-Brieuc Unit composed of (1) the high metamorphic complex of the Yf®niac Formation, (2) Pentevrian orthogneisses including the Port-Morvan Formation, (3) dominantly basic metavolcanics of the Lanvollon Formation, (4) clastic metasediments of the Binic Formation overlying the metavolcanics, (5) basic plutonic intrusions (Squif®ec, St-Quay and Fort-Lalatte plutons). (2) The Saint-Malo±Guingamp Unit, originally composed of clastic brioverian sediments, characterized by the widespread development of migmatites. (3) Paleozoic granites intruding the Brioverian sediments in the southern part of the map. (4) Ordovician overstep sequence, composed of clastic sediments. Cadomian Units (Saint-Brieuc and Saint-Malo± Guingamp) are bounded by major faults along which they were accreted to the Central Armorican continental block.
4.3. Data Geometric modelling is based on a geological map and nine cross-sections oriented radially with respect to the arcuate shape of the belt (Fig. 7). The different units are discretized in terms of points belonging to each unit so that they can be suitable for VoronoõÈ reconstruction. In the particular case of the studied area, the major faults can be considered as geological interfaces since they form major boundaries between distinct units. Traces of geological contacts are transformed into points on each side of the contact. Assuming that sections are nearly perpendicular to the main structures, these points are generated within the section normal to the contact trace. They are then assigned their respective geological formation. The way contacts are discretized is illustrated in Fig. 8. Point density is globally equal to one point per 1 km along contacts but the distribution may vary according to the size and complexity of structures. Special care must be taken that point distribution re¯ects the characteristic features of the sections as, for example, triple points or thin narrow elongated bodies. Where a contact has an uncertain position or is poorly constrained, as is the case with most of the contacts at depth, the points are generated farther away from the drawn contact. 4.4. Results Starting from a given data set, the VoronoõÈ method allows one to automatically model the volume geometry of different geological units, while ensuring
Fig. 8. Example of discretization on Section 9. Points density is globally equal to one point per 1 km along contacts but their distribution may vary according to the size and complexity of structures. Care must be taken that points distribution re¯ects the characteristic features of the sections as for example triple points or thin narrow elongated bodies. Where a contact has an uncertain position or is poorly constrained, as is the case with most of the contacts at depth, the points are generated farther away from the drawn contact.
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
191
Fig. 9. 3D VoronoõÈ diagram of the Cadomian belt. (a) Volumic construction of VoronoõÈ diagram with no smoothing; (b) smoothing boundaries; and (c) extraction and visualization of different volumes from the model. The model is made up of different volumes that share common boundaries. 1 Ð Ordovician overstep sequence; 2 Ð Paleozoic granite intrusion; 3 Ð the Saint-Malo±Guingamp units; 4 Ð Pentevrian orthogneiss including the Port-Morvan Formation; 5 Ð undifferentiated Lanvollon and Binic formations.
192
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
topologic consistency between data and volumes. About 5 min on a current work-station are required to provide a model with 2000 data points. The principle of the method is illustrated through modelling the bulk architecture of the belt. In a ®rst step, ®ve volumes were constructed, corresponding to the Palaeozoic granites, the Saint-Malo±Guingamp Unit, Saint-Brieuc Unit where the volcano-sedimentary basin (Lanvollon±Binic) has been distinguished, and Ordovician sediments. Fig. 9a shows the resulting VoronoõÈ diagram without any smoothing. Though this representation ®ts the data, it remains geologically unrealistic. Fig. 9b shows the same volumes after smoothing, resulting in a more realistic model. Smoothing is controlled in order to respect the topological relations between data points and volumes. The arcuate shape of the belt is well illustrated. Major contact units can be extracted from the volumes, showing that the model is not only a set of
single volumes (Fig. 9c) but can be considered as a whole entity containing volumes sharing common boundaries. Fig. 10 shows a more detailed representation of the Saint-Brieuc Unit. The volumic imagery of these different formations makes the model comprehensive and communicable, helping to pose problems by highlighting inconsistencies. As example, southward dips of Port-Morvan Orthogneiss in Eastern part of the model brings some torsion in the shape of the model that suggests that Eastern and Western parts of the model are not continuously and simply correatable. Adding or correcting data in the model is easily done by simply inserting new sites in the diagram. As an example, Fig. 11 shows two models of the Binic Formation and the Ordovician overstep sequence. The ®rst model was poorly constrained in the north-eastern part while some fuzzy information on the estimated depth of the interface between the
Fig. 10. Modelling geological formations of the Cadomian belt. 1 Ð Ordovician overstep sequence; 2 Ð Paleozoic granite intrusion; 3a Ð the Saint-Malo±Guingamp units including 3b Ð the Quessoy orthogneiss; the following formations of the Saint-Brieuc Unit are represented: 4a Ð Pentevrian orthogneiss including the Port-Morvan Formation; 4b Ð basic intrusions including the Saint-Quay gabbrodiorite; 4c Ð the Yf®niac Formation; 4d Ð clastic sediments of the Binic Formation.
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
Fig. 11. Binic Formation and Ordovician overstep sequence. (a) The model for the north-eastern part was poorly constrained. Fuzzy information on the estimated depth of the interface has been added to the second model. 1 Ð Ordovician overstep sequence; 2 Ð clastic sediments of the Binic Formation.
two formations has been added in the second model, leading to a more geologically reasonable model. It is clear that the method only re¯ects the geometry of the underlying geological model contained in the cross-sections. However, some quanti®cation on the three-dimensional consequences of a hypothesis can be performed. For example, the simulation of geophysical effect can be computed based on the assumed physical properties of the medium. Since modelled volumes are closed and de®ne a whole partition of the space, the classical method of Okabe (1979) has been used to compute the gravity effect of the model (Fig. 12). When compared to the observed anomaly, discrepancies are analyzed in order to discuss the geological uncertainties and eventually to converge
193
to the best satisfying model. It can be seen that, compared to the observed one (Truffert et al., 2001), the ®rst-order variations are in agreement but that negative anomalies are overestimated. This is probably due to an insuf®cient precision in Hercynian granite de®nition. The depth and thickness of Port-Morvan orthogneiss are probably under-estimated. Important discrepancies exist at the north-east of the model owing to the lack of data. Great uncertainty exists with regard to the extension of Binic Serie and Lanvollon Serie in the North East under the Ordovician overstep sequence. This area clearly needs to be better constrained. Another reason is that the chosen formation regroupments may not be the best ones with regard to gravity properties. So, new sections calibrated both on geological constraints and on gravity anomaly pro®les have been redrawn. These were used to reconstruct the model and led to a more satisfying model. This latter model and the comparison between the computed gravity effect and observed anomaly are discussed in Truffert et al. (2001). However, a remaining problem concerns how and where to interact on the initial sections or on the model to reduce discrepancies. Some solutions are explored using the Monte-Carlo inverse method to randomly walk around the model space solution (Mosegaard and Tarantola, 1995). 5. Discussion The ®rst interesting aspect of the method is mainly its global character: because all volumes are built simultaneously, it is not necessary to construct each geological interface separately and independently to achieve a complete model, as with most surface methods. Taking into account data on the interior of the formations and not only on their contacts is a great advantage. The simplicity of the method allows the geologist to focus on data speci®cation and validity rather than on fastidious 3D modelling reconstruction methods. Only 3D points and their geological assignment are necessary to build a model. The second interesting aspect is the automatic character of the method. The current methodology for geological modelling consists of testing geological hypotheses and comparing their gravimetric or
194
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
Fig. 12. Gravity anomaly computed from the geometric model. Coordinates are expressed in kilometres (kilometric Lambert2 coordinates system).
magnetic effect against the recorded anomalies. It is clear that an automatic model greatly facilitates testing many hypotheses. The general idea that 3D modelling allows as many hypotheses as required to be tested is therefore given real meaning. However, several points still need to be improved for this idea to become a true reality. ² The ®rst point concerns the data management facilities used for building the models. Testing many hypotheses requires tools that allow easy modi®cation and editing of geological cross-sections and maps. Work is currently in progress on this aspect. ² The second point concerns the method itself. The present model is much too sensitive to the way data have been discretized; for units that are too ªthinº in places, as for example the Yf®niac Formation, the volume is discontinuous where it is expected to be
continuous. This is due to the algorithm that assigns the formation to the nearest neighbour. In order to improve the quality of the results, two approaches are currently being investigated: ± incorporating an anisotropic rule to allow the direction and strain of geological formations to be taken into account; ± improving the discretization methods in order to respect geological constraints (hypotheses on layer continuity, thickness, strain data, anisotropy). Another improvement of the approach will be an automatic modi®cation of the geometric model using the inverse problem to ®t geophysical measurements. This can be done easily because modifying cells or adding new sites in the VoronoõÈ diagram is rapid and so there is no need to rebuild the whole VoronoõÈ.
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
Consequently, there is no need to recompute the whole anomaly for each incremental modi®cation. 6. Conclusions The principal result of this work is the demonstration that volume reconstruction using VoronoõÈ diagrams is suitable for geological objects. The volume reconstructions are automatic, reproducible and give a consistent partition of space according to the data speci®cations. This can make the geological modelling process more ef®cient by providing the possibility of testing many hypotheses rapidly. In the case of the Cadomian belt of northern Brittany, the model that was built illustrates one of the possible hypotheses on the architecture of the belt. Other models based on different interpretations also have been tested using VoronoõÈ reconstruction in this issue (Truffert et al., 2001). Moreover, in the light of this work, 3D inversion algorithms for ®tting geophysical data onto the data structure of the VoronoõÈ diagrams can be engaged positively. Acknowledgements This work was funded by the GeÂofrance3D programme and by a BRGM research project. We are grateful to Jean Chantraine and Emmanuel Egal for their involvement in this work. We acknowledge Delphine Rouby and Jan-Diederik van Wees for their constructive review of the manuscript. GeoFrance 3D publication 19. References Auerbach, S., Schaeben, H., 1990. Computer-aided geometric design of geologic surfaces and bodies. Math. Geol. 22 (8), 957±987. Aurenhammer, F., 1991. VoronoõÈ diagrams: a survey of a fundamental geometric data structure. ACM Comput. Surv. 23, 345±405. BaleÂ, P., Brun, J.P., 1989. Late Precambrian thrust and wrench zones in Northern Brittany (France). J. Struct. Geol. 11 (4), 391±405. Bertrand, P., Dufourd, J.F., FrancËon, J., Lienhardt, P., 1992. ModeÂlisation volumique aÁ base topologique: actes MICAD 92 Paris 1, 59±74. Blanchin, R., Chiles, J.P., 1993. The Channel tunnel: geostatistical
195
prediction of the geological conditions and its validation by the reality. Math. Geol. 25 (7), 963±974. Boissonnat, J.D., 1988. Shape reconstruction from planar crosssections. Comput. Vis. Graph, Image Process. 44 (1), 1±29. Boissonnat, J.D., Nullans, S., 1996. Reconstruction of geological structures from heterogeneous and sparse data. Research report no. 3069, INRIA, France, 24 pp. Boissonnat, J.D., Teillaud, M., 1993. On the randomized construction of the Delaunay tree. Theor. Comput. Sci. 112, 339±354. Chantraine, J., Chauvel, J.J., BaleÂ, P., Denis, E., Rabu, D., 1998. Le BrioveÂrien (ProteÂrozoõÈque supeÂrieur aÁ terminal) et l'orogeÁne cadomien en Bretagne (France). Bull. Soc. GeÂol. Fr. IV, 815±829. Chantraine, J., Egal, E., ThieÂblemont, D., Le Goff, E., Guerrot, C., BalleÁvre, M., Guennoc, P., 2001. The Cadomian active margin (North Armorican Massif, France): a segment of the North Atlantic Panafrican Belt. Tectonophysics 331 (1±2), 1±18. Edelsbrunner, H., Tan, T.S., 1992. An upper bound for conforming Delaunay triangulations. Proc. 8th Annu. ACM Symp. Comput. Geom., pp. 53±62. Egal, E., Guerrot, C., Le Goff, E., ThieÂblemont, D., Chantraine, J., 1996. The Cadomian orogeny revisited in northern Brittany. Avalonian and Related Peri-Gondwana Terranes of the Circum-North Atlantic, Geol. Soc. Am. Spec. Publ. 304, 281±318. Farin, G., 1988. Curves and Surfaces for Computer Aided Design, A Practical Guide. Academic Press, San Diego 401 pp. Foley, J., van Dam, A., Feiner, S., Hughes, J., 1995. C Edition Ð Interactive Computer Graphics: Principles and Practice. Addison-Wesley, Reading, MA 1174 pp. Halbwachs, Y., Courrioux, G., Renaud, X., Repusseau, P., 1996. Topological and geometric characterisation of fault networks using 3-dimensional generalized maps. Math. Geol. 28, 625±656. Mallet, J.L., 1992. Discrete smooth interpolation in geometric modelling. Comput. Aided Des. 24, 178±191. Mallet, J.L., 1997. Discrete modeling for natural objects. Math. Geol. 29 (2), 199±219. Mayoraz, R., 1993. ModeÂlisation et visualisation infographiques tridimensionnelles de structures et proprieÂteÂs geÂologiques. TheÁse, Sciences, EPF Lausanne, No. 1127, XX, 218 pp. Meyers, D., Skinner, S., Sloan, K., 1991. Surfaces form contours: the correspondence and branching problems. Graphics Interface 91, 246±254. Mosegaard, K., Tarantola, A., 1995. Montecarlo sampling of solutions to inverse problems. J. Geophys. Res. 100, 12,431±12,447. Okabe, M., 1979. Analytical expressions for gravity anomalies due to homogeneous polyhedral bodies and translations into magnetic anomalies. Geophysics 44, 730±741. Renard, P., Courrioux, G., 1994. 3-D geometric modeling of a faulted domain: the Soultz horst example (Alsace, France). Comput. Geosci. 20, 1379±1390. Seidel, R., 1994. The nature and meaning of perturbations in geometric computing. Proc. 11th Symp. Theor. Aspects Comput. Sci. (STACS), Lecture Notes in Computer Science, vol. 775. Springer, Berlin, pp. 3±17.
196
G. Courrioux et al. / Tectonophysics 331 (2001) 181±196
Siehl, A., Ruber, O., Valdivia-Machego, M., Klaff, J., 1992. Geological maps derived from interactive spatial modeling. Geol. Jahrb. 122, 273±279. Strachan, R.A., D'Lemos, R.S., Dallmeyer, R.D., 1996. Neoproterozoic evolution of an active plate margin: North Armorican Massif, France. Avalonian and Related Peri-Gondwana
Terranes of the Circum-North Atlantic, Geol. Soc. Am. Spec. Publ. 304, 319±332. Truffert, C., Egal, E., Le Goff, E., Courrioux, G., Guennoc, P., 2001. Gravity modelling of the Cadomian active margin of northern Brittany. Tectonophysics 331 (1±2), 81±97.