三维地质建模的数学模型与显示方法.pdf

工程地质计算机应用 20062006 年 第 3 3 期 总 4343 期 1 1 专题论述 问题探讨 三维地质建模的数学模型与显示方法 曾钱帮 1 何小萍2 1 中国科学院地质与地球物理研究所 北京 100029 2 北京软通动力科技有限公司 北京 100027 摘要 基于离散数据集的曲面插值拟合方法 精确通过工程勘察数据点 获得光滑连续的地质 界面的数学模型 可以用于表达地形 地下水位面 岩层面 构造面等各种地质界面和岩土体物 理力学参数的空间分布 单值界面的数学模型中的插值型滑动最小二乘法是局部插值方法 可避 免全局插值方法的缺憾 编制程序生成 AutoCAD 脚本文件的地质界面计算机显示方法的优点是对 编程技术要求不高 简单实用 可充分利用 AutoCAD 软件强大的图形显示功能 关键词 工程地质 三维建模 勘测数据 单值曲面 拟合函数 地层曲面 计算机显示 1 前言 工程地质三维建模与可视化是应用计算机图形学和图像处理技术 是将工程地质勘测数 据和工程地质岩土体力学数值模拟分析的计算结果转换为图形图像在计算机屏幕上显示出 来并进行交互处理的理论 方法和技术 复杂地质体中的各类地质信息都可以被看作是三维 空间的函数 利用各种野外勘测数据分别建立相应的曲面拟合函数 进而利用计算机建立三 维地质模型 逼真反映地质结构全貌 达到直观地表达地质信息的分布规律 提高对于地质 规律的认识 指导地质工程项目的勘测施工的目的 因此 工程地质三维建模与可视化研究 有其重要的理论和现实意义 工程地质三维建模与可视化研究中 地质界面的数学模拟是基础 由于地质界面必须精 确通过控制点 工程勘测数据点 通过离散的工程勘测数据点建立三维地质模型及其计算 机图形可视化显示属于重构问题 与计算机图形学中对于机械设计特别有效的形体的构造问 题是有根本区别的 所以针对机械设计发展起来的自由曲线曲面造型技术 1 在地质层面模拟 中的应用受到极大限制 另外 现今热门的通过 Voronoi 图和 Dulauny 三角剖分 2 在空间构 造不规则三角网 Triangular Irregular Net TIN 方法构造层面 一方面由于原始数据点一般 相距较远 常需要进行三角网的插值加密 另一方面 插值曲面不光滑 无法求出层面上某 点沿坐标轴的坡向 坡度 曲率和产状 TIN 方法无数学曲面的解析表达式 就无法求得对 于 x 或 y 的一阶和二阶偏导 而地质层面的走向 倾向和倾角又与曲面方程的偏导数有一定 的关系 3 用于地质界面的模拟也是不合适的 离散数据集的单值曲面拟合插值法是较为 成熟 使用频率高的方法 与距离成反比的加权方法 径向基函数插值方法 样条曲面方法 等 拟合过程大都是全局性的 当增加 删除和修改数据点时都要重新计算 有些方法还必 须要求解维数很大的线性方程组 这就使得这些方法的效率和稳定性大大降低 寻求局部插 值拟合方法就成为今后地质界面数学模拟的发展目标 对于用来模拟复杂褶皱的多值曲面拟 2 2 工程地质计算机应用 20062006 年 第 3 3 期 总 4343 期 合插值函数也是一个公认的难题 地质体通常是不规则形体 由多个各种成因类型形状各异的结构面 4 围限而成 有一定 的物质组成 结构和赋存于一定的地质环境中 遭受过多期次的变形和破坏 5 其工程地质 条件复杂 千变万化 对于单个单值层面和地质条件简单的多层连续地层的模拟和计算机图 形显示比较简单 但是 当把多个这样的单值层面在空间叠加 考虑出现断层错断岩层 地 层不整合和结构面的组合 相互穿插等地质现象 组成地质体的三维模型 并能进行切剖 其复杂度和建模技术难度会爆炸性地增加 一些理论问题 如三维拓扑结构分析 曲面求交 运算 形体相互遮挡的消隐等 需亟待解决 2 工程地质层面的数学模型 离散数据拟合插值所构造的层 曲 面模型是对地质信息在复杂地质体中分布的数学抽 象描述 为绘制和显示地质信息的空间分布提供了重要的方法基础 地质信息的插值和拟合 函数要根据实际勘测数据建立 实测数据越丰富精确 得到的地质模型越能够真实描绘出这 些信息的空间分布规律 另外 由于地质信息数据的特殊性 在进行空间数据的插值时 必 须考虑许多约束条件及相关的地质学原理 对于不同特点的地质信息 需采用不同的拟合函 数 才能形成准确可靠的模型 2 1 单值层面的拟合方法 空间离散数据的插值和拟合是构建三维模型的基础 地表地形测量数据 X 坐标 Y 坐标 和地表高程 Z 地下水位埋深测量信息 地下水位测点地表 X 坐标 Y 坐标和水位埋深 h 等 的单值曲面图形生成可归结为双自变量离散数据的插值和拟合 即 假设二维平面上有一组 n 个点 xi yi i 1 2 n 并有 fi f xi yi 插值问题就是要构造一个具有 C1连续的函数 F x y 使其在 xk yk k 1 2 n 点的函数值为 fk 并可根据该函数推求出区域范围内其他 任意点的函数值 从而重构一个具有连续特征的量值在三维空间的变化情况 空间单值曲面 插值函数有以下构造方法 如与距离成反比的加权方法 Shepard 方法 径向基函数插值法 Multiquadric 方法 基于平面弹性理论的曲面样条插值法 曲面样条函数法 插值型滑动最 小二乘拟合方法等 它们同样适用于单个连续地层界面 地球物理勘探数据 地球化学勘探 数据以及岩土体物理力学参数在地质体空间的分布 2 1 1 与距离成反比的加权方法 Shepard s Method 此方法又称为最小二乘距离加权插值算法 其基本思路是将插值函数 F x y 定义为各实 测数据点 fk的加权平均 即点 xk yk 的值 fk对于 F x y 的影响与 xk yk 至 x y 的距离成反 比 6 令 2 k 2 kk y y x xd 表示由 x y 点到 xk yk 点的距离 则插值函数 F x y 可表示为 n k kk n k k n k k k yxWf yxd yxd f yxF 1 1 2 1 2 1 工程地质计算机应用 20062006 年 第 3 3 期 总 4343 期 3 3 式中权函数 n k n kj j n kj j k yxd yxd yxW 1 2 2 有如下性质 y x Wk非负 y x Wk是C0连续 的 y x Wk kj 当k j时 kj 1 否则 kj 0 y x Wk 1 具有加权性质 Shepard方法的插值结果只能是 0 C连续的 当增加 删除或改变一个点时 权函数 y x Wk均需重新计算 因而该方法是一个全局插值算法 图1是根据离散数据利用与距离 成反比的加权方法绘制的地层曲面 2 1 2 径向基函数 Radial Basis Function 插值法 Multiquadric 方法 该方法采用的插值函数为 6 njcyyxxayxF n j jjj 2 1 1 2 1 222 其中c为常数 一般取1 将 n 个点 y x ii 的实测值 i f代入上式建立联立方程 nicyyxxaf jiji n j ji 2 1 2 1 222 1 求解此线性方程组可获得待定系数 aj j 1 2 n 将这些值代回插值函数式 即为通过 各实测数据点且处处连续光滑的曲面方程 在数据点数量不太大的情况下 Multiquadric 法计算不太麻烦 这一方法提出的近 20 年 间 在水文测量 大地测量 地质及采矿 地球物理等领域得到广泛应用 效果良好 2 1 3 基于平面弹性理论的曲面样条插值法 对地表地形或地下水位等离散数据进行函数拟合可借助平面弹性理论 把计算域视为无 限延伸平面 在平面 y x 11 点上给定一垂直位移 1 A时 平面各处将均会随着产生位移 若 限制距离 y x 11 点半径为 R 以外位移为零时 则由 1 A引起距离 y x 11 点为 r r 8 式中 ri x xi 为xi与x的距离 rmi为结点i的影响半径 为一正的小值 k为正整数 由式 8 容易得出 wi x 在ri 0 内存在关于坐标的k 1阶连续偏导数 则由式 3 得出的形函数ni x 也存在关于坐标的k 1阶连续偏导数 值越小 wi在i节点自身取值越大 而在远离i点处取值越小 并在距离rmi以外为0 当 0时 权函数具有奇异性 这种情况下 滑动最小二乘拟合符合插值条件 奇异权虽然在理论上可取 但在数值计算时造成滑动最小二乘法正规函数的病态 所以 实际计算中应取 0 此时 近似函数并不精确通过每一个已知点 但若 足够小 权函数 接近于奇异 滑动最小二乘法接近于插值型滑动最小二乘法 及k的选取具有一定程度的任 意性 但选取得好可提高计算精度 一般选用k 4 1 0 计算结果较好 rmi需适当选取 在尽量减少计算量的同时满足式 3 中A的非奇异性 在结点均匀分布时可取为 mi n r c 9 式中 n 3 线性基 6 二次基 10 三次基 c为结点分布密度 为大于1的 系数 一般取 4 6 2 2 多层连续地层曲面的拟合函数 地层面在地质体中按一定的层序排列 每一地层有相应的厚度 上面讲述的单值曲面拟 合函数也可用于描述连续多层地层界面 但要对其进行相应的修正 计算域中各地层厚度可 能差别很大 在对地层面进行函数拟合时 为使拟合函数更好地逼近原始地层面 需对各层 面由上到下作层序编号 且对不同地层面指定不同的函数值常量来描述 假定地质体中有L 个地层面 则相应有 L 1 个地层 同一地层各处厚度也可能不同 各个地层的参照厚度 H1 H2 HL 1可取各勘测点相应地层厚度的平均值 利用各地层面参照厚度计算出各地层面的 函数值常量 令最上部第一层面定义函数值V1 1 则第i层面的函数值常量Vi可表示为 7 12 15 1 1 max i j i ji H H VV i 2 3 L 式中Hmax为H1 H2 HL 1中的最大值 为了拟合各地层面 需要一组N点的实测层面定位数据 xi yi zi 和k 其中 xi yi zi 为地 层面空间测点坐标 k为该测点i所在地层面的序号 则 xi yi zi k与Vk满足关系式 Ni R r R r R r AV ijij N i ij jk 2 1 1ln 2 2 2 2 1 2 2 式中rij2 xi xj 2 yi yj 2 zi zj 2 对N个实测数据点解此线性方程组可得待定系数Aj j 1 2 n 将这些值代回插值函数式 即为通过各实测数据点且处处连续光滑的地层界面 工程地质计算机应用 20062006 年 第 3 3 期 总 4343 期 7 7 这种多层连续地层的插值拟合方法的优点是 当在受条件限制 数据分布极不均匀 某 些层位的数据较少的情况下 可以通过相邻层位的形态 数据的相似性及相关关系 或者利 用相邻层位在空间的延展趋势的协调一致性 来修正或改善这些数据稀疏的地层面的模型 3 地质界面的计算机显示方法 科学计算可视化技术是在计算机图形学基础上发展起来的一个崭新的领域 用形象 具 体的图形图像表达大量复杂 抽象和多维的测量数据 工程计算数据和科学计算数据所蕴涵 的内容 它将图形生成技术 图形处理技术和人机交互技术结合在一起 从而可以更快地分 析计算结果 更有效地控制计算过程 在地质勘探和工程地质力学分析中经常遇到大量数据 的处理和解释问题 尤其是地球物理数据采集和处理 工程地质勘测数据的处理 空间分布 的地球物理场和属性参数 包括岩土体物理力学参数等 地质目标层位 如地层 断层面 剖面 矿藏 的位置 延展和交切情况以及工程地质力学的数值模拟分析计算和处理结果的 解释等 利用科学计算可视化技术不仅丰富了体现数据的手段 而且为工程地质力学分析的 方法解释提供可靠依据 3 1 基于 OpenGL 的三维地层曲面模型可视化图形显示方法 OpenGL是美国高级图形和高性能计算机系统公司SGI所开发的一套高性能三维图形处 理系统 已被设计成适合于各种计算机操作系统下的三维图形应用程序编程接口 Application Programming Interface API 目前它已成为开放的国际图形标准 运用OpenGL 通过一系列基本的几何图元 Geometry Primitives 16 点 直线 三角片 带 四边形片 带 来建立地质界面模型 在绘制地质曲面时 通过对已知离散的勘测数据的 拟合插值形成空间规则网格 再划分成一系列三角片进行渲染 17 对于相邻地层之间的空间 由于相邻地层界面上对应的网格点组成的四边形片位于同一平面上 可以直接画出四边形带 以缝合上下相邻层面 从而形成侧面 18 根据地层岩性进行不同的颜色填充 形成地质模型 OpenGL提供了大量的图形变换函数 在编制程序时无需进行复杂的矩阵运算 就可方 便地将三维地质体模型显示在屏幕窗口并进行平移 缩放 旋转等操作 为了增强图形的真 实感 OpenGL还提供了隐藏线面消除 着色和光照 纹理映射和反走样等技术 简化了编 程 另外 OpenGL还提供了用双缓存技术和显示列表技术 可充分利用硬件加速功能 并 且利用MFC 微软基础类库 Microsoft Foundation Class 编写对鼠标和键盘的消息响应函数 19 实现对三维地质模型的动态显示和动画效果 3 2 编制程序生成 AutoCAD 脚本文件的图形显示方法 AutoCAD脚本文件 AutoCAD Script File 类似于DOS操作系统中的批处理文件 它可以将 不同的AutoCAD命令组合起来 并按确定的顺序自动连续地执行 22 脚本文件是文本文件 扩展名为 SCR 用户可使用任一文本编辑器来创建脚本文件 因为脚本文件可使一些命令 序列自动执行 所以常用来产生 编辑或观看图形 如幻灯放映 初始的图形设置等 三维地质模型的计算机显示 可以采用任何一种高级语言 如VC VB Delphi Java等 设计用户交互界面 对绘图所需要参数进行计算生成 然后确定AutoCAD命令 命令选项 8 8 工程地质计算机应用 20062006 年 第 3 3 期 总 4343 期 命令序列等 最后生成扩展名为SCR的 AutoCAD脚本文件 在AutoCAD中用SCR 命令来执行脚本文件 完成计算机显示 该 法的优点是对编程技术要求不高 简单实 用 可充分利用AutoCAD软件强大的图形显 示功能 表1列出了一组n 20的地表测点数据 采用2 1 3节的基于平面弹性理论的曲面样 条插值法和2 2节多层连续地层曲面的拟合 函数 通过编制程序生成AutoCAD脚本文件 得到的这组数据的地表曲面和多个地层曲 面的计算机图形 图2 由图2可以看出 曲面为连续光滑面 且函数通过各实测点 表 1 地表离散点实测数据 i xi yi zi i xi yi zi i xi yi zi i xi yi zi 1 0 0 406 180 50 35111001003016150 250 60 2 50 0 657 200 50 50121501006017200 250 50 3 80 0 408