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

摘要

本发明提供了一种顾及褶皱的三维地质结构模型角点网格剖分方法，在三维地质结构模型的基础上，充分考虑了褶皱的空间形态特征和兰姆赛褶皱分类中等斜线的规定，以角点网格为空间体元数据模型，以地层结构体模型为层向控制，进行地质结构模型的角点网格剖分。通过趋势线和地质结构模型在XOY平面的投影范围，进行平面网格划分。对地质模型所包含褶皱的多值面地层进行空间旋转。垂向网格划分由曲面上投影点沿倾角等倾斜线方向按照一定步长生长得到，剖分完成后经逆矩阵变换复位。通过横纵向网格划分，计算得到对应点的空间坐标，依据坐标和拓扑信息构建空间体元模型并可视化显示，实现了顾及褶皱的三维地质结构模型到角点网格体元模型的转换。

说明书

技术领域

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

背景技术

三维地质建模技术是地质信息科技领域的核心技术，也是三维空间信息系统建设 的关键技术和热点问题。地质空间数据具有多源、异构、多尺度等特征，这对具有真三维数 据结构的三维地质模型的构建和可视化提出了更高的要求，所建模型不仅要满足复杂地质 结构的三维可视化展示要求，更应支持后续地质空间分析。三维地质建模可分为三维地质 结构建模和地质属性建模。三维地质结构建模是指包含地层、岩体、地质构造等地质结构的 空间形态建模；地质属性建模是指包括密度、电阻率、矿化度、孔隙度等各种物化探参数在 内的地质属性场的空间分布特征建模。长期以来，人们更加关注三维地质结构模型的构建 方法，特别对尖灭、断层、褶皱等复杂地质结构模型的构建方法都有较深入的研究。但是，对 于地质体属性建模的理论和方法还不是很成熟，没有统一的数据模型来支持复杂地质结构 内部属性的空间分布及表达，难以实现地质结构模型与精细体元模型之间的转换方法。

