Orginal Article

High Quality Geological Surface Reconstruction from Planar Contours

  • YANG Yang , 1 ,
  • PAN Mao , 1, * ,
  • WU Gengyu 1 ,
  • SUN Ying 2 ,
  • LI Kuixing 3
Expand
  • 1. Key Laboratory of Orogenic Belts and Crustal Evolution,Ministry of Education, Peking University, Beijing 100871, China
  • 2. Beijing Institute of Hydrogeology and Engineering Geology, Beijing 100195, China
  • 3. Creatar Co., Ltd., Beijing 100085, China
*Corresponding author: PAN Mao, E-mail:

Received date: 2014-05-09

  Request revised date: 2014-07-27

  Online published: 2015-03-10

Copyright

《地球信息科学学报》编辑部 所有

Abstract

In three-dimensional geological modeling, because the geological section and the medical section are essentially similar, the three-dimensional surface reconstruction methods based on contour lines that have been widely used in biomedical modeling are now introduced to geological modeling. But most algorithms from medical modeling that have been applied to geological modeling do not fit specifically for geological data. They only consider the rationality of reconstructed surface, however they concern little about the quality of surface geometry. In geological modeling, there are challenges including the complex and changeable shapes of geological objects, sparse geological section data and various origins of data sources; therefore, the methods that simply connect data could not meet the quality requirements of three-dimensional geological modeling. To overcome it, we consider solving the problem of how to reconstruct high quality triangular surfaces from sparse contour lines in this paper. According to the characteristics of geological data, a new algorithm that improves the quality of geological modeling by interpolating transitional sections is presented. The algorithm deals with geological sections that are stored as vectors and integrate them through a series of manipulations, including: matching features between contour lines that have the same geological property; building a mapping function between each matched pair; generating transitional contour lines; constructing surfaces from transitional contours and original contours. The main process can be summarized as follows: (1) matching geological feature; (2) building mapping function; (3) generating transition model; (4) constructing surface. A standard mathematical model which can resolve contour lines from parallel sections and a extended model that can resolve contour lines from cross sections are both defined in this paper. Some illustrative examples and analytical data are given at the end of this paper to prove that the surface reconstructed by this algorithm has much better quality than traditional ones.

Cite this article

YANG Yang , PAN Mao , WU Gengyu , SUN Ying , LI Kuixing . High Quality Geological Surface Reconstruction from Planar Contours[J]. Journal of Geo-information Science, 2015 , 17(3) : 253 -259 . DOI: 10.3724/SP.J.1047.2015.00253

1 引言

在三维地质信息系统中,要直观和准确地表达地下空间,就需要生成便于显示或分析的三维数据。目前,对于不规则对象的三维数据获取和生成的方法有许多,归纳起来主要有轮廓线的三维重建、图像数据的三维重建、点云的三维重建、专业领域的三维重建等。这些三维重建方法得到的主要是表面数据。在三维地质建模中,有很大一部分数据来源是地质剖面,针对地质剖面轮廓线进行地质表面的重构,一直是三维地质建模的研究热点。实际上,在计算机图形学研究中[1],轮廓线重建三维表面已经被普遍研究,相关方法广泛应用于医学、生物、生产制造等领域,是科学计算可视化的主要内容之一[2]。虽然国内外对于此类算法研究很多,但是大都围绕医学图像三维重建,如利用CT或MRI断层切片构造三维人体或物体表面[3]
在地质领域,由于地质剖面与医学CT有本质上的相似性,因此,轮廓线的三维重建在矿体或矿山三维模型的构建[4-5]、平行剖面建模[6]、交叉剖面建模[7-8]等技术中有广泛的应用,目前,较为成熟的算法主要有最大体积法[9]、最小表面积法[10]、边长最小法[11]、最短对角线法[12]、同步前进法[13]和切开缝合法[14]等,但是,这些成果并没有考虑到地质数据的特殊性,将普适的算法应用于地质建模中,只是解决了重建曲面的合理性,并没有关注曲面的几何质量[15]
在计算机图形学中对曲面的质量有明确的定义,其包括美化度(fairness)和光滑度(smoothness)2个部分,美化度关注曲率的分布,而光滑度关注曲率的变化[16]。当然,美化度与光滑度并不是绝对独立的,往往提升其中一者,另一者也会有所改善。轮廓线建模中使用的不规则三角网(TIN),其美化度有更直接的定义:三角形的形态[17]。对于TIN,美化度比光滑度更为重要,因为当遇到退化三角形(图1)时,最基本的几何计算(如法向、曲率)精度会大大降低甚至失效,消除这些三角形可能付出极大的代价[18]。因此,在生成曲面的时候,必须对曲面的质量有基本的保证。
Fig. 1 Classifications of skinny triangles

图1 3种典型的退化三角形[18]

三维地质建模中,地质概念的表达与几何模型有密不可分的关系,目前的轮廓线建模算法将关注点放在轮廓线连接的合理性上[9-14],解决了地质概念表达的正确性,但是忽略了几何模型的质量问题。与医学、生物等领域不同,地质数据存在3个特殊性[15]:(1)信息不完全;(2)数据获取难;(3)数据稀疏。这些数据特性导致重构的曲面在质量上存在缺陷,过于稀疏的数据会产生大量退化三角形(图2),这将降低建模后期算法的准确性,甚至限制了算法的应用。
Fig. 2 This tessellation of a fold from traditional contour model has a large number of skinny triangles

图2 传统轮廓线建模会产生大量退化三角形

本文针对地质数据的特点,提出了一种新的轮廓线方法。在不改变轮廓线连接合理性与正确性的情况下,通过生成过渡剖面轮廓线,提高了生成曲面的质量,在数据稀疏的情况下达到了良好的建模效果,并且其扩展的模型对人工辅助信息也有很好的支持。

2 高质量轮廓线建模

地质剖面之间往往具有非常大的间隔尺度,剖面之间直接进行三角形构网,生成的三维地质模型比较粗糙图3(c),所以,可考虑能否在两两剖面之间生成更多的过渡剖面,以提高生成模型的质量。
由于地质剖面数据往往是以矢量的形式存储的,因此,首先找到起始剖面和目标剖面中用于表达相同地质曲面的轮廓线,将其进行匹配,在每一对轮廓线之间构造映射函数,然后生成过渡轮廓线,最后连接轮廓线形成整个曲面。其步骤为:(1)特征匹配;(2)构造映射方法;(3)生成过渡模型。
设{S1,S2,S3SK}是一组有序的近似平行的剖面,其中,Sk={C1,k,C2,kCL,k},即剖面Sk由有限条轮廓线 Ci,k组成。现考虑{S1,S2,S3SK}中表示相同地质界面的轮廓线集合{Ci,1,Cj,2Cw,K}如图3(a)所示。若在两两轮廓线Ci,k-1Cj,k间构造一种变形映射函数f,对给定的轮廓线采用该映射函数进行相关的变换处理,以生成所需的中间过渡轮廓线(图3(b)),然后在过渡轮廓线参与的情况下,就可生成高质量的曲面(图3(d))。
Fig. 3 An example of geological surface reconstruction from planar contours

图3 轮廓线三维地质建模应用示例

如果将剖面中的其他弧段也做如此处理,就可生成过渡剖面,进而生成高质量的三维地质体。下面按照这种思路,首先将建立一种标准模型,并构建相应的算法,然后,在标准模型的基础上再建立扩展的标准模型及相应算法。

2.1 标准模型

可将图3(a)的模型分解为两两轮廓线间的构网处理问题,现分别定义为源地质界线A={a1,a2an}和目标地质界线B={b1,b2bm},那么,AB就组成了标准模型(图4),其具有以下特点:
(1)A={a1,a2an}共面于Pa,B={b1,b2bm}共面于Pb
(2)PaPb近似平行。
标准模型的求解可分为3步:
① 建立源轮廓线和目标轮廓线的对应关系,折线的对应关系其实就是点的对应关系,按照建立点对应关系的算法,可以得到AB之间的对应关系。
② 构造映射函数f,该映射函数实际上就是将源轮廓线上的每一个线段根据一定的数学规律映射到目标轮廓线上的线段或者点。
③ 插值确定坐标,根据空间地理信息生成中间的曲线。
2.1.1 建立点对应关系
建立点对应关系[19]的过程实际上是连接轮廓线的过程,因此,在建立点的对应关系算法中,借鉴轮廓线连接法的思想加以改造,构造了建立点对应关系的算法。本文同步前进法[13]进行改造,形成点对应关系建立的算法。
该方法的基本思想是在连接2条相邻轮廓线上的点时,连接操作在两条轮廓线上尽可能同步进行。对轮廓线上的每个线段赋以权值 φ,其含义是该线段长度除以该线段所在轮廓线的总长度所得的值。
根据标准模型(图4)定义,在图5中可以假设:ΦA表示源地质界线A中已经存在的权值之和,ΦB表示目标地质界线B中已经存在的权值之和,AB的连接分别进行到点ai和点bj,线段aiai+1bjbj+1的权值分别为 φ A φ B 。则此时下一步的对应关系有:ai+1bj,aibj+1。如果 Φ A + φ A - Φ B < Φ B + φ B - Φ A ,则选择ai+1bj,即沿A前进一步;反之,则选择aibj+1,沿B前进一步。这样最多可通过m+n步操作之后,则建立了的对应关系CB,定义如下:CB:A→B。
Fig. 4 Illustration of standard model

