CN111968231B - 一种基于地质图切剖面的三维地层建模方法 - Google Patents
一种基于地质图切剖面的三维地层建模方法 Download PDFInfo
- Publication number
- CN111968231B CN111968231B CN202010819474.4A CN202010819474A CN111968231B CN 111968231 B CN111968231 B CN 111968231B CN 202010819474 A CN202010819474 A CN 202010819474A CN 111968231 B CN111968231 B CN 111968231B
- Authority
- CN
- China
- Prior art keywords
- stratum
- dimensional
- point
- contour
- boundary
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 80
- 238000005520 cutting process Methods 0.000 title claims abstract description 7
- 238000005070 sampling Methods 0.000 claims abstract description 35
- 238000004364 calculation method Methods 0.000 claims abstract description 17
- 238000000605 extraction Methods 0.000 claims abstract description 6
- 230000015572 biosynthetic process Effects 0.000 claims description 27
- 230000008569 process Effects 0.000 claims description 16
- 238000005516 engineering process Methods 0.000 claims description 11
- 239000011159 matrix material Substances 0.000 claims description 10
- 238000012545 processing Methods 0.000 claims description 7
- 230000008030 elimination Effects 0.000 claims description 5
- 238000003379 elimination reaction Methods 0.000 claims description 5
- 238000006243 chemical reaction Methods 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 238000012887 quadratic function Methods 0.000 claims description 2
- 229910002056 binary alloy Inorganic materials 0.000 claims 1
- 238000009826 distribution Methods 0.000 abstract description 5
- 230000006870 function Effects 0.000 description 15
- 238000010586 diagram Methods 0.000 description 7
- 230000000694 effects Effects 0.000 description 6
- 229910052500 inorganic mineral Inorganic materials 0.000 description 4
- 239000011707 mineral Substances 0.000 description 4
- 238000012800 visualization Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 239000000284 extract Substances 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 206010063659 Aversion Diseases 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 238000013499 data model Methods 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000003708 edge detection Methods 0.000 description 1
- 238000005429 filling process Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000012856 packing Methods 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 230000008439 repair process Effects 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
- G06T17/205—Re-meshing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4023—Scaling of whole images or parts thereof, e.g. expanding or contracting based on decimating pixels or lines of pixels; based on inserting pixels or lines of pixels
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Graphics (AREA)
- Remote Sensing (AREA)
- Geophysics And Detection Of Objects (AREA)
- Investigation Of Foundation Soil And Reinforcement Of Foundation Soil By Compacting Or Drainage (AREA)
- Processing Or Creating Images (AREA)
Abstract
本发明公开了一种基于地质图切剖面的三维地层建模方法,包括:S1、对图切剖面进行编辑,完成建模所需要的基本参数的设定;S2、对地层边界及内外区域进行离散采样,获得地层属性在空间的离散分布;S3、对地层三维空间离散点进行插值计算,生成三维空间地层属性数据场;S4、从三维空间地层属性数据场中进行等值面提取,获得地层的三维轮廓表面;S5、重复步骤S2~S4,获得所有地层的三维轮廓表面;S6、将所有地层三维轮廓表面进行空间拼接,完成三维地层建模。本发明既避免了二维轮廓重建方法中存在的复杂轮廓匹配、分叉及对应三角网连接等问题,又可以得到较高精度的三维矢量地质模型;本发明适用于任意复杂情况下的地层建模,建模效率高、适用范围广。
Description
技术领域
本发明涉及地质、矿产及工程领域,具体涉及的是一种基于地质图切剖面的三维地层建模方法。
背景技术
三维地质建模一直是工程、地质、矿产行业研究的技术热点和难点问题。通过有限的地质数据还原地下的地质构造形态和岩石物性参数,对于工程地质分析、油气及矿产勘探与开发具有非常重要的作用。由于地质数据的多元性、多解性及数据获取手段的有限性,使得三维地质建模难度大、实现过程复杂。其中,采用地质图切剖面进行三维地质建模是一种常见的建模方式。
图切剖面建模即是在地质图上绘制一系列的平行(或近似平行)剖线,采用地质分析的方法,将三维地质目标分割成一系列的二维地质剖面,然后采用三维重建手段构建三维目标地质体模型。利用图切剖面建模的技术主要有两种方式,一种是基于二维轮廓的三维重建技术;一种是对二维剖面进行栅格化的三维网格建模技术。
基于二维轮廓的三维表面重建技术基本实现思路是从三维物体中截取一系列的二维平行截面,然后针对特定的目标提取出每个截面的边界轮廓线,在三维空间中对相邻的轮廓线进行匹配和对应的三角连接,从而完成三维物体表面的建模。基于二维剖面的网格化建模是根据地质体的属性分类对二维剖面进行栅格化,然后再映射到三维空间中,生成三维网格化地质模型。相比较网格化的三维地质模型,基于二维轮廓的三维表面重建方法生成是一种三维矢量数据模型,模型精度高,可视化效果好,是最常用的三维建模技术。
在基于平行轮廓线重建三维表面的研究中,一般只考虑两个相邻的平行面上的轮廓线重建问题。因为多层轮廓线的重建问题可以简化为多个相邻平行轮廓线重构表面的叠加,所以相邻截面的轮廓线的连接是完成建模的关键。轮廓线的连接主要包含三个问题:轮廓对应、轮廓拼接和分叉问题,其中轮廓对应和分叉问题是难点问题。现有的这类算法对相邻的轮廓线有诸多限制,如要求相邻切平面之间的间距小、对应轮廓线之间的覆盖程度高、形状相似等,而有关分支问题的算法则侧重于解决一条轮廓线与多条轮廓线间正确连接的问题。
如果同一层上的轮廓线数超过一个,则必须决定轮廓线的对应问题(不同层上的哪些轮廓线是属于同一表面的)。如果层与层之间的距离特别小,对应问题可以通过轮廓线的重合程度来解决[1];然而,当距离较大而重合程度不能得到正确结果时,对应问题就成了一个很难解决的问题。即使解决了对应问题,剖分时还需要解决分支问题(通过轮廓线对应问题的检查和分析,来判定一层上的一条轮廓线和另外一层上的多条轮廓线是否属于同一表面,要通过分支结构将一条轮廓线与多条轮廓线通过三角面片连接成物体表面)。许多文献都对分支问题和对应问题进行了探讨,但是至今还没有十分令人满意的解决方法[2][3][4]。Barequet提出一个基于动态规划算法的最优化三角剖分策略,即BPLI(Barequetps Piecewise Linear Interpolation)方法[5][6]。对于输入的两层轮廓线,该方法不直接致力于解决分支问题和对应问题,而是设法将它们连接成一个不出现自相交的三维表面。这一方法对于复杂的表面能够得到拓扑结构正确的剖分结果,从而可以巧妙地绕开对应问题和分支问题。
如何从众多的可接受表面组合中,确定一种需要的组合,目前采用的方法有:重建后的三维形体体积最大法[7]、表面积最小法[8]、边长短法[9]、最短对角线法[10],相邻轮廓线同步前进法[11],基于领域知识的构建算法[12]、切开-缝合法[13]等。对于解决形状和顶点数目差异较大的相邻轮廓线重构问题,这些方法在适应性、连接效果和效率方面方法均有不同程度的局限性。
何金国[14](2003)等对BPLI做出了改进尝试研究,提出了轮廓线分段匹配算法及求解空间多边形三角剖分的新算法,以消除退化区域,提高匹配效果;祁伟丽[15](2007)等提出了在中间层插入一个新的轮廓线的一种方法,可解决简单情况下分支对于问题;何畏[16](2010)等在基于最小表面积和最大体积法的基础上,提出了最优路径(OP)算法,对轮廓连接进行了优化尝试;陈晓青[17](2011)等提出了多边形曲线演化的方法,通过简化多边形并提取其特征点来实现相邻轮廓线的顶点匹配及连接问题。不过这种方法也是针对单轮廓对应的连接问题,对于较复杂的轮廓可能会出现狭长三角形;高士娟[18](2015)等介绍了借助GOCAD软件,以离散光滑插值(DSI)算法实现了单轮廓的地质体三维形态重构;杨洋[19](2015)等提出了轮廓线间插入过渡轮廓线的方式,以提高建模质量。荆永滨[20](2016)等,提出了一种轮廓线贴面算法,利用小波滤波器组将轮廓线简化为与原始轮廓线现状相似的低分辨率轮廓线,然后对低分辨率轮廓线进行最优贴面,然后构建轮廓之间的连接三角形。这种算法对于一定程度上降低了有向环路图算法的复杂程度,但是对于差异较大的相邻轮廓线则很难通过滤波得到相似的贴面。
由于BPLI方法建模拓扑结构复杂,对于复杂的轮廓构造,特别是对于某些切片距离比较大的复杂轮廓数据很可能会出现相交的三维表面模型,得不到正确的三维模型结果。采用栅格化的建模方式可以避免这些问题。
栅格化建模实现技术比较简单,只需要对二维剖面进行栅格化,然后通过三维插值算法生成三维空间的网格模型,最后从网格模型中提取出需要的地层网格。网格模型中三维地质体采用块体模型来进行表达,三维地质体被离散成空间上的规则几何形体(如六面体模型、四面体模型、棱柱模型等)。该建模方法不用考虑空间上的复杂拓扑结构和三角连接,实现方法简单,常用于矿产建模和储量估算等。但是该方法需要大量的二维轮廓截面,图切剖面工作量大,其次该方法生成的三维模型可视化效果较差。
王丹[21](2012)探讨了利用二维剖面,栅格化,映射到三维空间中,并采用八叉树生成三维栅格模型;周良辰[22](2013)等提出了采用加密了研究区域的产状数据,用三次样条拟合与人工修编结合的方法,对褶皱等复杂构造作处理,以提高栅格模型的质量;张宝一[23](2017)等介绍了利用大量的平行图切剖面,通过对图切剖面栅格化构建三维网格模型的过程;李陈[24](2108)提出了采用距离函数生成二维剖面的距离场,然后采用B样条插值方法生成三维属性场,在此基础上采用Marching Cubes算法提取物性等值面的方法生成地质体的轮廓表面。该方法避免了对复杂轮廓线的对应连接问题,可适应各种复杂情况下的地质建模。但是该方法通过距离场生成的三维物性等值面,很难找到一个准确的边界值,造成构建的三维模型边界误差较大,结果不准确。其次,由于剖面距离场是均匀剖分的网格,剖面的网格数据量会很大。如采用全局插值计算三维属性场,则需要耗费大量的计算机内存。按200×200剖分网格、20个剖面,反距离加权插值计算,则计算距离矩阵耗费计算机内存可高达2384G而采用局部插值算法虽然可以避免这种情况,但是局部插值方法(或线性插值方法)具有很明显的局限性。局部插值方法不会考虑数据分布的空间的整体关联性。在二维轮廓线稀疏的情况下,建模结果在轮廓之间会出现断裂,空间分布不连续的情况。
综上所述,现有的基于二维轮廓进行三维建模的方法中,基于二维轮廓的三维表面重建方法可以获得较高质量的三维地质模型。但是实现技术复杂,对于二维轮廓线要求较高。对应轮廓线要求简单凸多边形,非简单多边形需要进行多边形简化处理,增加了建模工作量。而对于复杂多轮廓的对应问题、差异较大、相隔较远的轮廓线会出现三角形相交、狭长三边形等连接问题,现有技术仍然没有一种很好的解决方案。而基于栅格化的剖面网格建模方法虽然可以避免以上问题,实现技术简单。但是由于栅格化网格模型受剖分影响,模型精度较低,可视化效果较差,在应用中具有一定局限。而基于三维属性场中提取等值面建模的方法,研究较少。现有的方法存在剖面采样精度较低、采样数据量大、三维属性场插值困难、等值面提取不准确等问题。
参考文献:
[1]Meyers D,Skinner S.Surfaces from Contours[J].ACM Transactions onGraphics.1992,11(3):228~258.
[2]Klein R,Schilling A,StraBer W.Reconstruction and Simplification ofSurfaces from Contours[J].Graphical Models.2000,62:429~443.
[3]Barequet G,Shapiro D,Tal A.Multilevel Sensitive Reconstruction ofPolyhedral Surfaces from Parallel Slices[J].The Visual Computer.2000,16:116~133.
[4]Levin D.Multidimensional Reconstruction by Set 2valuedApproximation[J].IMA J Numer Anal.1986,6:173~184.
[5]Barequet G,SharirM.Piecewise-linear Interpolation betweenPolygonal Slices[C].Proc 10th Ann ACM Symp onComputational Geometry(SoCG).1994.
[6]Barequet G,SharirM.Piecewise-linear Interpolation betweenPolygonal Slices[C].ComputerVision ImageUnderstanding.1996.63:251~272.
[7]Keppel E.Approximating complex surfaces by triangulationof contourlines[J].IBM Journal of Research and Development.1975,19(1):2-11.
[8]Fuchs H,Kedem Z M,Uselton S P.Optimal surfaces reconstruction fromplanar contours[J].Communication of the A CM.1977,20(10):693-702.
[9]Christiansen H N,Sederberg T W.Conversion of complex contour linedefinitions into polygonal element mosaics[J].Computer Graphics.1978,12(3):187-192.
[10]刘刚,胡远来,邓林.层迭三维体表面重构的算法探讨[J].成都理工大学学报:自然科学版,2003,30(5):537–540.(Liu Gang,Hu Yuan2lai,Deng Lin.Discussion onthe algorithm of reconstructing the superposition 32D surface[J].Journal ofChengdu University of Technology:Science&Technology Edition,2003,30(5):537-540.).
[11]唐泽圣.三维数据场可视化[M].北京:清华大学出版社.2000.
[12]王恩德,李艳,鲍玉斌.矿体三维可视化建模技术[J].东北大学学报:自然科学版,2005,26(9):890–892.
[13]马洪滨,郭甲腾.一种新的多轮廓线重构三维形体算法:切开-缝合法[J].东北大学学报.2007,1.Vol 28).No.1.
[14]何金国,查红彬.基于BPLI从二维平行轮廓线重建三维表面的新算法[J].北京大学学报(自然科学版).2003(03)
[15]祁伟丽,秦新强,王溪,宋丽平.基于二维平行轮廓线重建三维表面的算法研究[C].全国第18届计算机技术与应用(CACIS)学术会议,2007.8,p891-895.
[16]何畏,吴文鹂.基于截面轮廓线人机交互三维地质体建模[J],物探化探计算技术,2010,32(4):P433-436.
[17]陈晓青,任凤玉,张国建,丁航行.一种复杂矿体相邻断面的匹配算法[J].东北大学学报(自然科学版),2011,32(4):579-582.
[18]高士娟,毛先成,张宝一等.基于平面地质图的地质体三维建模[J].地质找矿论丛,2015,30(4):594-601.
[19]杨洋,潘懋,吴耕宇等.一种新的轮廓线三维地质表面重建方法[J].地球信息科学学报,2015,17(3):253-259.
[20]荆永滨王公忠孙光中.复杂矿体三维模型二维轮廓线重建方法[J].金属矿山,2016,285(11):124-127.
[21]王丹.基于平面地质图的三维地质建模方法研究[D].2012,南京师范大学硕士论文.
[22]周良辰,林冰仙,王丹等.平面地质图的三维地质体建模方法研究[J].地球信息科学学报,2013,15(1):46-54.
[23]张宝一,杨莉,陈笑扬等.基于图切地质剖面的区域成矿地质体三维建模与资源评价——以桂西南地区锰矿为例[J].吉林大学学报(地球科学版),2017,47(3):933-948.
[24]李陈.李陈.基于剖面的三维复杂地质体建模技术研究[D].成都理工大学硕士论文,2018.
[25]夏艳华,白世伟.基于水平集的复杂三维地层模型建模及在地下工程中的应用研究[J].岩土力学,2012,33(5):1445-1450.
[26]李程.基于MT算法的地质剖面绘制[J].数字通信,2014,41(4):35-41.
[27]帅仁俊,陈书晶.一种改进的MC三维重建算法[J].中国数字医学,2016,11(3):83-86.
[28]王铮,李瑞明.移动立方体算法面二义性问题研究[J].软件工程,2017,20(9):6-9.
[29]王铮.剖面数据场中的三维表面建模[J].软件工程,2018,21(3):10-15.
发明内容
本发明针对现有建模方法存在的问题,提供了一种基于地质图切剖面的三维地层建模方法,其将栅格化建模和地质模型表面重建结合在一起,既避免了复杂的轮廓对应、三角网连接等问题,又可以得到较高精度的三维矢量地质模型;并且还解决了建模过程中采样精度和数据量、全局插值方法三维数据场的计算以及等值面提取的精度等问题。
为实现上述目的,本发明采用的技术方案如下:
一种基于地质图切剖面的三维地层建模方法,包括以下步骤:
S1、对图切剖面进行编辑,完成建模所需要的基本参数的设定;对图切剖面进行编辑的流程包括:(1)通过标定空间定位点将剖面平面坐标转换为地质空间三维坐标;(2)通过地层属性填充获得地层参数信息;(3)通过地层边界追踪获得地层的轮廓边界;
S2、对地层边界及内外区域进行离散采样,并分别设定地层内部、地层边界以及地层外各自的采样点属性值,然后将所有图切剖面采样点转换成空间三维坐标点,组合成空间离散点集合;
S3、对三维空间离散点进行插值计算,生成三维空间地层属性数据场;
S4、从三维地层属性数据场中进行等值面提取,获得地层的三维轮廓表面;
S5、重复步骤S2~S4,获得所有地层的三维轮廓表面,然后执行步骤S6;
S6、将所有地层三维轮廓表面进行空间拼接,完成三维地层建模。
具体地,所述步骤S1中,采用基于线性队列的快速泛洪填充算法对地层进行填充,并记录其边界点,得到地层多边形,获得地层参数信息。
具体地,所述步骤S1中,对于追踪得到地层边界点,采用了投影处理,以消除由于栅格图像显示分辨率影响带来的边界点追踪误差,确保追踪得到的地层轮廓边界与原始地层轮廓边界的一致性。
进一步地,所述的投影处理包括以下步骤:
(1)将地层轮廓线进行简化,删除轮廓直线上的冗余点,然后将每两个点组成轮廓线段{Lk}(k=1,2,..N),轮廓线段数为N;
(2)选定一个追踪得到的地层边界点Pi(i=1,2,...,M),计算该点到所有轮廓线段的距离集合{Dik}(k=1,2,..N);
(3)从距离集合{Dik}中找出最小距离的线段Lj,将边界点Pi投影到线段Lj上,得到Pij,将Pi替换为Pij作为新的地层边界点;
(4)重复步骤(2)、(3),将所有地层边界点投影到地层轮廓线上;
(5)对投影后的地层边界轮廓线进行简化,删除直线上的冗余点,完成地层边界追踪。
具体地,所述步骤S2包括以下步骤:
S2-1、设定一个采样间隔step,并沿地层轮廓线边界进行边界点采样,得到地层边界点集合B1,并对集合中所有点属性值设为1;
S2-2、设定一个指定大小的剖分网格M×N,将图切剖面进行二维网格剖分,并记录每个网格中心点坐标,得到剖点集合B2;
S2-3、对集合B2中的点进行地层内外判别,若点在地层多边形内,则该点的属性值设置为1,否则置为0;
S2-4、对集合B1、B2中的所有点进行空间坐标转换,将平面坐标点转换成地质空间坐标点;
S2-5、将转换后的空间坐标点及其属性值放入散乱点数据集合{pi}(i=1,2,...,N),,得到地层采样后的离散空间点数据,完成采样。
作为优选,所述步骤S3中,以逆多元二次径向基函数作为插值函数,采用径向基函数插值算法,进行三维离散点空间插值,生成三维空间地层属性数据场;径向基函数形式如下:
其中Z0是待估测的插值点值,wi是相应点的权值,||pi-p0||是pi与p0之间的距离;
插值计算过程如下:
(2)采用高斯消元法结合OpenCL多线程加速技术,计算逆距离矩阵D-1;
作为优选,所述步骤S4中,采用MarchingCubes算法从地层三维属性数据场中提取等值面作为地层的三维表面轮廓,其过程如下:
(1)从三维网格数据中选取相邻8个点组成六面体;
(2)给定等值面值V0,对六面体中的每个节点状态值Vi进行判别,若Vi>V0,该节点状态设为1;否则设为0;i=1,2,...,8;
(3)根据节点状态序列,按照顺序组合成二进制字节,得到六面体单元构型表索引值Ti(0≤Ti≤255);
(4)根据构型索引值Ti,从构型表中检索出对应的棱边节点状态组合,计算棱边节点状态为0和1的等值面在棱边的交点坐标P(x,y,z),计算公式为:
p(x,y,z)=p1(x,y,z)+(v-v1)/(v2-v1)*(p2(x,y,z)-p1(x,y,z))
式中,p(x,y,z)是待求交点坐标,v是待计算等值面值;p1,p2是棱边两节点坐标,v1,v2是两节点属性值;
(5)根据构型表索引值Ti,从构型表中检索出对应的三角形顶点及连接顺序,组合成封闭的三角形面;
(6)遍历所有网格单元,重复步骤(1)~(5),得到目标地层表面轮廓的所有三角面。
具体地,所述步骤S6包括以下步骤:
(1)将所有顶点进行组合,并进行重复点剔除,得到新的顶点集合;
(2)根据新的顶点集合,修改所有模型的三角面片索引,完成三维地层建模。
与现有技术相比,本发明具有以下技术效果:
1、本发明能够适用于任意复杂情况的地层或地质体建模,相对于现有的基于二维轮廓线的建模方法,适应范围更广,原因在于:
(1)现有的基于二维轮廓线的三维建模方法,大都采取对空间相邻轮廓表面进行连接的建模方法,对图切剖面要求平行或近似平行,对于图切剖面的制作要求更高。而本发明对图切剖面的制作没有特别要求,可以采用任意方向、任意角度的图切剖面进行建模,甚至可以采取非直线(折线)的图切剖面。
(2)由于采用了全局插值算法对地层离散采样点进行了插值计算,因而本发明对二维图切剖面的间距(或密度)没有要求。剖面只要能控制待建模型的三维地质体的空间形态即可,可以是稀疏轮廓线。并且本发明以逆多元二次径向基函数作为插值函数,采用径向基函数插值算法进行三维离散点空间插值,可以获得更好的模型光滑性和空间连续性。
2、本发明建模效率更高,主要表现在:
(1)可以采用任意方向的图切剖面进行建模,大大减少了图切剖面的工作量,且更容易控制模型的精度。
(2)几乎不需要对现有的图切剖面轮廓线作额外的拆分、简化、编辑修改等工作,只需要对地层进行填充和属性赋值即可。
(3)采用了基于线性队列的快速泛洪算法,创建了基于线程的核函数,以多线程的方式对地层边界进行填充,提高了填充效率。
3、本发明建模精度更高,主要表现在:
(1)由于填充过程受轮廓线的显示分辨率影响,图形显示分辨率越高,则追踪的边界点越接近真实的地层边界。而本发明在填充及地层边界检测过程中,对边界点和地层轮廓线进行了投影处理,将追踪的地层边界点自动转换成了轮廓线边界点,以实现地层追踪边界点与原始地层轮廓线的一致性,减少地层边界点计算误差,从而提高模型边界检测的准确性,解决了边缘检测对于显示像素分辨率依赖的问题。
(2)由于采用了边界采样和网格采样两种采样方式结合,既保障了地层原始模型的精度,又不至于产生过大的采样数据量,可以运用三维全局插值算法生成三维属性数据场,并结合Marching Cubes算法从地层三维属性数据场中提取等值面作为地层的三维表面轮廓,充分保障了后续地层三维轮廓表面拼接的精度,建模结果更为合理和准确。
综上,本发明基于三维地层属性场的思想,改进了地层边界追踪算法,提高了地层边界追踪的精度;采用特殊的剖面采样方式,既保障了地层边界采样的精度,又不会产生过大的采样数据量,为三维全局插值算法的应用提供了保障;采用径向基函数和基于GPU的多线程计算技术,保证了地层模型的空间连续性和光滑性,计算效率高。
附图说明
图1为本发明的流程示意图。
图2为本发明-实施例中图切剖面示意图。
图3为本发明-实施例中图切剖面直剖面与折剖面示意图。
图4为本发明-实施例中地层填充示意图。
图5为本发明-实施例中采用泛洪填充算法追踪地层边界点结果示意图。
图6为本发明-实施例中地层边界点投影及冗余点剔除结果示意图。
图7为本发明-实施例中地层边界追踪、边界点投影与冗余点剔除结果示意图。
图8为本发明-实施例中图切剖面平面坐标和空间坐标关系示意图。
图9为本发明-实施例中地层边界及网格离散采样结果示意图。
图10为本发明-实施例中图切剖面三维空间地层离散采样结果示意图。
图11为本发明-实施例中径向基函数三维插值生成地层属性数据场结果示意图。
图12为本发明-实施例中属性值为1的地层边界三维网格示意图。
图13为本发明-实施例中采用Marching Cubes算法计算地层三维轮廓表面的结果示意图。
图14为本发明-实施例中拼接后的三维地质模型结果示意图。
具体实施方式
本发明提供了一种基于地质图切剖面的三维地层建模方法,主要利用地质图切剖面或者其他二维轮廓线进行三维地质模型建模,其建模过程如图1所示,主要包括:
S1、对图切剖面进行编辑,完成建模所需要的基本参数的设定;对图切剖面进行编辑的流程包括:(1)剖面空间定位;(2)地层属性填充;(3)地层边界追踪;其中:通过标定空间定位点将剖面平面坐标转换为地质空间三维坐标;通过地层属性填充获得地层参数信息;通过地层边界追踪获得地层的轮廓边界;
S2、对地层边界及内外区域进行离散采样,并分别设定地层内部、地层边界以及地层外各自的采样点属性值,然后将所有图切剖面采样点转换成空间三维坐标点,组合成空间离散点集合;
S3、对三维空间离散点进行插值计算,生成三维空间地层属性数据场;
S4、从三维地层属性数据场中进行等值面提取,获得地层的三维轮廓表面;
S5、重复步骤S2~S4,获得所有地层的三维轮廓表面,然后执行步骤S6;
S6、将所有地层三维轮廓表面进行空间拼接,完成三维地层建模。
下面结合附图说明和实施例对本发明作进一步说明,本发明的方式包括但不仅限于以下实施例。
实施例
一、对图切剖面进行编辑
图2示出了一种图切剖面示意图,支持矢量图形格式(如DXF)或栅格图形格式(JPG,GIF,BMP等)。图切剖面编辑包括以下步骤:
(1)载入DXF图切剖面图形,在图切剖面上选定四个角点(至少需要2个定位点),并设定4个角点的三维空间坐标,将图切剖面映射到空间坐标系统中。本发明支持直剖面或折剖面的图切剖面。直剖面表示在制作图切剖面时采用直线制作图切剖面,折剖面表示制作图切剖面是采用折线制作图切剖面。折剖面对应的空间定位点应标记在折线的转折点上,见图3;
(2)选定一个地层,采用泛洪填充算法对地层进行填充,并记录其地层边界点,得到地层多边形,并设定地层的名称及属性值(地层唯一识别数字),见图4;
(3)由于泛洪填充算法追踪得到的地层边界点和实际地层轮廓边界存在误差(如图5所示),所以需对地层边界点进行投影处理,以得到准确的地层边界,投影处理步骤如下:
a、将地层轮廓线进行简化,删掉轮廓直线上的冗余点。然后将每两个点组成轮廓线段{Lk}(k=1,2,..N),轮廓线段数为N;
b、选定一个追踪得到的地层边界点Pi(i=l,2,...,M),计算该点到所有轮廓线段的距离集合{Dik};
c、从距离集合{Dik}中找出最小距离的线段Lk,将边界点Pi投影到线段Lk上,得到Pik,将Pi替换为Pik作为新的地层边界点;
d、重复b、c步骤,将所有地层边界点投影到地层轮廓线上;
e、对投影后的地层边界轮廓线进行简化,删除线上冗余点,完成地层边界追踪(如图6、7所示)。
二、离散采样及组合成空间离散点集合
(1)设定一个采用间隔step,沿地层轮廓线边界进行边界点采样,得到地层边界点集合B1,并对集合中所有点属性值设为1;
(2)设定一个指定大小的剖分网格M×N(例如M=100,N=80),将图切剖面进行二维网格剖分,并记录每个网格中心点坐标,得到剖分点集合B2;
(3)对集合B2中的点进行地层内外判别,若点在地层多边形内,则该点的属性值为1,否则为0;
(4)对集合B1、B2中的所有点进行空间坐标转换,将平面坐标点转换成地质空间坐标点;转换过程如下(如图8所示):
其中,P(x,y)为待转换平面点坐标,P′(x,y,z)为转换后的空间坐标;P1,P2为点P所在的定位点区间平面坐标,V1(x,y,z),V2(x,y,z)为P1,P2对应的空间坐标;Y1是定位点在平面Y方向上的最小值,Y2是定位点在平面Y方向上的最大值;Zmin和Zmax为对应的空间上的最小高程和最大高程值;
平面点P(x,y)坐标转换为空间坐标P′(x,y,z)的计算公式如下:
p′x=V1x+(px-p1x)*(V2x-V1x)/(p2x-p1x)
p′y=V1y+(px-p1x)*(V2y-V1y)/(p2x-p1x)
p′z=Zmin+(Zmax-Zmin)*(py-Y1)/(Y2-Y1)
图切剖面采样后的采样点分布如图9所示;
(6)将转换后的空间坐标及其属性值(0或1)放入散乱点数据集合S,得到地层采样后的离散空间点数据(如图10所示)。
三、生成三维空间地层属性数据场
采用径向基函数插值算法进行三维离散点空间插值,生成三维空间地层属性数据场。径向基函数形式如下:
其中Z0是待估测的插值点值,wi是相应点的权值,||pi-p0||是pi与p0之间的距离。基函数采用逆多元二次函数,其函数形式如下:
插值计算过程如下:
(2)采用高斯消元法结合OpenCL多线程加速技术,计算逆距离矩阵D-1;
完成三维插值后的三维属性数据场如图11所示,其中属性值为1的网格分布如图12所示。
四、获得地层的三维轮廓表面
采用MarchingCubes算法从地层三维属性数据场中提取属性值为l的等值面作为地层的三维表面轮廓,完成该地层的三维建模,如图13所示。
MarchingCubes算法提取等值面的过程如下:
(1)从三维网格数据中选取相邻8个点组成六面体;
(2)给定等值面值V0,对六面体中的每个节点状态值Vi进行判别,若Vi>V0,该节点状态设为1;否则设为0;i=1,2,...,8;
(3)根据节点状态序列,按照顺序组合成二进制字节,得到六面体单元构型表索引值Ti(0≤Ti≤255);
(4)根据构型索引值Ti,从构型表中检索出对应的棱边节点状态组合,计算棱边节点状态为0和1的等值面在棱边的交点坐标P(x,y,z),计算公式为:
p(x,y,z)=p1(x,y,z)+(v-v1)/(v2-v1)*(p2(x,y,z)-p1(x,y,z))
式中,p(x,y,z)是待求交点坐标,v是待计算等值面值;p1,p2是棱边两节点坐标,v1,v2是两节点属性值(坐标满足矢量运算法则);本实施例采用了HashMap线性表对棱边交点进行存储,确保等值面与每条棱边交点不被重复计算和重复存储;
(5)根据构型表索引值Ti,从构型表中检索出对应的三角形顶点及连接顺序,组合成封闭的三角形面;
(6)遍历所有网格单元,重复步骤(1)~(5),得到目标地层表面轮廓的所有三角面。
五、重复上述流程,创建所有地层的三维表面模型。
六、完成三维地层建模
将所有地层模型在空间上进行拼接,拼接过程为:
(1)将所有顶点进行组合,并进行重复点剔除,得到新的顶点集合;
(2)根据新的顶点集合,修改所有模型的三角面片索引。
如此即可生成一个三维模型,并输出为三维图形格式(如obj或ply),从而完成三维地层建模,如图14所示。
本发明设计合理、严谨,既避免了复杂的轮廓对应、三角网连接等问题,又可以得到较高精度的三维矢量地质模型。本发明很好地突破了现有技术的限制,顺应科技发展的潮流,很好地匹配了当前阶段工程上对三维地质建模的需求。因此,与现有技术相比,本发明技术进步十分明显,具有突出的实质性特点和显著的进步。
上述实施例仅为本发明的优选实施方式之一,不应当用于限制本发明的保护范围,凡在本发明的主体设计思想和精神上作出的毫无实质意义的改动或润色,其所解决的技术问题仍然与本发明一致的,均应当包含在本发明的保护范围之内。
Claims (3)
1.一种基于地质图切剖面的三维地层建模方法,其特征在于,包括以下步骤:
S1、对图切剖面进行编辑,完成建模所需要的基本参数的设定;对图切剖面进行编辑的流程包括:(1)通过标定空间定位点将剖面平面坐标转换为地质空间三维坐标;(2)通过地层属性填充获得地层参数信息;(3)通过地层边界追踪获得地层的轮廓边界;该步骤中,对于追踪得到地层边界点,采用了投影处理,以消除由于栅格图像显示分辨率影响带来的边界点追踪误差,确保追踪得到的地层轮廓边界与原始地层轮廓边界的一致性,其中,所述投影处理包括以下步骤:
(a)将地层轮廓线进行简化,删除轮廓直线上的冗余点,然后将每两个点组成轮廓线段{Lk}(k=1,2,..N),轮廓线段数为N;
(b)选定一个追踪得到的地层边界点Pi(i=1,2,…,M),计算该点到所有轮廓线段的距离集合{Dik}(k=1,2,..N);
(c)从距离集合{Dik}中找出最小距离的线段Lj,将边界点Pi投影到线段Lj上,得到Pij,将Pi替换为Pij作为新的地层边界点;
(d)重复步骤(b)、(c),将所有地层边界点投影到地层轮廓线上;
(e)对投影后的地层边界轮廓线进行简化,删除直线上的冗余点,完成地层边界追踪;
S2、对地层边界及内外区域进行离散采样,并分别设定地层内部、地层边界以及地层外各自的采样点属性值,然后将所有图切剖面采样点转换成空间三维坐标点,组合成空间离散点集合;该步骤包括以下步骤:
S2-1、设定一个采样间隔step,并沿地层轮廓线边界进行边界点采样,得到地层边界点集合B1,并对集合中所有点属性值设为1;
S2-2、设定一个指定大小的剖分网格M×N,将图切剖面进行二维网格剖分,并记录每个网格中心点坐标,得到剖点集合B2;
S2-3、对集合B2中的点进行地层内外判别,若点在地层多边形内,则该点的属性值设置为1,否则置为0;
S2-4、对集合B1、B2中的所有点进行空间坐标转换,将平面坐标点转换成地质空间坐标点;
S2-5、将转换后的空间坐标点及其属性值放入散乱点数据集合{pi}(i=1,2,…,N),,得到地层采样后的离散空间点数据,完成采样;
S3、对三维空间离散点进行插值计算,生成三维空间地层属性数据场;该步骤中,以逆多元二次径向基函数作为插值函数,采用径向基函数插值算法,进行三维离散点空间插值,生成三维空间地层属性数据场;径向基函数形式如下:
插值计算过程如下:
(g)采用高斯消元法结合OpenCL多线程加速技术,计算逆距离矩阵D-1;
(i)计算属性场值,三维空间中任意位置的场值S4、从三维地层属性数据场中进行等值面提取,获得地层的三维轮廓表面;该步骤中,采用Marching Cubes算法从地层三维属性数据场中提取等值面作为地层的三维表面轮廓,其过程如下:
(j)从三维网格数据中选取相邻8个点组成六面体;
(k)给定等值面值V0,对六面体中的每个节点状态值Vi进行判别,若Vi>V0,该节点状态设为1;否则设为0;i=1,2,…,8;
(l)根据节点状态序列,按照顺序组合成二进制字节,得到六面体单元构型表索引值Ti(0≤Ti≤255);
(m)根据构型索引值Ti,从构型表中检索出对应的棱边节点状态组合,计算棱边节点状态为0和1的等值面在棱边的交点坐标P(x,y,z),计算公式为:
p(x,y,z)=p1(x,y,z)+(v-v1)/(v2-v1)*(p2(x,y,z)-p1(x,y,z))
式中,p(x,y,z)是待求交点坐标,v是待计算等值面值;p1、p2是棱边两节点坐标,v1、v2是两节点属性值;
(n)根据构型表索引值Ti,从构型表中检索出对应的三角形顶点及连接顺序,组合成封闭的三角形面;
(o)遍历所有网格单元,重复步骤(j)~(n),得到目标地层表面轮廓的所有三角面;
S5、重复步骤S2~S4,获得所有地层的三维轮廓表面,然后执行步骤S6;
S6、将所有地层三维轮廓表面进行空间拼接,完成三维地层建模。
2.根据权利要求1所述的一种基于地质图切剖面的三维地层建模方法,其特征在于,所述步骤S1中,采用基于线性队列的快速泛洪填充算法对地层进行填充,并记录其边界点,得到地层多边形,获得地层参数信息。
3.根据权利要求1所述的一种基于地质图切剖面的三维地层建模方法,其特征在于,所述步骤S6包括以下步骤:
S6-1、将所有顶点进行组合,并进行重复点剔除,得到新的顶点集合;
S6-2、根据新的顶点集合,修改所有模型的三角面片索引,完成三维地层建模。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010819474.4A CN111968231B (zh) | 2020-08-14 | 2020-08-14 | 一种基于地质图切剖面的三维地层建模方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010819474.4A CN111968231B (zh) | 2020-08-14 | 2020-08-14 | 一种基于地质图切剖面的三维地层建模方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111968231A CN111968231A (zh) | 2020-11-20 |
CN111968231B true CN111968231B (zh) | 2023-05-30 |
Family
ID=73366107
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010819474.4A Active CN111968231B (zh) | 2020-08-14 | 2020-08-14 | 一种基于地质图切剖面的三维地层建模方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111968231B (zh) |
Families Citing this family (13)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113268900B (zh) * | 2021-04-02 | 2022-09-16 | 中国人民解放军战略支援部队信息工程大学 | 面向任务的空降场选址方法及装置 |
CN113192207B (zh) * | 2021-04-29 | 2024-06-14 | 西安恒歌数码科技有限责任公司 | 基于osg的物体外轮廓表面重建方法及系统 |
CN113593030B (zh) * | 2021-08-13 | 2024-04-19 | 长春工程学院 | 基于三维模型的地质剖面图生成方法、系统、终端及介质 |
CN114255188B (zh) * | 2021-12-23 | 2022-06-24 | 中国矿业大学(北京) | 三维地质岩性网格模型体绘制的边界线性平滑方法及装置 |
CN114612318B (zh) * | 2022-02-16 | 2024-04-19 | 西北大学 | 基于文物ct图像轮廓线的三维建模方法、系统及设备 |
CN114635681B (zh) * | 2022-03-22 | 2022-12-06 | 成都理工大学 | 一种高砂地比厚层辫状河三角洲前缘砂体构型构建方法 |
CN114677481B (zh) * | 2022-05-31 | 2022-09-13 | 中国飞机强度研究所 | 空天飞机地面测试的理想加热曲面等效逼近模型构建方法 |
CN115719411B (zh) * | 2023-01-10 | 2023-04-18 | 东华理工大学南昌校区 | 三维地质建模方法、系统、计算机及可读存储介质 |
CN117392335B (zh) * | 2023-09-26 | 2024-08-20 | 深圳市地质环境研究院有限公司 | 三维面元地质模型向三维体元地质模型的转换方法 |
CN117115391B (zh) * | 2023-10-24 | 2024-01-12 | 中科云谷科技有限公司 | 模型更新方法、装置、计算机设备及计算机可读存储介质 |
CN117541741B (zh) * | 2024-01-10 | 2024-03-19 | 自然资源部第三地理信息制图院 | 一种地质体三维模型的构建方法及电子设备 |
CN118332812A (zh) * | 2024-04-29 | 2024-07-12 | 浙江大学 | 一种优化二维块体识别结果的小块体处理方法 |
CN118505921B (zh) * | 2024-07-17 | 2024-10-11 | 山东省地矿工程勘察院(山东省地质矿产勘查开发局八〇一水文地质工程地质大队) | 一种地质档案的模型三维可视化方法及其系统 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109326002A (zh) * | 2018-11-27 | 2019-02-12 | 中南大学 | 基于钻孔数据的矿体建模方法、装置、系统及存储介质 |
CN109584357A (zh) * | 2018-11-27 | 2019-04-05 | 中南大学 | 基于多轮廓线的三维建模方法、装置、系统及存储介质 |
Family Cites Families (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2005010797A2 (en) * | 2003-07-23 | 2005-02-03 | Lee Wook B | Improved 3d veloctiy modeling, with calibration and trend fitting using geostatistical techniques, particularly advantageous for curved-ray prestack time migration and for such migration followed by prestack depth migration |
CN104715506B (zh) * | 2015-04-01 | 2017-08-25 | 中国地质大学(北京) | 双体式地质模型的构建方法 |
CN106709987B (zh) * | 2015-11-13 | 2020-01-17 | 星际空间(天津)科技发展有限公司 | 一种三维地质剖面模型动态构建方法 |
CN106097448A (zh) * | 2016-06-08 | 2016-11-09 | 江西理工大学 | 一种多特征约束下的盐腔围岩地质三维建模方法 |
CN109753707B (zh) * | 2018-12-25 | 2023-10-24 | 核工业北京地质研究院 | 一种利用勘探线剖面提取地层界线开展三维建模的方法 |
CN110599594B (zh) * | 2019-07-29 | 2021-07-20 | 成都理工大学 | 一种岩石物性结构三维建模的方法 |
CN111415415B (zh) * | 2020-03-31 | 2023-05-09 | 南京师范大学 | 一种三维图切地质剖面的自动构建方法 |
-
2020
- 2020-08-14 CN CN202010819474.4A patent/CN111968231B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109326002A (zh) * | 2018-11-27 | 2019-02-12 | 中南大学 | 基于钻孔数据的矿体建模方法、装置、系统及存储介质 |
CN109584357A (zh) * | 2018-11-27 | 2019-04-05 | 中南大学 | 基于多轮廓线的三维建模方法、装置、系统及存储介质 |
Also Published As
Publication number | Publication date |
---|---|
CN111968231A (zh) | 2020-11-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111968231B (zh) | 一种基于地质图切剖面的三维地层建模方法 | |
Gong et al. | Three-dimensional modeling and application in geological exploration engineering | |
CN106934860B (zh) | 一种基于t样条的三维地质建模方法 | |
CN102194253B (zh) | 一种面向三维地质层面结构的四面体网格生成方法 | |
Pellerin et al. | Automatic surface remeshing of 3D structural models at specified resolution: A method based on Voronoi diagrams | |
CN104966317B (zh) | 一种基于矿体轮廓线的三维自动建模方法 | |
Chen et al. | Uniform offsetting of polygonal model based on layered depth-normal images | |
US20120166160A1 (en) | Block model constructing method for complex geological structures | |
CN111797555B (zh) | 一种基于有限元模型的几何重构方法 | |
CN103903061B (zh) | 三维矿产资源预测评价中信息综合处理装置及其方法 | |
CN109658431B (zh) | 基于区域生长的岩体点云平面提取方法 | |
CN106844871A (zh) | 基于bim的三维岩土工程勘察信息模型构建方法 | |
CN110322547B (zh) | 一种储层自适应四面体剖分方法 | |
Allili et al. | Cubical homology and the topological classification of 2D and 3D imagery | |
CN102867332A (zh) | 基于复杂边界约束的多级细分网格曲面拟合方法 | |
CN113674388A (zh) | 一种基于机器学习的三维地质体纹理贴图方法 | |
CN113971718B (zh) | 一种对三维点云模型进行布尔运算的方法 | |
CN113658333B (zh) | 一种基于等值面抽取的地质体建模方法 | |
de Oliveira Miranda et al. | Finite element mesh generation for subsurface simulation models | |
CN117496077A (zh) | 一种三维地质模型可视化分析方法及系统 | |
CN117173357A (zh) | 一种矿山三维地质建模与分层切割方法 | |
Chen et al. | A point-based offsetting method of polygonal meshes | |
Pellerin et al. | Topological control for isotropic remeshing of non-manifold surfaces with varying resolution: application to 3D structural models | |
CN106097447A (zh) | 一种大规模地震数据的曲面重建方法 | |
Zhang et al. | An automatic unified modeling method of geological object and engineering object based on tri-prism (TP) |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |