基于似三棱柱的三维水文地质体建模技术研究

基于似三棱柱体的三维水文地质体建模技术研究目 录引言 .21 孔隙地下水含水层系统概念模型 .31.1 含水层的空间分布特点 .31.2 承压含水层系统建模数据源 .31.3 三维水文地质体模型特点 .32 水文地质钻孔模型构建方法 .42.1 钻孔数据结构设计 .42.2 水文地质钻孔建模方法 .43 分层界面三角网（TIN） .73.1 三角网模型数据结构 .73.2 三角形生成算法 .74 基于三角形区域的线性插值 .94.1 加密点生成算法 .94.2 线性插值算法 .105 基于矩形区域的双三次样条插值函数构造 .125.1 双三次 B 样条函数插值函数求解方法 .125.2 二次加密点的生成 .136 基于似三棱柱体的三维地质体模型构建方法 .147 结论与展望 .17摘要：以盐城市水文地质钻孔实际检测含水层数据为基础，研究了盐城地区孔隙地下水承压含水层系统三维地质模型的构建的一般技术路线，并在 VC+环境下使用计算机可视化技术（VTK）设计实现了地质体建模系统原型的开发。首先设计钻孔数据模型，根据水文地质钻孔标准化卡片绘制水文地质钻孔模型,使用逐点插入法构建初始三角网，以此为基础，建立各含水层分层界面的初始 DEM；之后再根据三角形区域的线性插值算法实现对初始 DEM 的加密，为了满足二元插值函数的计算要求，本文设计了一种生成各分层界面计算节点规则矩形格网的算法；最后利用该加密后的计算节点矩形格网，使用二元双三次 B 样条函数进行二次加密，构建具有最优模拟效果的基于似三棱柱的三维含水层模型。关键字：水文地质钻孔；孔隙地下水；逐点插入法；线性插值；三维地质体建模；引言三维水文地质建模，往往是为了后续开展孔隙地下水流场数值模拟与可视化工作服务的。只有在一个准确、合理的水文地质体模型的基础上，才能使地下水的数值模拟误差控制在有效的范围以内；另外由于地质体模型的构建所使用到的初始采样数据一般比较稀疏，根据原始采样数据建立起的地层模型模拟效果较差。因此如何能够根据实际探测到的地层采样数据，通过不同的插值算法，构建起精度满足要求、模拟效果较好的最佳地质体模型，便成为了首先要考虑定的问题。一般说来，如果研究区域的采样点数据是以矩形点阵排列的，可以通过采用基于矩形区域代数插值函数、矩形区域的样条插值函数、矩形区域的最小二乘插值函数等对地质体表面进行逼近。但在实际问题中，采样点往往并不能够以矩形点阵的形式排列，此时，可以采用基于三角形区域的插值、康斯曲面以及移动曲面拟合法等插值算法对地质体表面进行加密。由于三维模型表面常常采用三角网来表示，因此本文选择基于三角形区域的线性插值函数作为地质体表面的逼近算法。然而基于三角形区域的线性插值实际上是以平面去逼近曲面。虽然在整个区域上可以保证“连续” ，但是在连接处却是不光滑的，即一阶导数不相等。因此需要使用样条函数对地层表面模型进行平滑。在三角形插值的基础上构建规则矩形格网，之后在此矩形格网上使用二元样条函数进一步内插，最终形成整个区域上光滑、连续的地质体表面模型。在加密后的地质体表面模型基础上，构建不同含水层界面的 TIN 模型。针对水文地质钻孔数据特点，考虑到含水层存在缺失的情况，设计可以识别不同含水层情况的三维地质体建模算法。水文地质体建模所采用的空间数据通常有似三棱柱、四面体、六面体等数据模型。由于水文地质钻孔被可近似表示为垂直的，且含水层分层界面 DEM 采用的是不规则三角形格网（TIN）结构；又由于地层中经常出现缺失的情况，直三棱柱往往会退化为四面体或者四棱锥。而直三棱柱、四面体或者是四棱锥，都可以统一采用似三棱柱的数据模型加以表示。似三棱柱可以精确模拟承压含水层对象的表面与内部结构，灵活性较强、生成算法较为简单。因此对于水文地质体的三维建模而言，似三棱柱是具有较高的实用价值和应用前景的一种空间数据模型。1 孔隙地下水含水层系统概念模型1.1 含水层的空间分布特点孔隙地下水含水层在空间上分布相对均匀，连续性较好。从上至下依次可划分为潜水含水层、第一承压含水层、第二承压含水层、第三承压含水层、第四承压含水层。不同的含水层之间一般由透水性很弱的隔水层加以间隔。地下水赋存于岩层的孔隙中，其分布范围与含水层的分布范围是基本一致的。1.2 承压含水层系统建模数据源由于孔隙地下水系统埋于地下，只有通过水文地质钻探与地球物理勘探的手段揭示含水层与隔水层的空间分布特征。其中水文地质钻探获取到的水文地质钻孔数据是描述孔隙地下水系统水文地质结构的基础的数据，也是孔隙地下水系统三维水文地质体模型构建的主要数据源。水文地质钻孔数据中，描述水文地质结构的数据有：钻孔地面高程、含水层的顶底板高程、隔水层的顶底板高程、含水层与隔水层的厚度、含水层与隔水划分数据、含水层岩性数据等。1.3 三维水文地质体模型特点水文地质钻孔数据的精度是由一定区域范围的水文地质勘查精度决定，在空间上呈离散分布，对于那些勘察盲区的水文地质结构只有通过水文地质专家根据第四系地层的沉积规律与自身的经验推断获取。另外，由水文地质钻探获到的反映水文地质结构的数据是一成不变的，仅反映勘探时刻的含水层与隔水层的空间分布情况，不能动态描述含水层与隔水层随着地下水位下降其空间结构也发生变化的特点。因此水文地质钻孔数据具有时间性的特点，构建的三维水文地质模型是一个静态模型，反映了某一水文地质钻探时间的隔水层与含水层的空间分布情况。2 水文地质钻孔模型构建方法水文地质钻孔模型可近似视作为具有一定直径和高度的分段圆柱体。不同的分段代表不同的含水层，不同的含水层种类，应以不同的颜色加以区别。由于钻孔模型是规则的几何体，因此可以采用规则格网数据模型来表示水文地质钻孔。规则格网数据模型，即结构化网格数据模型，内部是由具有规则拓扑结构的单元组成，而其外部的几何形状则是不规则的。它本质上可以看作是一系列等间距规则分布的点集，按照指定的结构（规则）分布形成的三维网格数据。2.1 钻孔数据结构设计设计的水文地质钻孔数据结构中应包含钻孔的位置信息、含水层信息，以及钻孔编号信息。水文地质钻孔数据结构Class CBoreInt id; /钻孔编号 double3 positionCoord; /钻井坐标float height; /钻孔深度int layerNumber; / 包含含水层层数vectorhydroCategory; /包含含水层的种类vectorhydroDepth; /包含的各含水层的层底深度;2.2 水文地质钻孔建模方法为了区分出不同含水层之间的隔水层，将潜水含水层和第一承压含水层之间的隔水层记为第一隔水层，依此往下分别记作第二隔水层、第三隔水层、第四隔水层和第五隔水层。这样一来，含水层从上往下正常的顺序应为：潜水含水层、第一隔水层、第一承压含水层、第二隔水层、第二承压含水层、第三隔水层、第三承压含水层、第四隔水层、第四承压含水层、第五隔水层。为了便于管理，自上往下依次编号为 0、1、2、3、4、5、6、7、8、9。含水层与编号的对应关系如表 2.1 所示。表 2.1 含水层编号含水层编号 含水层名称0 潜水含水层1 第一隔水层2 第一承压含水层3 第二隔水层4 第二承压含水层5 第三隔水层6 第三承压含水层7 第四隔水层8 第四承压含水层9 第五隔水层每一个水文地质钻孔是由一系列的分段圆柱拟合而成的，每一个分段圆柱对应于一种含水层介质。上一层的圆柱底面，对应于下一层的圆柱底面。生成时首先需要确定钻孔的分段数，也即包含的含水层层数 k，之后指定圆柱体的圆周拟合数 i，该拟合数决定了圆柱体圆周的圆滑程度。确定了这两个参数之后，依据如下算法，即可绘制出三维水文地质钻孔：（1）从数据库中读取钻孔位置信息和含水层信息，创建一个钻孔对象CBore aBore。（2）对于每一个钻孔对象，分段数 k = aBore.layerNumber；圆周拟合数 i需要手动指定，可以设置为 10；另外还需要指定圆周半径 rMax 和 rMin。此处应令 rMin = 0 deltaRad = (rMax - rMin)/(2- 1);（3）对于每一个圆柱分段 ki，获取层底圆周面中心点位置坐标， (x0,y0,z0 aBore. hydroDepth ki);（4）根据层底圆周面中心位置坐标，计算圆周拟合点坐标：for(int j=0 ; j RadiansFromDegrees(i *angel );x0 = radius * cos(theta) + x0;x1 = radius *sin(theta) + y0;Index = i + kIndex + jIndex;points-InsertPoint(Index ,x);（5）将每一层对应的的含水层编号赋给分段圆柱，最终形成 VTK 可识别的规则格网数据文件。导入到 VTK 可视化管线中显示。如图 2.1 所示。图中为了较为清楚地展示模型的属性与特征，在 z 方向上做了拉伸处理。图 2.1 三维水文地质钻孔建模3 分层界面三角网（TIN）3.1 三角网模型数据结构三角网由一系列的三角形单元组成，每一个三角形由可拆解为三个有向边。每个有向边又是有起点和终点两个端点确定的。基于此拓扑关系，设计数据结构如下：（1 ） 点数据结构Class CPointInt id /点索引号float x y z; /点的位置坐标int layerCategory; /所属含水层种类int zkBelong; /点所属钻孔的编号;（2 ） 有向边数据结构Class Cline Int id /有限变索引号long p0; /边起始点索引long p1; /边终点索引int useCount; /边被使用次数（3 ） 三角形单元数据结构Class CTriangleint id; /三角形索引号long L0; /第一边索引long L1; /第二边索引long L2; /第三边索引TIN 三角网，即将三角形数据 CTriangle 对象依次排列并编号形成的动态数组对象。3.2 三角形生成算法含水层分层界面可以采用 TIN 方法来进行建模。TIN 方法利用所有的采样点离散数据，按照优化组合的原则，把这些离散点连接成相互临街但不相交的三角形单元。对于 TIN 模型来说，Delaunay 三角网在地形拟合方面表现得最为出色。由 Delaunay 三角形构成额 TIN 三角网可写作 D-TIN。D-TIN 的生成算法很多，包括分治算法、三角网生长法、逐点插入法等。本文采用三角网生长法生成水文地质体含水层分层界面三角网。建模时各分层界面上的离散点，实际上为钻孔的原始采样点、与各地层的交点以及内插点。根据钻孔数据获取到每一层面所包含的离散点集合。采用三角网生长算法得地层界面三角网。本文三角形生长算法设计实现过程如下：（1）在同一层的离散点中任取一点 p0，选择离它最近的一点 p1 连接成有向边 L0,作为基准边。（2）遍历剩下的所有点，搜寻与基边两端点连线组成夹角最大的点。组成第一个三角形。将第一个三角形推入 TIN 三角网数组。三角形三边的可拓展属性均初始化为 1。该属性值为 1 表示该边只为一个三角形的边，如果为 2，则表示为两个三角形的公共边。一个边最多被两个三角形共享，因此，该属性值只可能为 1 或者 2。（3）对于三角网数组中的三角形 k，遍历其每一个边。对于每一个有向边，首先判断是否属于可扩展边，即其扩展属性是否为 1。如果可扩展，将与该边共三角形而不共边的那个点作为基准点。（4）遍历剩余所有离散点，如果与基准点在该有向边异侧，则求取其与该有向边两端点连线夹角。找到该夹角最大的离散点，形成新的三角形，并将该三角形推入数组。（5）每一个三角形三边依照（3） （4）两步扩展完毕后，三将网数组当前三角网索引 k 自增 1。重复步骤（3） （5） 。当 k 的值与三角网数组中包含三角形个数相同时，退出循环。此时所有的点都已经纳入到三角网中。三角网生成完毕。TIN 的拓扑结构是：三角网中包含的其实是每一个三角形的索引，每个三角形包含了三个有向边的索引，有向边包含了两个离散点的索引。离散点作为基础的几何数据，包含了每个点的位置和属性数据。最终生成的三角网效果如图 3.1 所示：图 3.1 分层界面的 Delaunay 三角网生成4 基于三角形区域的线性插值由于每一层的离散点个数较少且分布不均匀。无法满足孔隙地下水流模型数值模拟计算的要求，而且生成出来的地质体较为突兀，模拟效果较差。因此就需要对整个研究区域做插值处理。由于本文的侧重点在于地质体模型的建立以及美化，因此可以将加密点按照规则格网的形式加以排列，以便于平滑。为了方便后续要进行的地下水流模型数值计算过程中，在根据规则格网的点进一步做内插加密操作时，应将代表地下水位监测井和开采井的离散点纳入加密点集合中。由于已经建立起了初始三角网，且根据初始离散点的分布特征，所以采用基于三角形区域的线性插值方法。这一步骤，是为了得到规则的含水层分层界面 DEM 规则格网。4.1 加密点生成算法在对三角网做插值操作之前，首先需要生成定义规则格网的一系列格网点。格网点的生成原则应遵循如下两点：（1 ） 考虑到研究区域范围，点与点的间距定为 10 公里，即 10000 米。（2 ） 离散点应呈规则分布，取平面坐标位于三角形内部的离散点做内插。可以选取一个水位地质钻孔离散点作为起始点，采取延展法向四个方向扩充加密点。算法如下：（1） 以基准点（X0，Y0，Z0）为起点，X 坐标增加 10000 米，Y 坐标不变，创建新的点对象（X ， Y，0） ；（2） 遍历三角网中每一个初始三角形，寻找包含该点对象的三角形单元；（3） 搜寻结果不为空，则将该点推入格网点数组，线性插值出该点的 Z 值，更新点坐标为（X ，Y ，Z） ；（4） X 增加 10000 米，Y 不变，重复（2） （3）两步。（5） 如果搜寻结果为空，则 X = X0,Y 坐标增加 10000 米。创建新的点对象（X， Y，0） ；（6） 同样方式，遍历初始三角形数组，寻找包含该点对象的三角形单元；（7） 搜寻结果不为空。则将点（X， Y， 0）推入各网点数组，线性插值获取z 值，更新点坐标为（X , Y ,Z）;（8） Y 增加 10000 米，X 不变。 。重复（5） （6） ，直到搜寻结果为空，表明基准点东北方所有各网点生成完毕，令 X = X0 10000；Y = Y0，开始生成第二个方向上的所有离散各网点（9） 依此生成基准点东南、西北、西南方向以及东南西北方向线上的所有离散格网点。并经过线性插值操作，最终计算出加密后的含水层分层界面DEM。如图 4.1 所示图 4.1 规则格网离散点的生成4.2 线性插值算法基于三角形区域的线性插值算法如下:设有三角形区域 如图 4.2 所示。其顶点记为 Ai =(xi , yi );三个顶点对应的值记为 u( xi , yi );图 4.2 三角形单元示意图则线性插值函数可写为：其中 U(x,y)在空间表示一个平面方程。使用平面来逼近曲面。该方法又称为“板法” 。由于在整个区域上，相邻三角形公共边对应有两个公共节点，两点决定一条直线，而这个直线即为相邻平面的交线。因此该方法可以保证相邻单元在连接处连续。但是无法保证平滑。因此如果想要得到模拟效果较为光滑的水文地质体模型，就需要在线性内插得到所有规则格网点的基础上，利用样条函数，进行第二次内插。最终插值效果如图 4.2 所示：图 4.2 初始加密格网点的生成5 基于矩形区域的双三次样条插值函数构造5.1 双三次 B 样条函数插值函数求解方法插值的方法主要可以分为全局插值、分块插值以及移动曲面拟合法两类。由于格网 DEM 差值单元为规则矩形区域，所以可以使用基于矩形区域的代数插值或者分片插值来获得全局函数。但是由于此处格网数目较多，逼近多项式会出现高次项而引起龙格（Rung）现象，导致函数值出现震荡；而利用分片矩形的小区域代数插值逼近，即分片插值又使得整体区域虽然可以消除龙格现象，且有较好的连续性，但光滑性较差。一种较好的解决办法就是使用基于矩形区域的样条函数。由于三次样条函数具有次数低、二阶导数连续，连接处较为光滑等优点，本文通过构造二元双三次 B 样条函数作为规则格网矩形区域上的插值函数。基于 DEM 的样条函数插值问题，可等价为如下问题：即对于给定的 N + 1个规则分布的离散格网点， （x i,yj）,in，固定 yt,则有取 t = 0,1,2m，对于每一个 yt,都可写出一个包含 n+1 个方程的方程组,最终求解得到 ij第二步：记第三步：求解：第四步：建立方程组解出 Cij 的值，从而确定差值曲面 U(x,y)，求解完毕。根据该插值函数。计算各二次加密点的值。5.2 二次加密点的生成对于规则格网点额进一步与加密，即“二次加密” ，应尽可能均匀全面地分布在整个区域中。除了要求要有足够的密度以保证三维地质模型的光滑性外，还应考虑到后续工作中地下水流场模型的验证与校核工作，因此二次加密点应遵循如下原则生成：（1 ） 首先应包含地下水水位监测井采样数据、地下水开采井采样数据。（2 ） 应对格网进行适当加密，在每个格网内部，延 X 和 Y 方向均匀地划分二级加密点。最终形成经过二次加密后的离散采样点集合。对每一层含水层分界面重复上述插值操作，记录每个采样点所属的分层界面，最终得到所有分层界面离散采样点的集合，通过调用三角网生成算法，重新生成各含水层分界面 TIN 三角网。该三角网即为平滑连续的分层界面模型。6 基于似三棱柱体的三维地质体模型构建方法似三棱柱体建模是以已经生成好的各层 TIN 三角网为基础，在相邻的分层面之间构建一系列的似三棱柱单元。由于本文将钻孔数据近似成为了竖直的圆柱体，因此似三棱柱的形态仅包括如图 6.1 所示的三种类型。、图 6.1 基于竖直钻孔的似三棱柱的三种类型对于上图中的几种类型三棱柱，为了统一数据结构，均视作由六个顶点p0， p1,p2,p3,p4,p5 组成的数据结构。其中 p0,p1,p2 构成下三角形，p3,p4,p5 构成上三角形。三角形顶点排列方向应一致。本文均已顺时针排列。如图 6.2 所示。图 6.2 似三棱柱单元组成结构为了便于绘制，首先设计似三棱柱单元数据结构：Class CPrism Int id /三棱柱索引号int p0,p1,p2; /上三角形顶点索引int p3,p4,p5; /下三角形顶点索引int category; /三棱柱所属含水层编码每一个含水层，即是由一系列这样的三棱柱对象组成。一个含水层对象，可以用一个以三棱柱对象为元素的动态数组对象来表示。似三棱柱的生成过程，是把每个三角形的三个顶点所属分层界面种类进行比较，并按照一定的规则向下扩展而形成的。地层编码已经设计成为了自上往下顺序排列。建模的主要过程如下：（1） 由于不同层加密离散点在平面位置上是相同的，因此可以根据每一层具有相同平面位置属性的加密点序列，生成一个具有钻孔数据结构 CBore的虚拟钻孔。经过该步骤操作之后，整个区域密集均匀分布着很多虚拟钻孔。（2） 生成第 li(i=0)层含水层：获取第一层地层分界面三角网，对于每一个小三角形单元，向下寻找地层层底深度。寻找算法如下a） 设置当前搜寻地层编码 int p，如果 p aBore.hydroCategory.size()。表示搜寻的地层编码已经超出该虚拟钻孔所包含的的地层数，那么此时只需将地层层底深度设置为顶层深度 减去 aBore.hyDrodepth 数组中的最后一个元素即可;b） 如果 p aBore.hyDroCategoryp , 那么将钻孔地层深度数组中最后一个层底深度作为该含水层层底深度。d） 如果 li aBore.hyDroCategoryp;层底深度 z = 层顶深度 - aBore.hyDroDepthp - 1；e） 如果 li = aBore.hyDroCategoryp