图4 标准模型示意图

Fig. 5 Illustration of the correspondences between points

图5 点对应关系示意图

2.1.2 构造映射函数
在建立源与目标的匹配C之后,接下来需要构造变形映射函数f,这是本算法中最重要的步骤,若将A到B的变形视为一种内插,那么f可以有任意形式的表达,可以是线性,也可以为非线性。由于线性插值在地质建模中最为普适,本文根据线性内插函数的表达形式,对源U及其匹配C(U)的映射函数 f 作如下定义:
f ( U , t ) = tU + ( 1 - t ) C ( U ) , t [ 0,1 ] (1)
用类似的函数结构,可以构造剖面轮廓线间的映射函数,其中,t表示其与源U的近似率。若如图6所示,假设轮廓线A为源,轮廓线B为目标,建立2条轮廓线的匹配关系CB后,根据式(1),AB间的过渡轮廓线P可以表示如下:
P = tA + ( 1 - t ) C B ( A ) (2)
Fig. 6 Illustration of interpolation of the transition line in standard model

图6 标准模型过渡轮廓线求取示意图

假设需要构造n条过渡轮廓线,可认为n条过渡轮廓线与源的近似率逐条递减,即近似率tx的值与过渡轮廓线的序数x有关,由此可以得出第x条过渡轮廓线对于源的近似率tx的计算公式如下:
t x = ( n - x ) n (3)
根据式(2)则可求解 t x 时的过渡轮廓线Px={p1,p2,⋯}(图6),即对于任意ai(xi,yi,zi)及其匹配点CB(ai)=bj(xj,yj,zj),利用式(2)可知pk(xk,yk,zk)的计算公式如下:
p k = t x a i + ( 1 - t x ) b j (4)
图6所示,对匹配CB中所有匹配点求解,可得到过渡轮廓线Px,利用相同的方法,最终可构造出所需要的n条过渡轮廓线。对包括过渡轮廓线在内的所有轮廓线间应用点对应关系算法,则可产生一个质量好的地质曲面。

2.2 扩展模型

三维地质建模时,在2条平行剖面生成过渡剖面时,往往需要考虑其他的约束信息的影响,如地质专家的交互解译信息[4]或者多组交叉剖面的影响[20-21],只有将这些信息作为控制条件带入构建过程中,才能生成足够精确的地质界面。
若将加入的约束地质界线E,F分别视为源与目标,那么,可产生2个对应关系:CB:AB;CF:EF
图7所示,对于约束线E,F假设E={a1,e2b1}, F = { a n , f 2 b m } ,对于CF,exCF(ex)=fy是一对匹配点,那么,所求的过渡轮廓线Q,就是以ex(xe,ye,ze)为起点, f y ( x f , y f , z f ) 为末点产生既满足exfy趋势,又满足AB趋势的曲线,由于要满足2个趋势,因此,需要2次处理过程来构造所求轮廓线Q
Fig. 7 Illustration of interpolation of the transition line in extended model

图7 扩展模型过渡轮廓线求取示意图

首先,需要先构造对于CB趋势线P,与标准模型类似,但是,由于存在exfy的约束,P上各点对于A的近似率t会随着与ex,fy的远近而变化,因此,需重新定义近似率公式。
假设A长度为DA,在 a i 处的累积长度为dA,那么可以定义tA=dA/DA,类似可以定义tB,tE,tF。那么,式(2)中的近似率t可以定义如下:
t = t A t F + ( 1 - t A ) t E (5)
由式(2)、(5)可求得P,P满足了CB的趋势。下面需要将P变换到Q以满足exfy的趋势。
CP:PQP,Q间的一一对应,那么,对于任意点 p k ( x p , y p , z p ) 及其点对 C F ( p k ) = q k ( x q , y q , z q ) ,即定义:
Δ p k = q k - p k (6)
由于 q 1 = e x , q n = f y ,可得:
Δ p 1 = e x - p 1 (7)
Δ p n = f y - p n (8)
Δ p k 满足:
Δ p k = t A Δ p n + ( 1 - t A ) Δ p 1 (9)
将式(6)、(9)结合,可得:
q k = p k + t A Δ p n + ( 1 - t A ) Δ p 1 (10)
由于 p k 已由式(4)求得,根据式(10),最终可在满足 C F 的趋势下将P变换到Q
C F 中的每一对匹配都进行以上运算,可构造出所求的一系列过渡轮廓线。最后,利用轮廓线连接点对应关系算法,最终形成满足约束的地质界面。

