一种顾及褶皱的三维地质结构模型角点网格剖分方法与流程

本发明涉及一种顾及褶皱的三维地质结构模型角点网格剖分方法，属于三维图形学及三维地质建模领域。

背景技术：
三维地质建模技术是地质信息科技领域的核心技术，也是三维空间信息系统建设的关键技术和热点问题。地质空间数据具有多源、异构、多尺度等特征，这对具有真三维数据结构的三维地质模型的构建和可视化提出了更高的要求，所建模型不仅要满足复杂地质结构的三维可视化展示要求，更应支持后续地质空间分析。三维地质建模可分为三维地质结构建模和地质属性建模。三维地质结构建模是指包含地层、岩体、地质构造等地质结构的空间形态建模；地质属性建模是指包括密度、电阻率、矿化度、孔隙度等各种物化探参数在内的地质属性场的空间分布特征建模。长期以来，人们更加关注三维地质结构模型的构建方法，特别对尖灭、断层、褶皱等复杂地质结构模型的构建方法都有较深入的研究。但是，对于地质体属性建模的理论和方法还不是很成熟，没有统一的数据模型来支持复杂地质结构内部属性的空间分布及表达，难以实现地质结构模型与精细体元模型之间的转换方法。体元模型基于三维空间的体元分割对地质体进行三维实体表达，有利于三维空间的地质属性的可视化表达与分析。按体元的面数，体元模型可以分为四面体(Tetrahedron)、五面体(Pentrahedron)、六面体(Hexahedron)、多面体(Polyhedron)和棱柱体(Prism)等5种基本类型，也可以根据体元的规整性分为规则体元和非规则体元两个大类。角点网格模型是体元模型的一个代表性数据模型，是目前应用较广的一种结构化网格类型，本质上是不规则六面体网格，网格位置能用行列层的编号i、j、k来定义，单元网格的长、宽大小可变，垂向连接顶底网格点的网格面也可以是倾斜的，因此可以很好的表达出复杂地质构造的形态。空间体元剖分是地质属性建模的基础，只有对地质空间按照地质构造规律进行合理的体元剖分，才能在此基础上利用空间插值等方法对每个空间体元赋予不同的地质属性，从而实现地质结构内部非均质属性特征的三维可视化表达。然而现有的空间体元剖分方法没有考虑地质结构的空间分布特征，没有以地质结构模型为依据进行空间体元的划分，这使得所建立的属性模型缺乏地质构造条件的约束，为后续的属性赋值以及其对真实地质属性表达的准确性带来了困难。褶皱是地质结构中最基本的构造形式，普遍存在且结构复杂多样。褶皱形态多种多样，根据轴面产状和两翼产状特点可分为：直立褶皱、斜歪褶皱、倒转褶皱和平卧褶皱(如图2所示)。兰姆赛的经典专著“现代构造地质学方法(第二卷)”中对褶皱的特征及分类给出了详细的解释(Ramsay，1967)，该分类法已经成为描述褶皱的标准方法。如图3(a)所示，是兰姆赛给出的确定等倾斜线和正交厚度的方法，tα为褶皱等倾角线之间的距离，即正交厚度。以t′α＝cosα为标准化正交厚度，据此褶皱的基本类型可描述为：1A类褶皱翼部的正交厚度大于枢纽处的正交厚度；1B类褶皱的正交厚度处处相等，为平行褶皱；1C、2和3类褶皱的翼部都变薄，与真正的相似褶皱(2类)相比，1C类翼部减薄少一些，而3类多一些(如图3(b))。如图3(c)所示为5种基本褶皱类型的倾角等斜线，通过这种图解方法能够迅速的确定岩层的基本褶皱类型，并能反映褶皱形成的构造应力情况和岩体的基本特征。因此，在对包含褶皱的三维地质结构模型进行空间体元剖分时，一定要考虑褶皱的空间展布形态，使其符合兰姆赛褶皱等斜线分类的基本要求，而不是做简单的空间体元划分。鉴于上述分析的顾及褶皱的复杂地质结构的特征，以及现有体元模型的研究现状，有必要提出一种专门针对褶皱的三维地质结构模型的体元剖分方法。