体元模型基于三维空间的体元分割对地质体进行三维实体表达，有利于三维空间 的地质属性的可视化表达与分析。按体元的面数，体元模型可以分为四面体 (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，通过以下公式计算其 上表面或下表面的平均法线方向

$\overrightarrow{F}=\underset{1}{\overset{N}{\Sigma }}\frac{s_k}{S}\times {\overrightarrow{f}}_k$

其中，S为PHi曲面的表面积，即组成PHi曲面的所有三角面面积之和，s_k为第k个三 角面的面积，N为三角面的数目，表示PHi曲面的所有三角面中第k个三角面的法线方向的 单位向量，则最终求得的向量所表示的方向即为PHi的空间轴向；

(3)确定旋转变换矩阵：运用右手系坐标，将PHi的空间轴向表示为将通过旋转变换矩阵M旋转至YOZ平面，所述旋转变换矩阵M为：

其中，$\begin{array}{cc}cos\theta =\frac{b}{\sqrt{a^2+b^2}},& sin\theta =\frac{a}{\sqrt{a^2+b^2}}\end{array},$

(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)确定各投影点位置的倾角等斜线：以行优先开始计算，过每一行第一个投影点 p_u,1做上表面的垂线，其中u∈[1,m]，求取与下表面的交点p′_u,1，初始化p_u,1与p′_u,1的连线为 投影点p_u,1位置的倾角等斜线，对于其它任意一个有效的投影点p_u,v，v∈(1,n]，执行以下过 程：确定投影点p_u,v处的法线方向，获取投影点p_u,v-1位置的倾角等斜线与地层体模型PHi′的 下表面的交点p′_u,v-1，求取投影点p_u,v和交点p′_u,v-1的连线与投影点处法线的夹角α，从投影 点p_u,v和交点p′_u,v-1的连线处开始，以设置步长的旋转角ω求取与地层体模型PHi′下表面的 交点，总的旋转角度ω不超过2α，选取交点处法线与投影点p_u,v处法线最相近且距离最短的 交点p′_u,v，投影点p_u,v与该交点p′_u,v的连线为该投影点p_u,v位置的倾角等斜线；

(8)垂向格网划分：计算投影点p_u,v与交点p′_u,v之间的距离D，设PHi′的细分层数为 t，则投影点p_u,v处的垂向格网划分时每一步的生长距离d＝D/t，然后从投影点p_u,v开始沿着 直线p_u,vp′_u,v方向以距离d进行剖分，并依次记录新得到的点的空间坐标；对所有有效投影 点执行完毕得到PHi′垂向格网划分后，按顺序存储的所有空间点坐标；

(9)逆向旋转变换：该PHi′的所有垂向格网划分后的按顺序存储的所有空间点进 行逆矩阵变换，以还原地层的真实形态和网格点的真实坐标，其逆向旋转的变换矩阵表示 为：

各层地层体模型执行完步骤(2)至步骤(9)后，按照空间坐标及拓扑序列，构建空 间不规则六面体体元，得到对应于地质结构模型的精细体元模型。

步骤(2)中对于上表面或下表面的选取，是由剖分顺序来定，如果是按照从上到下 的地层序列，则选取地层体模型的上表面，反之则为下表面。

步骤(6)中格网的有效性决定了地层PHi′在对应格网处的有效性，即对于同一套 地层，细分后每一个小层的有效性是一致的，所以对某一套地层的剖分只需要计算一次。

步骤(7)中投影点p_u,v处法线的计算，具体过程为：以该投影点p_u,v为中心，取邻接 的四个平面网格点，所述四个平面网格点均未超出网格边界范围并且有效；以投影点p_u,v为 中心将四个平面网格点组成的四边形划分为四个三角面，计算四个三角面的平均法线方 向，确定投影点p_u,v的法线方向的单位向量

本发明基于其技术方案所具有的有益效果在于：

(1)本发明充分考虑了褶皱的空间形态特征，采用角点网格为基本数据模型，通过 平面趋势线和垂向细划分数目的控制，实现了空间体元的局部加密，可对需详细研究区域 进行更加细致的剖分；

(2)本发明以地层结构模型为层控因素，使得剖分后的空间体元依然具备结构模 型的特征，实现了从地质结构模型到精细体元模型的转换；

(3)本发明通过空间矩阵变换的方式，将因褶皱而产生的多值曲面转换为单值曲 面，进行空间体元剖分后在通过逆矩阵还原，有效解决了褶皱剖分中因投影的多值性带来 的困难；

(4)本发明采用沿着褶皱面倾角等斜线方向生长的剖分方式，不仅使剖分得到的 体元的整体形态与褶皱形态一致，满足兰姆赛褶皱倾角等斜线的要求，而且使得剖分后同 一层的体元具有更加相似的成因，即剖分得到的体元分层与地质沉积过程相一致；

(5)本发明可在各种三维地质建模、地质信息系统、地理信息系统、地质过程模拟 系统等软件中推广应用。

附图说明

图1是本发明的顾及褶皱的三维地质结构模型角点网格剖分方法流程示意图。

图2是褶皱基本形态示意图。

图3是褶皱基本类型及倾角等斜线特征，其中(a)表示等斜线及正交厚度的确定方 法，(b)表示标准化正交厚度t_α′对倾角α的曲线图及褶皱分类的基本类型，(c)表示基本褶 皱类型的倾角等斜线特征。

图4是本发明中平面格网剖分示意图。

图5是本发明中地层体模型空间轴向旋转过程示意图。

图6是本发明中投影点p_u,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以及地质结构模型 的趋势线对最大边界矩形进行平面网格划分，使划分的平面网格的横线和纵线分别与地质 结构模型的趋势线的方向和趋势线的垂直方向一致，相邻两条趋势线之间的网格数目相 同；

具体地，沿着模型走向方向绘制若干条趋势线，计算趋势线和地质结构模型的最 大边界矩形，即求得左下角坐标P_min(x_min,y_min)和右上角坐标P_max(x_max,y_max)；

确定趋势线的平均方向保持为平面格网的U方向(横向)，与趋势线的平均方向垂 直的方向为V方向(纵向)；

输入横纵网格数目(m,n)，以P_min(x_min,y_min)和P_max(x_max,y_max)确定的矩形范围为边 界，将U方向划分为m个格网，V方向划分为n个格网，并保证趋势线处网格划分必须与趋势线 走线一致，相邻两条趋势线之间网格数目相同，其平面格网划分如图4所示，其中l₁、l₂、l₃、 l₄、l₅为趋势线，L为模型范围边界；

(2)确定地层体模型的空间轴向：对于第i层地层体模型PHi，通过以下公式计算其 上表面或下表面的平均法线方向

$\overrightarrow{F}=\underset{1}{\overset{N}{\Sigma }}\frac{s_k}{S}\times {\overrightarrow{f}}_k$

其中，S为PHi曲面的表面积，即组成PHi曲面的所有三角面面积之和，s_k为第k个三 角面的面积，N为三角面的数目，表示PHi曲面的所有三角面中第k个三角面的法线方向的 单位向量，则最终求得的向量所表示的方向即为PHi的空间轴向；

(3)确定旋转变换矩阵：如图5所示，运用右手系坐标，将PHi的空间轴向表示为 将通过旋转变换矩阵M旋转至YOZ平面，所述旋转变换矩阵M为：

其中，$\left(\begin{array}{cc}cos\theta =\frac{b}{\sqrt{a^2+b^2}},& sin\theta =\frac{a}{\sqrt{a^2+b^2}},\end{array}\right)$

(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，确定各投影点位置的倾角等斜线：以行优先开始计算，过每一行第一 个投影点p_u,1做上表面的垂线，其中u∈[1,m]，求取与下表面的交点p′_u,1，初始化p_u,1与p′_u,1的连线为投影点p_u,1位置的倾角等斜线，对于其它任意一个有效的投影点p_u,v，v∈(1,n]，执 行以下过程：确定投影点p_u,v处的法线方向，获取投影点p_u,v-1位置的倾角等斜线与地层体模 型PHi′的下表面的交点p′_u,v-1，求取投影点p_u,v和交点p′_u,v-1的连线与投影点处法线的夹角 α，从投影点p_u,v和交点p′_u,v-1的连线处开始，以设置步长的旋转角ω求取与地层体模型PHi′ 下表面的交点，总的旋转角度ω不超过2α，选取交点处法线与投影点p_u,v处法线最相近且距 离最短的交点p′_u,v，投影点p_u,v与该交点p′_u,v的连线为该投影点p_u,v位置的倾角等斜线；

(8)垂向格网划分：计算投影点p_u,v与交点p′_u,v之间的距离D，设PHi′的细分层数为 t，则投影点p_u,v处的垂向格网划分时每一步的生长距离d＝D/t，然后从投影点p_u,v开始沿着 直线p_u,vp′_u,v方向以距离d进行剖分，并依次记录新得到的点的空间坐标；

可以设置当前生长步数c进行循环控制，若c≤t则返回继续求取投影点生长后的 空间坐标；

对所有有效投影点执行完毕得到PHi′垂向格网划分后，按顺序存储的所有空间点 坐标；

(9)逆向旋转变换：参照图7，该PHi′的所有垂向格网划分后的按顺序存储的所有 空间点进行逆矩阵变换，以还原地层的真实形态和网格点的真实坐标，其逆向旋转的变换 矩阵表示为：

各层地层体模型执行完步骤(2)至步骤(9)后，得到整个地质结构模型的空间体元 划分结果。按照空间坐标及拓扑序列，构建空间不规则六面体体元，得到对应于地质结构模 型的精细体元模型，如图8所示。

步骤(2)中对于上表面或下表面的选取，是由剖分顺序来定，如果是按照从上到下 的地层序列，则选取地层体模型的上表面，反之则为下表面。

步骤(6)中格网的有效性决定了地层PHi′在对应格网处的有效性，即对于同一套 地层，细分后每一个小层的有效性是一致的，所以对某一套地层的剖分只需要计算一次。

步骤(7)中投影点p_u,v处法线的计算，具体过程为：以该投影点p_u,v为中心，取邻接 的四个平面网格点，所述四个平面网格点均未超出网格边界范围并且有效；以投影点p_u,v为 中心将四个平面网格点组成的四边形划分为四个三角面，计算四个三角面的平均法线方 向，确定投影点p_u,v的法线方向的单位向量

以上实施例仅供说明本发明之用，而非对本发明的限制，有关技术领域的技术人 员，在不脱离本发明的精神和范围的情况下，还可以作出各种变换或变型，因此所有等同的 技术方案，都落入本发明的保护范围。