3 轮廓线三维地质表面重建实验及分析

轮廓线三维地质表面重建,系利用Creatar XModeling三维地质建模平台SDK在Creatar XModeling三维地质建模基础平台上实现的。软件开发环境使用Visual C++ 2008和OpenSceneGraph 3.2.1进行搭建,硬件配置为Intel Core2 Duo 2.26 GHz CPU,2.0 GB RAM。为验证算法适用性,本文挑选了5组测试模型(表1)进行实验:对简单矿体、复杂矿体和水平地层数据使用标准模型进行构建,并与传统轮廓线连接算法进行了比较(图3、8、9);对于一些复杂情况如断层剖面线结合地表走势线、带趋势的复杂褶皱使用扩展模型进行构建实验(图10、11)。
Tab.1 Description of the testing model

表1 测试模型说明表

测试模型 模型说明
简单矿体(图3 共6条剖面轮廓线数据,且均为“1对1”的情况,无尖灭情况
水平地层(图8 共3条剖面轮廓线数据,趋势较为平滑不存在复杂构造
复杂矿体(图9 共13条剖面轮廓线数据,存在“n对n”,以及“1对n”的情况,并且沿矿体展布方向存在尖灭(“1对0”)的情况
断层面(图10 由2条剖面轮廓线(图10中实线),以及地表地质图上该断层的走势线(图10中虚线)组成
复杂褶皱(图11 由4条交叉剖面上的轮廓线组成
Fig. 8 Illustration of a stratigraphic model

图8 水平地层模型

Fig. 9 Reconstruction of a complex ore body model

图9 复杂矿体模型

Fig. 10 Reconstruction of a fault from fault trace and fault sticks

图10 剖面断层轮廓线结合地表地质图上断层走势线进行断层面重构

Fig. 11 Reconstruction of a fold contours from the traces of four cross sections

图11 利用交叉剖面上的轮廓线进行褶皱面重建

实验中,曲面质量由美化度与光滑度来评价,美化度的计算参考文献[17],光滑度可以用曲面平均离散曲率[22]来评价。其中,美化度 B [ 0,0.5 ] 越大说明曲面三角形形态越好,平均离散曲率H越小说明曲面越光滑。对于实验的模型,进行曲面质量评估(表2)。从表2可明显看到,对于地质剖面间距过大,数据稀疏的情况,传统算法无法保证模型质量,而使用本算法的模型质量得到极大改善。
Tab. 2 Quality evaluation diagram of the testing models

表2 模型质量评估表

测试模型 传统算法美化度B 传统算法光滑度H 本算法美化度B 本算法光滑度H
图3 0.161608 5.847635 0.299777 1.258108
图8 0.075692 27.52082 0.422103 5.292667
图9 0.144637 20.63944 0.369782 3.956764
图10 0.110934 2.144750 0.408326 1.545399
图11 0.097412 6.380472 0.399181 0.075950
在实际建模过程中,对于矿体的分支情况,处理的方法在文献[4]、[5]、[14]、[17]中都有涉及,其中心思想是通过增加辅助信息对轮廓线进行处理(图9中左侧分支为利用文献[17]方法进行了处理),最终将分支转化为多个“1对1”情况进行处理,因此,无论用何种方法处理分支情况,都可采用本文提出的算法生成高质量表面模型。对于地质中的尖灭情况,实际上是一种特殊的“1对1”情况,本文算法依然适用。以上情况在图9算例中得到很好的验证。

4 结论

本文引入地质约束的概念,针对地质数据的复杂特点提出扩展模型。地质约束的来源没有固定的形式,可以是同源数据(图11的例子就是4条交叉剖面组成的扩展模型),也可以是不同源数据(图10的例子就是剖面断层线与其地表地质图上走势线结合而成的扩展模型),甚至可以是专家解译的人工干预数据。
综上所述,轮廓线三维地质表面重建算法,针对地质数据的稀疏性、复杂性、多源性等特点,极大提高了生成模型的质量;扩展模型的提出,解决了地质复杂数据的解译生成。算法复杂度低,易于编程,适用于有良好趋势的地质单体表面的构建,为整个地质建模过程提供了良好的质量基础。
本文算法不足之处在于:对剖面数据的假设与传统轮廓线建模算法相同,需要相邻建模剖面近似平行,这在一定程度上限制了其适用性;不能够提高传统轮廓线表面重建算法的准确性与合理性;在推导变形映射函数时,使用的是最简单的线性内插,导致产生的过渡剖面有可能并不符合地质认识。若能在今后的工作中针对以上不足,寻找符合地质现象的变形映射函数,使得构建的地质表面更加符合地质特征,提高地质曲面重建的合理性与正确性,同时考虑以折剖面或任意剖面数据作为输入,那将会在地质建模中发挥更大的作用。

The authors have declared that no competing interests exist.

[1]
孙正兴. 计算机图形学教程[M].北京:机械工业出版社,2006.

[2]
石教英,蔡文立.科学计算可视化算法与系统[M].北京:科学出版社,1996.

[3]
赖大进,余轮.在Windows环境中实现医学图像三维重建[J].福州大学学报:自然科学版,2002,30(2):201-205.

[4]
李梅,毛善君,马蔼乃.平行轮廓线三维矿体重建算法[J].计算机辅助设计与图形学学报,2006,18(7):1017-1021.

[5]
南格利. 矿体线框模型及其建立方法[J].有色矿山,2001,30(5):1-4.

[6]
屈红刚,潘懋,王勇,等.基于含拓扑剖面的三维地质建模[J].北京大学学报(自然科学版),2006,42(6):717-723.

[7]
陈国良,刘修国,盛谦,等.一种基于交叉剖面的地质模型构建方法[J].岩土力学,2011,32(8):2409-2415.

[8]
屈红刚,潘懋,明镜,等.基于交叉折剖面的高精度三维地质模型快速构建方法研究[J].北京大学学报:自然科学版,2009,44(6):915-920.

[9]
Keppel E.Approximating complex surfaces by triangulation of contour lines[J]. IBM Journal of Research and Development, 1975,19(1):2-11.

[10]
Fuchs H, Kedem Z M, Uselton S P.Optimal surface reconstruction from planar contours[J]. Communications of the ACM, 1977,20(10):693-702.

[11]
Christiansen H N, Sederberg T W.Conversion of complex contour line definitions into polygonal element mosaics[J]. ACM Siggraph Computer Graphics, 1978,12(3):187-192.

[12]
Ekoule A, Peyrin F, Odet C.A triangulation algorithm from arbitrary shaped multiple planar contours[J]. ACM Transactions on Graphics (TOG), 1991,10(2): 82-99.

[13]
Ganapathy S, Dennehy T.A new general triangulation method for planar contours[J]. ACM Siggraph Computer Graphics, 1982,16(3):69-75.

[14]
马洪滨, 郭甲腾. 一种新的多轮廓线重构三维形体算法: 切开-缝合法[J].东北大学学报:自然科学版,2007,28(1):111-114.

[15]
Pan M, Li Z L, Gao Z B, et al.3-D Geological Modeling - Concept, Methods and Key Techniques[J]. Acta Geologica Sinica-English Edition, 2012,86(4):1031-1036.

[16]
Botsch M, Pauly M, Rossl C, et al.Geometric modeling based on triangle meshes[C]. Proceedings of the ACM SIGGRAPH 2006 Courses, Boston, Massachusetts, 2006.

[17]
Mallet J L.Geomodeling[M]. New York: Oxford University Press, 2002.

[18]
Botsch M, Kobbelt L.A robust procedure to eliminate degenerate faces from triangle meshes[C]. Proceedings of the VMV, 2001:283-290.

[19]
邓敏,彭东亮,徐震,等.一种基于弯曲结构的线状要素Morphing方法[J].中南大学学报(自然科学版), 2012,43(7):2674-2682.

[20]
明镜. 三维地质结构建模系统及关键建模方法研究[D].北京:北京大学,2010.

[21]
明镜,颜玫.基于Morphing的三维地质界面生成[J].地理与地理信息科学,2014,30(1):37-40.

[22]
Taubin G.Geometric signal processing on polygonal meshes[C]. Eurographics State of the Art Reports, 2000,4(3):81-96.

Outlines

/