Liangfeng Zhu*, Chengjuan Zhang, Mingjiang Li, Xin Pan, Jianzhong Sun
Key Laboratory of GIS Science for Ministry of Education, East China Normal University, Shanghai 200062, China
* Corresponding author. E-mail: zhuliangfeng@163.com. Tel.: +86 21 62237221. Fax: +86 21 62232332.
Postal address: Department of Geography, Key Laboratory of GIS Science for Ministry of Education, East China Normal University, No.3663, N. Zhongshan Rd., Shanghai, P. R. China 200062
E-mail addresses: zhuliangfeng@163.com (L. Zhu), 51090801048@ecnu.cn (C. Zhang), lmj_fighting@sina.com (M. Li), xpan@admin.ecnu.edu.cn (X. Pan), sunjzh2007@gmail.com (J. Sun).
Abstract
3D solid models of geological structures are particularly useful to practical geological analysis and engineering design. The main difficulty raised by 3D geological modeling of sedimentary system is determining the geological genesis and geometrical boundaries of missing strata. For the lack of the comprehensive mechanism to handle missing strata, it is difficult to construct spatial geometric shapes of complicated strata with a desired accuracy in 3D utilizing the existing modeling methods. This situation limits the reliability and the practicality of the computer models. In order to construct the discontinuous geological surfaces induced by missing strata, an adapted and automatic approach for generating 3D solid models of sedimentary stratigraphic systems from borehole data, called the Boreholes-Surfaces-Solids method, is presented. The method first utilizes the topologic dimidiate data structure to discretize borehole data into a series of scatter points, then interpolates the initial elevations of the top and bottom surfaces for each stratum, and automatically deduces the genesis of the missing strata. Subsequently, according to different geological genesis, surfaces intersecting, elevations adjusting and consistency processing are performed automatically on the missing strata’s surfaces and their control surfaces. And finally, the solid model filled with 3D blocks or triangular prism meshes is built. The Boreholes-Surfaces-Solids method has higher automaticity and stronger adaptability, and overcomes limitations of the existing modeling methods. Two concrete examples of using this method to Shanghai’s construction projects show that the resulting models are natural, geologically reasonable and close to the actual stratigraphic distribution. In addition, during the implementing process, the geological laws, such as the different genesis controlling the spatial geometry of the missing strata, are converted skillfully into modeling rules that can be identified and programmed automatically by modelers, and this is helpful to promote further study on 3D modeling techniques for complex geological structures.
Key words
sedimentary stratigraphic system; borehole; 3D solid modeling; missing stratum; Boreholes-Surfaces-Solids method
1. Introduction
Computer modeling and visualization of geological objects in 3D is currently a topical research area both in Engineering Geology and Geo-information Science (Jones, 1988; Hack et al., 2006; Turner, 2006; De Rienzo et al., 2008; Royse et al., 2009). 3D solid models constructed utilizing 3D geoscience modeling techniques are particularly useful to practical geological analysis and engineering design. Solid models of geological objects in 3D can vividly define the boundaries of different geological phenomena and complex structures within geological units, and then enhance the visibility and accuracy of geological analysis and engineering design. Nowadays, 3D solid models are widely used in a number of fields such as the geometrical representation of geological structures, the visual analysis of spatial inhomogeneity for geological properties, the preprocessing and postprocessing of numerical simulation models, etc. The research objects of 3D geoscience modeling are geological objects buried in the subsurface of the Earth’s crust. In terms of their morphological characteristics, geological objects can be classified as “stratified” or “non-stratified” objects. The sedimentary stratigraphic system, which is composed of stratified objects, and common in earth surface system, not only contains energy, minerals, groundwater and other resources, but also provides vast areas of fertile farmland and construction sites. At present, many large cities and engineering facilities are built in the delta areas covered by sedimentary strata. Therefore, building solid models of sedimentary stratigraphic system in 3D is extremely important. In addition, the research of 3D modeling techniques for a sedimentary system still has great theoretical significance since it is a foundation for further study on more complicated geological structures.
Over the past three decades, a series of modeling theories and techniques have been presented to address the needs of 3D solid models for sedimentary stratigraphic system. The related research involves three essential aspects: standardization of modeling data, representation of 3D solids, and construction techniques of 3D geological objects. In the realm of the standardization of modeling data, the past research has focused on the design of data standards and code systems for geological data. For example, a typical standard form of borehole data was suggested and implemented in a web-based GIS system (Chang and Park, 2004). In addition, considerable effort has been given to develop robust management systems and processing procedures for various types of data, especially boreholes, cross-sections and contours (Nathanail and Rosenbaum, 1998; McCarthy and Graniero, 2006). Several software systems, such as Geotouch (Lees, 2000) and BoreIS (McCarthy and Graniero, 2006), were developed as tools to aid in the storage, manipulation, visualization, querying and analysis of borehole and other geological data. These essential advances provide tremendous supports for 3D geological modeling and visualization. In the realm of the representation of 3D solids, the past research has focused on 3D spatial data models and their data structures that are especially suitable for stratified geological objects. More than twenty data models were proposed and described for 3D geological modeling such as vector octree (Jones, 1989), grid, triangulated irregular network (TIN), tri-prism, generalized tri-prism (GTP) (Wu, 2004), tetrahedral network (TEN), supervoxel (Wu and Xu, 2004), Geocellular, NURBS-TIN-BRep hybrid (Zhong et al., 2006) and so on. In the realm of the construction techniques of 3D geological objects, the past research has focused on generating 3D solid models from different types of data like boreholes, cross-sections or geological maps. Several approaches to build 3D solids have been developed and applied (He et al., 2002; Lemon and Jones, 2003; Zhu et al., 2004; He et al., 2005; Zhu and Pan, 2005; Zhang et al., 2006; Smirnoff et al., 2008; Tremblay et al., 2010). However, up to now there are still no perfect methods or easy-to-handle software toolkits for the reconstruction of sedimentary stratigraphic systems. In brief, for the standardization of modeling data and the representation of 3D solids, all or most issues have been solved, common understandings and recognitions have been reached in academic and engineering fields; while for construction techniques of 3D geological objects, a comprehensive modeling approach, which can handle all types of geological objects simply, efficiently and automatically, is still lacking. Currently, it is an essential task to optimize, improve and validate the existing modeling approaches coupling with various types of geological settings.
Engineering drilling is a traditional technique to observe and sample the subsurface directly. Borehole data are simple, intuitive, exact and detailed for practical users. Geologists, engineering geologists and geotechnical engineers are all quite familiar with borehole data, and a variety of critical information like stratigraphical observational data can be obtained very easily from a normative borehole database (He et al., 2002; Turner, 2006). While several techniques have been developed for building 3D geological models, the most popular technique is still based on boreholes. A number of research teams have invested considerable effort in developing modeling methods from boreholes directly, and several methods to the construction, modeling and representation of sedimentary stratigraphic system were proposed and applied (He et al., 2002; Lemon and Jones, 2003; Zhu et al., 2004; He et al., 2005; Zhu and Pan, 2005; Zhang et al., 2006; Smirnoff et al., 2008; Tremblay et al., 2010). These methods mostly differ in the details of their interpolation of data and in their ability to represent complex structures in 3D (Wellmann et al., 2010).
Several experiments have highlighted a number of shortcomings and some serious limitations when using these existing methods (Zhu and Pan, 2005; Zhang et al., 2006). A critical problem is that the modeling result differs from the actual structure in geometric shapes. Thus, the computer model is geologically unreasonable in some situations, and the result model can not reflect the actual spatial distribution characteristics of the complex geological object. The reason why this happens is a comprehensive treatment of all types of missing strata in sedimentary system is still lacking when utilizing these existing methods. The discontinuous spatial distribution, frequently induced by missing strata, increases the complexity of the interpolating and fitting process for geological interfaces. In order to precisely control the shapes of missing strata in the resulting models, some research teams proposed that the additional data, like user-defined cross-sections (Lemon and Jones, 2003) or virtual boreholes (Zhu et al., 2006), can be merged into the modeling procedure. However, drawing numerous cross-sections or virtual boreholes manually in a complicated area requires significant user intervention, and mainly depends on the modeler’s judgments with great subjectivity, thus it is extremely tedious and time-consuming. Therefore, this method is difficult to use and not very well suited to automated modeling processes with computers.
This paper explores a new 3D solid modeling method which is directly from borehole data, suitable for sedimentary stratigraphic systems, and taking into account the influence of missing strata. The rest of the paper is organized as follows. Structure characteristics of sedimentary stratigraphic system are summarized in Section 2, which also considers the classification of missing stratum. Section 3 introduces the modeling methodology of Multi-layer DEM, and a novel implementation method, called the Boreholes-Surfaces-Solids method, is presented. Section 4 concentrates on the major steps and technical details of the Boreholes-Surfaces-Solids method. We demonstrate the application of the Boreholes-Surfaces-Solids method to Shanghai’s construction projects in Section 5. The conclusions of this paper are provided in Section 6.
2. Structure characteristics of sedimentary stratigraphic system
In recent sedimentary systems, the geometric shape of each stratigraphic unit is relatively simple and regular since it is not intersected with any faults, joints or other fracture structures. Strata in sedimentary system can be classified as “complete” or “missing” strata in terms of their integrities and spatial distribution characteristics (Zhang et al., 2006). A complete stratum is one that distributes continuously in a given study area, thus its top or bottom interface is a complete curved surface without any “void hole”. In contrast, a missing stratum is one that distributes discontinuously in a given study area, thus its top or bottom interface is composed of either a continuous curved surface with one or more “void holes”, or a combination of multiple disconnected continuous curved surfaces.
In this paper, the missing strata in sedimentary systems are also classified into three different types: Type 1, Type 2, or Type 3, depending on their stratigraphic settings and structure characteristics.
l Type 1: This type of missing stratum is commonly induced by non-deposition. That is to say, in a given area, there is no sediment all the time. In general, the missing stratum of Type 1 derives from the endogenic geological process of the earth, thus the geometrical boundaries of missing stratum are controlled by the top surface of the underlying stratum.
l Type 2: This type of missing stratum is commonly induced by erosion. That is, in a given area, the historical sedimentation had ever occurred, but it was eroded completely before the deposition of the overlying stratum. The missing stratum of Type 2 generally derives from exogenic geological processes, such as fluvial erosion and weathering denudation, thus its geometrical boundaries are controlled by the bottom surface of the overlying stratum.
l Type 3: This type of missing stratum is induced by the superposition, and/or the compound, of non-deposition and erosion. The geometrical boundaries of Type 3 missing stratum are controlled by both the bottom surface of the overlying stratum and the top surface of the underlying stratum.
A typical sedimentary system with a set of stratigraphic units is illustrated in Fig.1 (Turner, 2006). In this case, S1, the lowest stratum, is a complete stratum; S2, the middle stratum, is a missing stratum of Type 2 since its missing areas are induced by erosion, and its geometrical boundaries are controlled by the bottom surface of the overlying stratum; and S3, the topmost stratum, is a missing stratum of Type 3 since its missing areas are partly induced by non-deposition and partly induced by erosion.
(Fig.1. A typical sedimentary system with missing strata)
Based on large numbers of case studies, we summarize five major structure characteristics of sedimentary stratigraphic systems. These characteristics are:
l Stratified. In terms of a given criteria for classification, the sedimentary system can be divided into several stratigraphic units, and each unit is called as a stratum. Within every stratigraphic unit, the depositional age and the mechanical properties are assumed to be approximately uniform, thus each stratum can be regarded as being composed of the same soil or rock mass, and commonly denoted as a “geotechnical unit”.
l Sequential. For a given study area, a standard depositional sequence can be established in terms of given criteria or rules.
l Continuous. In a sedimentary system, the complete stratum distributes continuously while the missing stratum distributes discontinuously in local areas. Nevertheless, in the missing area of a missing stratum, the bottom interface can be regarded as coinciding with the top interface. Thus, the missing stratum can be treated specially as the complete stratum with “zero-thickness” units, and all strata in sedimentary system can be regarded as continuous stratified geological objects (Xu and He, 2004).
l Enclosed. Each stratum is enclosed by the top, the bottom and the side surface. The top and bottom surfaces can be imagined as two curved surfaces that can be completely projected onto the reference horizontal plane with coincident forms. During the modeling process, the modeler just needs to reconstruct the top and the bottom surface, while the side surface can be generated automatically.
l Single-valued. On the top or bottom surface of a certain stratum, a unique elevation value corresponds to a given planar point. Thus, the top or the bottom surface of each stratum can be assumed to be single-valued with respect to a 2D coordinate system.
The modeling result of sedimentary system is known as the “layer-cake” model (Turner, 2006). But in fact, the actual strata are more complicated than the layers of a cake. Missing strata and discontinuous surfaces are very common since stratigraphic interfaces may intersect each other in complicated areas. It is not easy to effectively, accurately construct a geologically reasonable model directly from boreholes. During the modeling process, the modeler must consider comprehensively the influence of missing strata in order to provide a means for dealing with special geological structures such as pinchouts, intrusions and lens.
3. Modeling methodology of Multi-layer DEM
The dominant characteristic of a sedimentary system is the sequential, regular stacking of sedimentary strata and their interfaces (Turner, 2006). Although methodologies for the description and modeling of a sedimentary system in 3D have been developed and explored for many years, the most important step of these techniques is still to define and simulate the top and bottom surfaces of each stratum. In recent years, Multi-layer DEM technique has been investigated by several research teams as a practical method to construct 3D stratigraphic models for sedimentary systems (He at al., 2002). Four steps are generally followed in this technique: first, according to the stratified information contained in geological boreholes, a clear and well-organized depositional sequence of all stratigraphic units in the study area should be obtained and determined; and then, based on the control sample points for each interface, a series of DEM surfaces can be interpolated independently utilizing a 2.5D method; after that, the modeler need to conduct the intersection and division operation of multiple DEM surfaces, and the spatial geometric framework of geological objects is formed in terms of the properties of strata; and finally, according to the fundamental geological framework, each stratigraphic unit is subdivided into a serious of structured or unstructured meshes by applying discretization methods, and a 3D voxel-based solid model incorporating the geometric framework information of each stratum is established.
The advantage of the Multi-layer DEM modeling technique is obvious. It requires simple input data, and the modeling process is fast and straightforward. Currently, there are several implementation algorithms for the Multi-layer DEM technique. Some algorithms that have been used more often are horizons-to-solids algorithm (Lemon and Jones, 2003), boreholes-interfaces algorithm (Zhu et al., 2004), strata-framework algorithm (Zhu and Pan, 2005), and vertical sub-block algorithm (Zhang et al., 2006).
To meet the requirements of 3D geological models in the construction projects of Shanghai, China, we have used these different algorithms for creating 3D solid models from boreholes. The reliability of these algorithms is tested with actual data, and some drawbacks are detected. The chief disadvantage of these algorithms is their difficulty in handling the missing strata of sedimentary system. For example, horizons-to-solids algorithm (Lemon and Jones, 2003), boreholes-interfaces algorithm (Zhu et al., 2004), and vertical sub-block algorithm (Zhang et al., 2006) are ideally suited for alluvial systems, as these just successfully deal with Type 1 missing strata induced by non-deposition. Strata-framework algorithm (Zhu and Pan, 2005) can handle Type 2 missing strata induced by erosion, but it is inappropriate for Type 1 missing strata. Type 3 missing strata are more difficult to handle. It is impossible to construct missing strata of Type 3 only using the existing algorithms and borehole data. In order to construct geologically reasonable missing strata in complicated areas, large numbers of cross-sections or boreholes need to be added as additional control data (Lemon and Jones, 2003; Zhu et al., 2006). Up to date, there is no algorithm which is based on Multi-layer DEM technique to comprehensively handle all types of missing strata in 3D geological modeling.
Whereas the existing algorithms are very effective for complete strata with clear depositional sequences and approximately horizontal spatial distributions, they are not suited for complicated sites which have strong geological tectonic activities and multiple types of missing strata. We believe this always applies due to the lack of a comprehensive consideration of geological settings for all types of missing strata in existing modeling algorithms. These algorithms fail to convert the geological laws into the modeling rules that can be identified and programmed automatically by computers. In order to produce geologically sound models, the existing modeling algorithms need to be improved by coupling with various types of missing strata and geological settings.
Based on Multi-layer DEM technique, a novel method, called the Boreholes-Surfaces-Solids method, is presented for modeling sedimentary systems in 3D, which not only effectively handles the missing strata of Type 1 and Type 2, but also automatically handles the missing strata of Type 3 simultaneously. The Boreholes-Surfaces-Solids method first utilizes the topologic dimidiate data structure (Zhu and Wu, 2005; Zhang et al., 2006) to discretize borehole data into a series of scatter points, then interpolates the initial elevations of the top and bottom surfaces for each stratum, and automatically deduces the genesis of missing strata. Subsequently, according to different geological genesis, surfaces intersecting, elevations adjusting and consistency processing on the missing strata’s surfaces and their control surfaces are performed automatically. And finally, the solid model filled with 3D blocks or triangular prism meshes is built. The Boreholes-Surfaces-Solids method overcomes limitations of the existing modeling methods, and enables the fast construction of a geologically reasonable model directly just using borehole data. This method has a higher level of automatic process and stronger adaptability.
4. Modeling steps involved in the Boreholes-Surfaces-Solids method
The Boreholes-Surfaces-Solids method is based on the recent geological modeling methods (Lemon and Jones, 2003; Zhu et al., 2004; Zhu and Pan, 2005; Zhang et al., 2006), but differentiates itself from the previous methods since we consider geological settings of the missing strata during the modeling process. The Boreholes-Surfaces-Solids method can be performed using an automatic, straightforward fashion that makes the shapes of missing strata to be effectively controlled just with borehole data. As Fig.2 shows, this method involves 11 steps, and the following is an explanation of the main steps.
(Fig.2. Modeling flow of the Boreholes-Surfaces-Solids method)
For a given study area, maybe numerous boreholes have been collected before 3D modeling process, and all these borehole data can be organized and stored into a GIS system based on a standard database format (Chang and Park, 2004). The first step in the modeling process is to extract borehole data of the site being modeled. A variety of critical information, such as the types of the boreholes, the locations of the boreholes in 2D space and the detailed stratified data, is extracted from the GIS database and can be used as the initial sample data for the subsequent processes.
First, each stratigraphic unit intersected with boreholes is assigned a stratum identifier in terms of the depositional sequence, and an ordered strata list that contained all strata in the site being modeled is created. The stratum identifier represents the order in the bottom-to-top depositional sequence. The stratum identifiers should start at 1 and increase from the bottom to the top. Therefore, the bottommost (also the oldest) stratigraphic unit is denoted as S1 with a stratum identifier = 1, and the ith (i ≥ 1) stratigraphic unit is denoted as Si with a stratum identifier = i.
And then, boreholes that extracted from the database are divided into “complete” and “partial” boreholes in terms of their integrity. If a borehole detects both the bottommost and the topmost stratum, it is a complete borehole; otherwise, it is a partial borehole.
And finally, strata of the site being modeled are broadly separated into two categories: complete strata and missing strata. If a stratum is detected by all complete boreholes as well as the partial boreholes which have chances to meet this stratum, it is a complete stratum; otherwise, it is a missing stratum.
It should be pointed out that the topmost and the bottommost stratigraphic units of the site being modeled are assumed as complete strata for the convenience of the subsequent processes. If the site being modeled cannot fulfill this requirement, a virtual complete stratum will be added automatically above the topmost stratum and/or below the bottommost stratum. And just after “Step 10: Build 3D solid model” (shown in Section 4.10), the additional virtual stratum will be removed automatically.
In this paper, borehole data are organized into contacts. A contact is defined as the interface between two adjacent strata (Lemon and Jones, 2003). Each contact has a borehole identifier, a location (x, y, z), an identifier for the stratum above the contact, and an identifier for the stratum below the contact. In the third step of the modeling process, contacts of boreholes are discretized into a series of scatter points.
We apply the topologic dimidiate data structure (Zhu and Wu, 2005; Zhang et al., 2006) to organize the stratified data of boreholes. In this data structure, each contact is discretized into a scatter point, and the topological properties of the scatter point are described by two identifiers of the contact (one for the stratum above and one for the stratum below). The data structure of the scatter point is described as follows:
Struct BoreholeContactPoint{
long m_lBoreholeID; // identifier of borehole
double m_dX; // x coordinate of the contact location
double m_dY; // y coordinate of the contact location
double m_dZ; // z coordinate of the contact location
int m_iAboveID; // identifier for the stratum above the contact
int m_iBelowID; // identifier for the stratum below the contact
}
A sample of discretization of borehole data is illustrated in Fig.3. A sedimentary system with four stratigraphic units and eight boreholes is shown in Fig.3A, and the scatter points resulting from the boreholes are shown in Fig.3B.
(Fig.3. Discretization of borehole data based on topologic dimidiate data structure)
After discretization, contacts of all boreholes are merged into one scatter point set P, which collects sample data for subsequent interpolating elevations of the top and bottom surfaces for each stratum. In addition, we must point out that not only boreholes but also other available geological data, such as cross-sections, contours and so on, can be discretized and merged into the scatter point set P as sample data for subsequent processes.
In sedimentary system, the horizontal projections of the top and bottom surfaces of each stratum is coincident with each other, then we can define a “primary TIN” as the reference triangular network to fit geometrical forms of each interface. The primary TIN is defined as a triangulated irregular network that based on the horizontal coordinates (x, y) of all boreholes, generated by constrained Delaunay triangulation with outer boundaries of the site as constrained-edges, and densified automatically by the subdivision operation (Lemon and Jones, 2003; Zhu et al., 2004). The primary TIN not only explicitly defines the outer boundary of the 3D solids, but also implicitly establishes the common topological and geometric relationship between the top and the bottom surface of each stratum. Using the primary TIN with a consistent topology for the stratum is a key to simplify the subsequent processes (Lemon and Jones, 2003), and it also effectively improves the robustness of the Boreholes-Surfaces-Solids method.
Each stratum is enclosed by the top TIN and the bottom TIN since the side surface can be generated automatically. In the fifth step of the modeling process, we extract 3D coordinates (x, y, z) of the sample data from the scatter point set P, and interpolate the initial elevations of each vertex in each top or bottom TIN for each stratum.
When Fi_Top, the top TIN of stratum Si that numbered as i, is being interpolated, we exact all scatter points with m_iBelowID = i from P as sample data; and when Fi_Bottom, the bottom TIN of Si, is being interpolated, we exact all scatter points with m_iAboveID = i from P as sample data.
Several often-used interpolation schemes, like the inverse distance weighted (IDW), natural neighbor, the nearest neighbor distance, radial basis function (RBF), and Kriging methods, can be used to interpolate the elevations. These methods are relatively simple, convenient and robust since they all support both interpolation and extrapolation, and produce excellent results (Lemon and Jones, 2003).
The sixth step is to estimate and deduce automatically the types of all missing strata in terms of the initial elevations of the top and bottom surfaces for each stratum. For a given missing stratum Si, which has an overlying stratum Si+1 and an underlying stratum Si-1, the top TIN of Si is denoted as Fi_Top, and the bottom TIN of Si is denoted as Fi_Bottom. If Si is absent in the location of borehole Bj, Bj is called as a missing borehole for Si. In the location of Bj, the contact between Si+1 and Si-1 is denoted as P0, the elevation of P0 is denoted as Z0, the vertex of Fi_Top is denoted as PTop, and the vertex of Fi_Bottom is denoted as PBottom. The initial elevations of PTop and PBottom, denoted as ZTop and ZBottom respectively, are interpolated by Step 5. The distance between PTop and P0 is denoted as d1, and the distance between PBottom and P0 is denoted as d2. The values of d1 and d2 can directly be calculated as follows:
d1=ZTop - Z0
d2=ZBottom - Z0
Since ZTop ≥ ZBottom, then d1 ≥d2. For each pair of d1 and d2, there are four cases to be considered:
l Case 1: d1=0, and d2=0. In this case, PTop and PBottom happen to coincident with P0. Therefore, in the location of borehole Bj, Si is pinch out and can be regarded as the continuous stratum with “zero-thickness” unit.
l Case 2: d1 ≤ 0, and d2<0. In this case, we regard Si is missing induced by non-deposition in the location of borehole Bj (as shown in Fig.4A).
l Case 3: d1>0, and d2 ≥ 0. In this case, we regard Si is missing induced by erosion in the location of borehole Bj (as shown in Fig.4B).
l Case 4: d1>0, and d2<0. In this case, we regard Si is missing induced by the superposition and compound of non-deposition and erosion in the location of borehole Bj (as shown in Fig.4C).
(Fig.4. Determine the types of missing strata)
Case 1 is not considered in the following since it has no impact on the missing stratum. If Case 2 occurs in every missing borehole for Si, we regard Si as a missing stratum of Type 1. If Case 3 occurs in every missing borehole for Si, we regard Si as a missing stratum of Type 2. If Case 2 and Case 3 occur together in missing boreholes for Si, or Case 4 occurs in any missing borehole for Si, we regard Si as a missing stratum of Type 3.
After the interpolating process performed in Section 4.5, the initial elevation of each vertex in the top and bottom TINs of each stratum is calculated. For a given complete stratum, its top TIN and bottom TIN cannot intersect with each other, and they also cannot intersect with the top or bottom TIN of other complete strata. However, for a missing stratum, its top TIN or bottom TIN may intersect with other TINs for the overlying or underlying strata (this depends on the type of missing stratum). Thus, we need to intersect the TIN surfaces and adjust the elevations for the missing strata and their control strata in terms of the types of missing strata.
If there is more than one missing stratum in the site being modeled, we need to define the priority of these missing strata to generate a processing sequence for intersecting TIN surfaces and adjusting elevations. In the seventh step of the modeling process, we arrange the orders of the missing strata in the light of the following rules:
(1) The missing stratum of Type 3 can be treated as the combination of a Type 1 and a Type 2 missing stratum. Consequently, if there are Type 3 missing strata in the site being modeled, we can decompose them into Type 1 and Type 2 missing strata. Therefore, a Type 3 missing stratum need to be recorded twice since it is both a Type 1 and a Type 2 missing stratum.
(2) The missing stratum of Type 1 has precedence over Type 2. Thus, we should first deal with all Type 1 missing strata, then deal with Type 2 missing strata.
(3) The older missing stratum of Type 1 has precedence over other relatively new missing stratum of Type 1. Therefore, if there is more than one Type 1 missing stratum in the site being modeled, we should sort them in the bottom-to-top sequence.
(4) The relatively new missing stratum of Type 2 has precedence over other older missing stratum of Type 2. Therefore, if there is more than one Type 2 missing stratum in the site being modeled, we should sort them in the top-to-bottom sequence.
After sorting process, we obtain a processing sequence for missing stratum which only contains Type 1 and Type 2 missing stratum.
The eighth step is to intersect the top and bottom TINs of each missing stratum with their control surface, and adjust elevations of the vertices on the top and bottom TINs of each missing stratum. Starting with the first missing stratum in the processing sequence generated by Step 7, the top and bottom TINs of each missing stratum are sequentially intersected with their control TIN and adjusted subsequently. Since the missing strata of Type 1 and Type 2 are controlled by different geological interfaces respectively, we need to apply different algorithms to deal with different types of missing strata. In addition, two essential characteristics to accelerate the intersecting process are also necessary to be described in detail.
If Si is a missing stratum of Type 1, its geometrical boundaries are controlled by the top surface of Si-1, the underlying stratum of Si. After the interpolating process performed in Step 5, the initial forms of Fi_Top, Fi_Bottom, F(i-1)_Top and F(i+1)_Bottom are shown in Fig.5A. Here Fi_Top is the top TIN of Si, Fi_Bottom is the bottom TIN of Si, F(i-1)_Top is the top TIN of Si-1, F(i+1)_Bottom is the bottom TIN of Si+1.
First is to handle Fi_Top, the top TIN of Si. Fi_Top is intersected with F(i-1)_Top, and all intersection points are calculated automatically. These intersection points are inserted into the primary TIN as new vertices. At the same time, we modify the vertices and the topology of Fi_Top,F(i-1)_Top and other TINs for each stratum. Normally, the elevation of a given vertex in Fi_Top cannot go below the elevation of the corresponding vertex in F(i-1)_Top. However, in the areas where Si is absent, the elevation of a given vertex in Fi_Top may go below the elevation of the corresponding vertex in F(i-1)_Top. In this case, the elevation of the vertex in Fi_Top is uplifted in order to set it equal to the elevation of the corresponding vertex in F(i-1)_Top (as the red arrows shown in Fig. 5A). The modified Fi_Top is denoted as F’i_Top.
Next is to handle Fi_Bottom, the bottom TIN of Si. Theoretically, the elevation of each vertex in Fi_Bottom should be equal to the elevation of corresponding vertex in F(i-1)_Top. However, since Si is a missing stratum, the sample points used to interpolate the elevations of the vertices in Fi_Bottom and F(i-1)_Top are not entirely identical. Thus, it leads to the elevations of corresponding vertices in these two TINs may not equal. In order to keep the consistency for these two surfaces, we need to intersect them and adjust elevations. Fi_Bottom is intersected with F(i-1)_Top, and all intersection points are calculated automatically. These intersection points are inserted into the primary TIN as new vertices, the vertices and the topology of Fi_Bottom,F(i-1)_Top and other TINs for each stratum are also modified. Then we compare the elevation of each vertex in Fi_Bottom with the elevation of corresponding vertex in F(i-1)_Top. If the elevation of Fi_Bottom is below (or above) F(i-1)_Top, we uplift (or depress) the elevation of Fi_Bottom and set it equal to F(i-1)_Top (as the blue arrows shown in Fig. 5A). The modified Fi_Bottom is denoted as F’i_Bottom.
After intersecting surfaces and adjusting elevations, the modified surfaces of Si are illustrated in Fig.5B.
(Fig.5. Intersect surfaces, adjust elevations and keep consistency for Type 1 missing stratum)
If Si is a missing stratum of Type 2, its geometrical boundaries are controlled by the bottom surface of Si+1, the overlying stratum of Si. After the interpolating process performed in Step 5, the initial forms of Fi_Top, Fi_Bottom, F(i-1)_Top and F(i+1)_Bottom are shown in Fig.6A. Here Fi_Top is the top TIN of Si, Fi_Bottom is the bottom TIN of Si, F(i-1)_Top is the top TIN of Si-1, F(i+1)_Bottom is the bottom TIN of Si+1.
First is to handle Fi_Bottom, the bottom TIN of Si. Fi_Bottom is intersected with F(i+1)_Bottom, and all intersection points are calculated automatically. These intersection points are inserted into the primary TIN as new vertices. Meanwhile, we modify the vertices and the topology of Fi_Bottom, F(i+1)_Bottom and other TINs for each stratum. Normally, the elevation of a given vertex in Fi_Bottom cannot go above the elevation of the corresponding vertex in F(i+1)_Bottom. However, in the areas where Si is absent, the elevation of a given vertex in Fi_Bottom may go above the elevation of the corresponding vertex in F(i+1)_Bottom. In this case, the elevation of the vertex in Fi_Bottom is depressed (as the red arrows shown in Fig. 6A), and it is set equal to the elevation of the corresponding vertex in F(i+1)_Bottom. The modified Fi_Bottom is denoted as F’i_Bottom.
Next is to handle Fi_Top, the top TIN of Si. Theoretically, the elevation of each vertex in Fi_Top should be equal to the elevation of corresponding vertex in F(i+1)_Bottom. However, since Si is a missing stratum, the sample points used to interpolate the elevations of the vertices in Fi_Top and F(i+1)_Bottom are not entirely identical. Thus, it leads to the elevations of corresponding vertices in these two TINs may not equal. In order to keep the consistency for these two surfaces, we need to intersect them and adjust elevations. Fi_Top is intersected with F(i+1)_Bottom, and all intersection points are calculated automatically. These intersection points are inserted into the primary TIN as new vertices, the vertices and the topology of Fi_Top, F(i+1)_Bottom and other TINs for each stratum are also modified. Then we compare the elevation of each vertex in Fi_Top with the elevation of corresponding vertex in F(i+1)_Bottom. If the elevation of Fi_Top is below (or above) F(i+1)_Bottom, we uplift (or depress) the elevation of Fi_Top and set it equal to F(i+1)_Bottom (as the blue arrows shown in Fig. 6A). The modified Fi_Top is denoted as F’i_Top.
After intersecting surfaces and adjusting elevations, the modified surfaces of Si are illustrated in Fig.6B.
(Fig.6. Intersect surfaces, adjust elevations and keep consistency for Type 2 missing stratum)
In general, intersecting two TINs is a complex and time-consuming process since each triangle of one TIN must be checked against each triangle of the other TIN (Lemon and Jones, 2003). However, since the top and bottom TINs of each stratum are all based on the primary TIN, they are identical in plan view with the same topology. Thus, a given triangle from the first TIN can only intersect the corresponding triangle from the second TIN. Additionally, when a TIN is intersected with another TIN, the new vertices are generated at the intersected locations and inserted into the primary TIN. We also need to calculate elevations of these new vertices. Since the new vertex surely lies on the edge of a TIN triangle, a simple linear interpolation can be used to compute the elevation of each new vertex (Lemon and Jones, 2003). During the modeling process, these two characteristics can be used to accelerate the intersecting process, and the performance of the Boreholes-Surfaces-Solids method can be improved greatly.
In this step, we are ready to handle the interfaces between the missing strata and their adjoining strata in order to keep the consistency of all surfaces. Starting with the first missing stratum in the processing sequence generated by Step 7, each missing stratum is treated as followings: If Si is a missing stratum of Type 1, as shown in Fig. 5C, the elevation of each vertex in F(i+1)_Bottom (the bottom TIN of Si+1) is set equal to the elevation of the corresponding vertex in F’i_Top (the top TIN of Si), and the modified F(i+1)_Bottom is denoted as F’(i+1)_Bottom; if Si is a missing stratum of Type 2, as shown in Fig. 6C, the elevation of each vertex in F(i-1)_Top (the top TIN of Si-1) is set equal to the elevation of the corresponding vertex in F’i_Bottom (the bottom TIN of Si), and the modified F(i-1)_Top is denoted as F’(i-1)_Top.
After above process, the consistency of the surfaces for all adjoining strata is ensured. That is, no matter the complete stratum or the missing stratum, the top TIN of each stratum is identical with the bottom TIN of its overlying stratum, and the bottom TIN of each stratum is identical with the top TIN of its underlying stratum.
After previous steps, construction of the top and bottom surfaces for each stratum is accomplished, and a set of TINs composed of triangular patches are generated. In this step, we need to establish a solid model by the topological relationship between stratum and each surface. This is a relatively easy task. For each stratum, the top and the bottom TINs are extruded and a block is built (Lemon and Jones, 2003). All these blocks can be combined into one solid model. A voxel-based solid model that filled with triangular prism meshes can also be established by applying 3D discretization methods.
Finally, the modeling result is used for 3D visualization and spatial analysis. Several operations for 3D-interaction of the solid model, such as 3D observation, slice up, arbitrary incision, virtual drilling, virtual roaming, spotting and measurement of property value in any spatial position, excavation of foundation pit or tunnel, etc., can be performed freely since the solid model is very suitable for spatial analysis and spatial query.
5. Verification and validation
The Boreholes-Surfaces-Solids method has been programmed in Visual C++ and the OpenGL graphics library on a PC platform, and has been integrated into a 3D Geological Modeling and Visualization System (Zhu et al., 2006) which is based on MapGIS, one of the widely used GIS softwares in China. Two case studies with different geological settings are given below to illustrate the feasibility and practicability of the presented method.
The first study area, which is located in Shanghai Pudong New District, China, is a part of the site area for the World Expo 2010 Shanghai (Shanghai Geotechnical Investigations & Design Institute Ltd., 2008). As Fig.7A shows, the data set for 3D geological modeling consists of 7 shallow boreholes in the area of about 250000 m2 (500×500m2), and 6 stratigraphic units are detected. The strata are denoted as S1, S2, S3, S4, S5 and S6 from the bottom to the top. S1, S4, S5 and S6 are complete strata, while S2 and S3 are incomplete. The modeling results utilizing the Boreholes-Surfaces-Solids method are shown in Fig.7B to Fig.7F. Fig.7B displays the spatial distribution of the top surface for each stratum. Fig.7C shows the solids filled with triangular prism meshes. Fig.7D shows the solids represented as 3D blocks. Fig.7E shows a 3D cross-section which is created from the solids and passes through borehole B2, B5 and B4. Fig.7F shows a fence diagram cut from the solids.
(Fig.7. Example of missing stratum induced by erosion)
This model is a typical example of sedimentary system controlled by fluvial erosion and aggrading action. From the solid model, we can see that S2 and S3 are missing strata of Type 2 that induced by erosion, and the surface of each stratum is natural and geologically reasonable.
In order to quantitatively evaluate the accuracy of the solid model, a set of additional borehole data and excavation data obtained from the practical construction are compared with the solids. The comparison result shows that the error between the computer model and the measured data is within 5cm to the surface of complete strata, while 8cm for the surface of missing strata. Therefore, we can say that the solid model has a higher accuracy and can be used for practical projects directly.
In addition, the implementation algorithm of the Boreholes-Surfaces-Solids method was tested with the borehole data in this site using several different primary TINs. The testing was performed on a PC with Intel Core i7-740QM 1.73GHz CPU and 4G Memory. Based on 7 boreholes in the study area, we created six primary TINs as the testing dataset, and all TINs had the same outer boundary, but differed in the number of triangles. The computation times for the various TINs are shown in Table 1. The corresponding relationship between the computation time and the number of triangles is shown in Fig.8. From Table 1 and Fig.8, we can see that the computation time of the Boreholes-Surfaces-Solids method is linear with the number of triangles in the primary TIN, and it is in line with our expectations.
(Table 1. Testing dataset and computation time)
(Fig.8. Computation time linear with the number of triangles)
The second study area is located on the riverside of Suzhou River in Shanghai Putuo District, China, and covers approximately 800×600m2 (Shanghai Geotechnical Investigations & Design Institute Ltd., 2010). As Fig.9A shows, there are 6 shallow boreholes detecting 4 stratigraphic units. The strata are denoted as S1, S2, S3 and S4 from the bottom to the top. S1 and S4 are complete strata, while S2 and S3 are incomplete. The modeling results are shown in Fig.9B to Fig.9F. Fig.9B displays the spatial distribution of the top surface for each stratum. Fig.9C shows the solids filled with triangular prism meshes. Fig.9D shows the solids represented as 3D blocks. Fig.9E shows a 3D cross-section which is created from the solids and passes through borehole C2, C6, C5 and C4. Fig.9F shows a fence diagram cut from the solids.
(Fig.9. Example of missing stratum induced by combined action)
Although there are only four stratigraphic units in this site, the geological setting is more complex than the previous one as the boundaries of the missing strata are controlled by superposition and compound of non-deposition and erosion. From Fig.9E, we can see that in the location of borehole C2, S2 and S3 are missing induced by non-deposition, while in the location of borehole C4, S2 and S3 are missing induced by erosion. Thus, S2 and S3 are missing strata of Type 3. In addition, the practical project also proves that the solid model is reliable. This example demonstrates that the Boreholes-Surfaces-Solids method is an effective, reasonable method to handle complex sedimentary system, especially in the areas controlled by superposition and compound of non-deposition and erosion.
6. Conclusions
The main difficulty raised by 3D geological modeling of sedimentary strata system is determining the geological genesis and geometrical boundaries of missing strata. For the lack of the perfect mechanism to handle missing strata, it is difficult to construct spatial geometric shapes with a desired accuracy of complicated strata in 3D utilizing the existing modeling methods. In this paper, an adapted method to 3D solid modeling of sedimentary system from borehole data, called the Boreholes-Surfaces-Solids method, is described and applied to the construction projects in Shanghai, China. The Boreholes-Surfaces-Solids method overcomes limitations of the previous modeling methods, and the most significant feature is that a comprehensive treatment of all types of missing strata is considered, thus the boundaries of missing strata can be precisely controlled just using borehole data. Compared with the previous methods, the substantial advantages and benefits of the Boreholes-Surfaces-Solids method are obvious:
(1) The modeling process is automatic, simple and intuitive. In this method, both explicit and implicit geological information that detected by boreholes are used effectively. We can automatically estimate and deduce the geological genesis of the missing strata, intersect surfaces, adjust elevations and keep consistencies for the missing strata and their adjoining strata. Thus, the geometrical boundaries of missing strata are extracted automatically.
(2) The implementation algorithm is robust, time efficient and suitable in both simple and complex geological setting. The Boreholes-Surfaces-Solids method can not only effectively handle the missing strata induced by non-deposition or erosion respectively, but also automatically handle the missing strata induced by the superposition and compound of non-deposition and erosion. This method is flexible since it can deal with complicated sedimentary system composed of any number of stratigraphic units.
(3) The modeling results are natural, reliable and geologically reasonable. During the modeling process, the geological laws, such as the different geological genesis controlling the spatial geometry of missing strata, are converted skillfully into modeling rules that can be identified and programmed automatically by modelers. Thus, the geometrical boundaries of the missing strata can be precisely controlled with borehole data, and this method can automatically simulate all types of missing strata no matter the sample data are sufficient or not.
At present, societal expectation for sustainable development and continued environmental protection raise demands for more complex and quantitative assessments of subsurface conditions (Turner, 2006). To address this subsurface characterization need, remarkable advances have been made in 3D geoscience modeling technique over the past three decades. However, there are still three major technical challenges in subsurface geological modeling. These challenges are: How to detect geological laws hidden under large numbers of geological data? How to convert these ambiguous geological laws into rigorous modeling rules that can be identified and programmed automatically by computers? How to develop specialized tools for exploring and modeling complex geological systems? Maybe the research approach in this paper can be helpful for promoting the further research and development on these subjects.
Financial support for this work, provided by the National Natural Science Foundation of China (Grant No.40902093), the National Science and Technology Program of China (Grant No.SinoProbe-08-02), the Open Foundation of Shanghai Key Lab for Urban Ecology and Sustainability (Grant No.2011-18) and the Morning Light Plan of the Shanghai Educational Development Foundation (Grant No.2007CG34), is gratefully acknowledged.
Chang, Y., Park, H., 2004. Development of a web-based Geographic Information System for the management of borehole and geological data. Computers & Geosciences 30(8): 887-897.
De Rienzo, F., Oreste, P., Pelizza, S., 2008. Subsurface geological-geotechnical modelling to sustain underground civil planning. Engineering Geology 96(3-4): 187-204.
Hack, R., Orlic, B., Ozmutlu, S., Zhu, S., Rengers, N., 2006. Three and more dimensional modeling in geo-engineering. Bulletin of Engineering Geology and the Environment 65(2): 143-153.
He, H.J., Bai, S.W., Zhao, X.H., Chen, J., 2002. Discussion on strata partition in three dimension strata model. Rock and Soil Mechanics 23(5): 637-639 (in Chinese).
He, M.C., Li, X.Y., Liu, B., Xu, N.X., 2005. Study on processing method of drilling data for three dimensional modeling of engineering rock mass. Chinese Journal of Rock Mechanics and Engineering 24(11): 1821-1825 (in Chinese).
Jones, C.B., 1989. Data structures for three-dimensional spatial information systems in geology. International Journal of Geographical Information Systems 3(1): 15-31.
Jones, T.A., 1988. Modeling geology in three dimensions. Geobyte 3(1): 14-20.
Lees, J.M., 2000. Geotouch: software for three and four dimensional GIS in the earth sciences. Computers & Geosciences 26(7): 751-761.
Lemon, A.M., Jones, N.L., 2003. Building solid models from boreholes and user-defined cross-sections. Computers & Geosciences 29(5): 547-555.
McCarthy, J.D., Graniero, P.A., 2006. A GIS-based borehole data management and 3D visualization system. Computers & Geosciences 32(10): 1699-1708.
Nathanail, C.P., Rosenbaum, M.S., 1998. Spatial management of geotechnical data for site selection. Engineering Geology 50(3-4): 347-356.
Royse, K.R., Rutter, H.K., Entwisle, D.C., 2009. Property attribution of 3D geological models in the Thames Gateway, London: new ways of visualising geoscientific information. Bulletin of Engineering Geology and the Environment 68(1): 1-16.
Shanghai Geotechnical Investigations & Design Institute Ltd., 2008. Geologic Survey Report of the World Expo Park in Shanghai Pudong New District, China. Final Report, Shanghai Geotechnical Investigations & Design Institute Ltd., China, pp.75-90 (in Chinese).
Shanghai Geotechnical Investigations & Design Institute Ltd., 2010. Geologic Survey Report of Changfeng CBD in Shanghai Putuo New District, China. Final Report, Shanghai Geotechnical Investigations & Design Institute Ltd., China, pp.34-42 (in Chinese).
Smirnoff, A., Boisvert, E., Paradis, S.J., 2008. Support vector machine for 3D modelling from sparse geological information of various origins. Computers & Geosciences 34(2): 127-143.
Tremblay, T., Nastev, M., Lamothe, M., 2010. Grid-based hydrostratigraphic 3D modelling of the Quaternary sequence in the Chateauguay River Watershed, Quebec. Canadian Water Resources Journal 35(4): 377-398.
Turner, A.K., 2006. Challenges and trends for geological modelling and visualization. Bulletin of Engineering Geology and the Environment 65(2): 109-127.
Wellmann, J.F., Horowitz, F. G., Schill, E., Regenauer-Lieb, K., 2010. Towards incorporating uncertainty of structural data in 3D geological inversion. Tectonophysics 490(3-4): 141-151.
Wu, L.X., 2004. Topological relations embodied in a generalized tri-prism (GTP) model for a 3D geoscience modeling system. Computers & Geosciences 30(4): 405-418.
Wu, Q., Xu, H., 2004. On three-dimensional geological modeling and visualization. Science in China Series D: Earth Sciences 47(8): 739-748.
Xu, N.X., He, M.C., 2004. 3D modeling methods and spatial data model of layered rock-mass. Journal of China University of Mining & Technology 33(1): 103-108 (in Chinese).
Zhang, F., Zhu, H.H., Ning, M.X., 2006. Modeling method of 3D strata suitable for massive data. Chinese Journal of Rock Mechanics and Engineering 25(Supp.1):3305-3310 (in Chinese).
Zhong, D.H., Li, M.C., Song, L.G., Song, L.G., Wang, G., 2006. Enhanced NURBS modeling and visualization for large 3D geoengineering applications: An example from the Jinping first-level hydropower engineering project, China. Computers & Geosciences 32(9): 1270-1282.
Zhu, H.H., Wu, J.B., 2005. 2D and 2.5D Modeling of strata based on Delaunay triangulation. Chinese Journal of Rock Mechanics and Engineering 24(22): 4073-4079 (in Chinese).
Zhu, L.F., Pan, X., 2005. Reconstruction of 3D stratigraphic model for fluvial erosion and aggrading action. Rock and Soil Mechanics 26(Supp.1): 65-68 (in Chinese).
Zhu, L.F., Wu, X.C., Liu, X.G., Shang, J.G., 2004. Reconstruction of 3D strata model based on borehole data. Geography and Geo-Information Science 20(3): 26-30 (in Chinese).
Zhu, L.F., Wu, X.C., Pan, X., 2006. Mechanism and implementation of error correction for 3D strata model. Rock and Soil Mechanics 27(2): 268-271 (in Chinese).