技术实现要素：
为了解决现有技术的不足，本发明提供了一种顾及褶皱的三维地质结构模型角点网格剖分方法，面向包含褶皱的复杂地质结构特征，及其在地质属性模型和后续地质空间分析中的应用需求，以角点网格作为基本数据结构模型，根据兰姆赛褶皱分类中等斜线的规定，实现了三维地质结构模型到精细体元模型的快速转换，而且其体元形态符合褶皱、尖灭等地质构造的空间特征，为后续空间体元的属性赋值提供了载体和模型基础。本发明为解决其技术问题所采用的技术方案是：提供了一种顾及褶皱的三维地质结构模型角点网格剖分方法，包括以下步骤：首先按空间轴向顺序依次获取地质结构模型的各层地层体模型，对于每层地层体模型，执行步骤(1)至步骤(9)；(1)确定平面网格：设置横向网格数目m和纵向网格数目n，根据地质结构模型在XOY平面上的投影的最大边界矩形，根据横向网格数目m、纵向网格数目n以及地质结构模型的趋势线对最大边界矩形进行平面网格划分，使划分的平面网格的横线和纵线分别与地质结构模型的趋势线的方向和趋势线的垂直方向一致，相邻两条趋势线之间的网格数目相同；(2)确定地层体模型的空间轴向：对于第i层地层体模型PHi，通过以下公式计算其上表面或下表面的平均法线方向其中，S为PHi曲面的表面积，即组成PHi曲面的所有三角面面积之和，sk为第k个三角面的面积，N为三角面的数目，表示PHi曲面的所有三角面中第k个三角面的法线方向的单位向量，则最终求得的向量所表示的方向即为PHi的空间轴向；(3)确定旋转变换矩阵：运用右手系坐标，将PHi的空间轴向表示为将通过旋转变换矩阵M旋转至YOZ平面，所述旋转变换矩阵M为：其中，(4)地层体模型旋转：按照旋转变换矩阵M，将PHi上的所有数据点进行旋转，则对于PHi上的任一点P(x,y,z)，旋转后的点P′(x′,y′,z′)表示为(x′,y′,z′)＝(x,y,z)M；PHi旋转后为PHi′；(5)坐标投影：对m×n个平面网格点，依次进行垂直坐标投影，得到所有平面网格点与PHi′的上表面或下表面的交点，即投影点，如果没有交点，则该投影点的Z坐标设置为无效值，该投影点为无效的投影点，该投影点对应的原平面网格点无效；(6)网格有效性判定：对m×n个平面网格点，依次取相邻4个平面网格点，如果4个平面网格点均为有效值，则说明4个平面网格点构成的网格是有效的，否则为无效网格；(7)确定各投影点位置的倾角等斜线：以行优先开始计算，过每一行第一个投影点pu,1做上表面的垂线，其中u∈[1,m]，求取与下表面的交点p′u,1，初始化pu,1与p′u,1的连线为投影点pu,1位置的倾角等斜线，对于其它任意一个有效的投影点pu,v，v∈(1,n]，执行以下过程：确定投影点pu,v处的法线方向，获取投影点pu,v-1位置的倾角等斜线与地层体模型PHi′的下表面的交点p′u,v-1，求取投影点pu,v和交点p′u,v-1的连线与投影点处法线的夹角α，从投影点pu,v和交点p′u,v-1的连线处开始，以设置步长的旋转角ω求取与地层体模型PHi′下表面的交点，总的旋转角度ω不超过2α，选取交点处法线与投影点pu,v处法线最相近且距离最短的交点p′u,v，投影点pu,v与该交点p′u,v的连线为该投影点pu,v位置的倾角等斜线；(8)确定垂向网格：计算投影点pu,v与交点p′u,v之间的距离D，设PHi′的细分层数为t，则投影点pu,v处的垂向网格划分时每一步的生长距离d＝D/t，然后从投影点pu,v开始沿着直线pu,vp′u,v方向以距离d进行剖分，并依次记录新得到的点的空间坐标；对所有有效投影点执行完毕得到PHi′垂向网格划分后，按顺序存储的所有空间点坐标；(9)逆向旋转变换：该PHi′的所有垂向网格划分后的按顺序存储的所有空间点进行逆矩阵变换，以还原地层的真实形态和网格点的真实坐标，其逆向旋转的变换矩阵表示为：各层地层体模型执行完步骤(2)至步骤(9)后，按照空间坐标及拓扑序列，构建空间不规则六面体体元，得到对应于地质结构模型的精细体元模型。步骤(2)中对于上表面或下表面的选取，是由剖分顺序来定，如果是按照从上到下的地层序列，则选取地层体模型的上表面，反之则为下表面。步骤(6)中网格的有效性决定了地层PHi′在对应网格处的有效性，即对于同一套地层，细分后每一个小层的有效性是一致的，所以对某一套地层的剖分只需要计算一次。步骤(7)中投影点pu,v处法线的计算，具体过程为：以该投影点pu,v为中心，取邻接的四个平面网格点，所述四个平面网格点均未超出网格边界范围并且有效；以投影点pu,v为中心将四个平面网格点组成的四边形划分为四个三角面，计算四个三角面的平均法线方向，确定投影点pu,v的法线方向的单位向量本发明基于其技术方案所具有的有益效果在于：(1)本发明充分考虑了褶皱的空间形态特征，采用角点网格为基本数据模型，通过平面趋势线和垂向细划分数目的控制，实现了空间体元的局部加密，可对需详细研究区域进行更加细致的剖分；(2)本发明以地层结构模型为层控因素，使得剖分后的空间体元依然具备结构模型的特征，实现了从地质结构模型到精细体元模型的转换；(3)本发明通过空间矩阵变换的方式，将因褶皱而产生的多值曲面转换为单值曲面，进行空间体元剖分后在通过逆矩阵还原，有效解决了褶皱剖分中因投影的多值性带来的困难；(4)本发明采用沿着褶皱面倾角等斜线方向生长的剖分方式，不仅使剖分得到的体元的整体形态与褶皱形态一致，满足兰姆赛褶皱倾角等斜线的要求，而且使得剖分后同一层的体元具有更加相似的成因，即剖分得到的体元分层与地质沉积过程相一致；(5)本发明可在各种三维地质建模、地质信息系统、地理信息系统、地质过程模拟系统等软件中推广应用。附图说明图1是本发明的顾及褶皱的三维地质结构模型角点网格剖分方法流程示意图。图2是褶皱基本形态示意图。图3是褶皱基本类型及倾角等斜线特征，其中(a)表示等斜线及正交厚度的确定方法，(b)表示标准化正交厚度t′α对倾角α的曲线图及褶皱分类的基本类型，(c)表示基本褶皱类型的倾角等斜线特征。图4是本发明中平面网格剖分示意图。图5是本发明中地层体模型空间轴向旋转过程示意图。图6是本发明中投影点pu,v倾角等斜线的确定方法处示意图。图7是本发明中顾及褶皱的地质结构模型垂向体元剖分示意图，其中(a)为某地质结构模型的垂向切面示意图，(b)为第i套地层PHi的轴向旋转示意图，(c)为旋转后的地层的垂向网格剖分示意图，(d)为轴向复位后地层的垂向网格剖分状态及其与上下地层的接触关系。图8是利用本发明方法对包含褶皱的复杂地质结构模型的角点网格剖分效果示意图，其中(a)为含褶皱的地层结构模型，(b)为整体剖分效果网格显示，(c)为褶皱层的剖分效果，(d)为剖分后的一个小层。具体实施方式下面结合附图和实施例对本发明作进一步说明。本发明提供了一种顾及褶皱的三维地质结构模型角点网格剖分方法，在三维地质结构模型的基础上，充分考虑褶皱的空间形态特征，以兰姆赛褶皱基本类型和倾角等斜线为依据，以角点网格为空间体元数据模型，以地层结构体模型为层向控制，通过趋势线和地质结构模型在XOY平面的投影范围，进行平面网格划分，然后计算地层体模型的空间轴向并通过矩阵变换将其旋转到与Z轴一致，垂向网格划分则是由曲面上投影点沿该点倾角等斜线方向按照一定步长生长得到，通过平面和垂向网格划分，计算得到对应点的空间坐标，依据坐标和拓扑信息构建角点网格模型并可视化显示，实现了三维地质结构模型到精细体元模型的快速转换，而且其体元形态符合褶皱、尖灭等复杂地质构造的空间特征，为后续空间体元的属性赋值提供了载体和模型基础。参照图1，顾及褶皱的三维地质结构模型角点网格剖分方法包括以下步骤：首先输入待剖分地质结构模型PH，其由N层地层体模型组成PH＝{PHi}(1≤PHi≤N)，并将地质结构模型按地层顺序依次放入不同图层，然后根据地层层序编号和对应图层名称定制地层序列定义表，如表1所示；表1地层序列定义表说明对于每层地层体模型，执行步骤(1)至步骤(9)；(1)确定平面网格：设置横向网格数目m和纵向网格数目n，根据地质结构模型在XOY平面上的投影的最大边界矩形，根据横向网格数目m、纵向网格数目n以及地质结构模型的趋势线对最大边界矩形进行平面网格划分，使划分的平面网格的横线和纵线分别与地质结构模型的趋势线的方向和趋势线的垂直方向一致，相邻两条趋势线之间的网格数目相同；具体地，沿着模型走向方向绘制若干条趋势线，计算趋势线和地质结构模型的最大边界矩形，即求得左下角坐标Pmin(xmin,ymin)和右上角坐标Pmax(xmax,ymax)；确定趋势线的平均方向保持为平面网格的U方向(横向)，与趋势线的平均方向垂直的方向为V方向(纵向)；输入横纵网格数目(m,n)，以Pmin(xmin,ymin)和Pmax(xmax,ymax)确定的矩形范围为边界，将U方向划分为m个网格，V方向划分为n个网格，并保证趋势线处网格划分必须与趋势线走线一致，相邻两条趋势线之间网格数目相同，其平面网格划分如图4所示，其中l1、l2、l3、l4、l5为趋势线，L为模型范围边界；(2)确定地层体模型的空间轴向：对于第i层地层体模型PHi，通过以下公式计算其上表面或下表面的平均法线方向其中，S为PHi曲面的表面积，即组成PHi曲面的所有三角面面积之和，sk为第k个三角面的面积，N为三角面的数目，表示PHi曲面的所有三角面中第k个三角面的法线方向的单位向量，则最终求得的向量所表示的方向即为PHi的空间轴向；(3)确定旋转变换矩阵：如图5所示，运用右手系坐标，将PHi的空间轴向表示为将通过旋转变换矩阵M旋转至YOZ平面，所述旋转变换矩阵M为：其中，(4)地层体模型旋转：按照旋转变换矩阵M，将PHi上的所有数据点进行旋转，则对于PHi上的任一点P(x,y,z)，旋转后的点P′(x′,y′,z′)表示为(x′,y′,z′)＝(x,y,z)M；PHi旋转后为PHi′；(5)坐标投影：对m×n个平面网格点，依次进行垂直坐标投影，得到所有平面网格点与PHi′的上表面或下表面的交点，即投影点，得到交点坐标集pointArray[m][n]，如果没有交点，则该投影点的Z坐标设置为无效值，该投影点为无效的投影点，该投影点对应的原平面网格点无效；(6)网格有效性判定：对m×n个平面网格点，依次取相邻4个平面网格点pointArray[i][j]、pointArray[i][j+1]、pointArray[i+1][j]、pointArray[i][j+1]，如果4个平面网格点均为有效值，则说明4个平面网格点构成的网格是有效的，标记数组A[i][j]＝1，否则为无效网格，标记数组A[i][j]＝0；(7)参照图6，确定各投影点位置的倾角等斜线：以行优先开始计算，过每一行第一个投影点pu,1做上表面的垂线，其中u∈[1,m]，求取与下表面的交点p′u,1，初始化pu,1与p′u,1的连线为投影点pu,1位置的倾角等斜线，对于其它任意一个有效的投影点pu,v，v∈(1,n]，执行以下过程：确定投影点pu,v处的法线方向，获取投影点pu,v-1位置的倾角等斜线与地层体模型PHi′的下表面的交点p′u,v-1，求取投影点pu,v和交点p′u,v-1的连线与投影点处法线的夹角α，从投影点pu,v和交点p′u,v-1的连线处开始，以设置步长的旋转角ω求取与地层体模型PHi′下表面的交点，总的旋转角度ω不超过2α，选取交点处法线与投影点pu,v处法线最相近且距离最短的交点p′u,v，投影点pu,v与该交点p′u,v的连线为该投影点pu,v位置的倾角等斜线；(8)确定垂向网格：计算投影点pu,v与交点p′u,v之间的距离D，设PHi′的细分层数为t，则投影点pu,v处的垂向网格划分时每一步的生长距离d＝D/t，然后从投影点pu,v开始沿着直线pu,vp′u,v方向以距离d进行剖分，并依次记录新得到的点的空间坐标；可以设置当前生长步数c进行循环控制，若c≤t则返回继续求取投影点生长后的空间坐标；对所有有效投影点执行完毕得到PHi′垂向网格划分后，按顺序存储的所有空间点坐标；(9)逆向旋转变换：参照图7，该PHi′的所有垂向网格划分后的按顺序存储的所有空间点进行逆矩阵变换，以还原地层的真实形态和网格点的真实坐标，其逆向旋转的变换矩阵表示为：各层地层体模型执行完步骤(2)至步骤(9)后，得到整个地质结构模型的空间体元划分结果。按照空间坐标及拓扑序列，构建空间不规则六面体体元，得到对应于地质结构模型的精细体元模型，如图8所示。步骤(2)中对于上表面或下表面的选取，是由剖分顺序来定，如果是按照从上到下的地层序列，则选取地层体模型的上表面，反之则为下表面。步骤(6)中网格的有效性决定了地层PHi′在对应网格处的有效性，即对于同一套地层，细分后每一个小层的有效性是一致的，所以对某一套地层的剖分只需要计算一次。步骤(7)中投影点pu,v处法线的计算，具体过程为：以该投影点pu,v为中心，取邻接的四个平面网格点，所述四个平面网格点均未超出网格边界范围并且有效；以投影点pu,v为中心将四个平面网格点组成的四边形划分为四个三角面，计算四个三角面的平均法线方向，确定投影点pu,v的法线方向的单位向量以上实施例仅供说明本发明之用，而非对本发明的限制，有关技术领域的技术人员，在不脱离本发明的精神和范围的情况下，还可以作出各种变换或变型，因此所有等同的技术方案，都落入本发明的保护范围。