|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
三维地质建模的数学模型与显示方法 曾钱帮1
何小萍2 （1
中国科学院地质与地球物理研究所 北京 100029） (2
北京软通动力科技有限公司 北京 100027) 【摘要】基于离散数据集的曲面插值拟合方法，精确通过工程勘察数据点，获得光滑连续的地质界面的数学模型，可以用于表达地形、地下水位面、岩层面、构造面等各种地质界面和岩土体物理力学参数的空间分布。单值界面的数学模型中的插值型滑动最小二乘法是局部插值方法，可避免全局插值方法的缺憾。编制程序生成AutoCAD脚本文件的地质界面计算机显示方法的优点是对编程技术要求不高，简单实用，可充分利用AutoCAD软件强大的图形显示功能。 【关键词】工程地质
三维建模 勘测数据 单值曲面 拟合函数 地层曲面
计算机显示 1
前言 工程地质三维建模与可视化是应用计算机图形学和图像处理技术，是将工程地质勘测数据和工程地质岩土体力学数值模拟分析的计算结果转换为图形图像在计算机屏幕上显示出来并进行交互处理的理论、方法和技术。复杂地质体中的各类地质信息都可以被看作是三维空间的函数，利用各种野外勘测数据分别建立相应的曲面拟合函数，进而利用计算机建立三维地质模型，逼真反映地质结构全貌，达到直观地表达地质信息的分布规律、提高对于地质规律的认识、指导地质工程项目的勘测施工的目的。因此，工程地质三维建模与可视化研究有其重要的理论和现实意义。 工程地质三维建模与可视化研究中，地质界面的数学模拟是基础。由于地质界面必须精确通过控制点（工程勘测数据点），通过离散的工程勘测数据点建立三维地质模型及其计算机图形可视化显示属于重构问题，与计算机图形学中对于机械设计特别有效的形体的构造问题是有根本区别的，所以针对机械设计发展起来的自由曲线曲面造型技术[1]在地质层面模拟中的应用受到极大限制。另外，现今热门的通过Voronoi图和Dulauny三角剖分[2]在空间构造不规则三角网（Triangular
Irregular Net, TIN）方法构造层面，一方面由于原始数据点一般相距较远，常需要进行三角网的插值加密，另一方面，插值曲面不光滑，无法求出层面上某点沿坐标轴的坡向、坡度、曲率和产状（TIN方法无数学曲面的解析表达式，就无法求得对于x或y的一阶和二阶偏导，而地质层面的走向、倾向和倾角又与曲面方程的偏导数有一定的关系[3]），用于地质界面的模拟也是不合适的。离散数据集的单值曲面拟合插值法是较为成熟、使用频率高的方法。与距离成反比的加权方法、径向基函数插值方法、样条曲面方法等，拟合过程大都是全局性的，当增加、删除和修改数据点时都要重新计算，有些方法还必须要求解维数很大的线性方程组，这就使得这些方法的效率和稳定性大大降低，寻求局部插值拟合方法就成为今后地质界面数学模拟的发展目标。对于用来模拟复杂褶皱的多值曲面拟合插值函数也是一个公认的难题。 地质体通常是不规则形体，由多个各种成因类型形状各异的结构面［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]。令：
式中权函数
有如下性质：
①Wk(x,y)非负；
②
Wk(x,y)是C0连续的；
③
Wk(x,y)
=
Shepard方法的插值结果只能是
C
0连续的，当增加、删除或改变一个点时，权函数
Wk(x,y)
均需重新计算，因而该方法是一个全局插值算法。图1是根据离散数据利用与距离成反比的加权方法绘制的地层曲面。 2.1.2
径向基函数（Radial Basis Function）插值法（Multiquadric方法） 该方法采用的插值函数为[6]：
其中c为常数，一般取1。将
n 个点
(xi, yi)的实测值
fi代入上式建立联立方程：
求解此线性方程组可获得待定系数aj
(j=1,2,…,n)，将这些值代回插值函数式，即为通过各实测数据点且处处连续光滑的曲面方程。 在数据点数量不太大的情况下，Multiquadric法计算不太麻烦。这一方法提出的近20年间，在水文测量、大地测量、地质及采矿、地球物理等领域得到广泛应用，效果良好。 2.1.3
基于平面弹性理论的曲面样条插值法 对地表地形或地下水位等离散数据进行函数拟合可借助平面弹性理论，把计算域视为无限延伸平面，在平面(x1,
y2)
点上给定一垂直位移
(A1)
时，平面各处将均会随着产生位移，若限制距离
(x1,
y2)
点半径为R以外位移为零时，则由
(A1)
引起距离(x1,
y
2)
点为r(r<R) 处的点(x,y)的位移
ZA(r)为：
式中r12
= (x - x1)2 + ( y - y1)2，R为平面影响半径，其影响函数zA(r)随r1的增大而逐渐减小；它在r
=0和r =R时dzA(r)/dr = 0。利用上式对N个实测数据点进行叠加可以得到离散点数据的曲面拟合函数[7,
20]：
为使
(x1,
y1)
点处F(x,y)等于实测数据点值
(F1)
，需建立联立方程组：
式中rij2=(xi－xj)2+(yi－yj)2。对N个实测数据点求解此线性方程组可得待定系数
(Aj)
(j=1,2, … ,n)，将这些值代回插值函数式，即为通过各实测数据点且处处连续光滑。 2.1.4
曲面样条插值法 该方法的优点是原始数据点可以任意排列，拟合出的光滑曲面通过已知数据点。其数学表达式[3,21]为：
式中，
a0, a1, a2, fi(i=1,2,...,n)为待定系数；ri2
= (x-xi)2+(y-yi)2；
改写成矩阵形式:
AX=B
其中，
X=(F1,F2,...Fn, a0, a1, a2)T B=(z1,z2,...zn,
0,
0,
0)T。 以上方法都是全局插值算法，当增加、修改和删除数据点时均需重新计算权函数或线性方程组的解，算法的稳定性和效率都明显下降。 2.1.5
插值型滑动最小二乘法 A、基本公式 设场函数为u(x)，其中x
= (x,y)T为场点，给定u(x)在n个节点上的值：
u(xi)=ui, i=1, 2, ..., n
则使用滑动最小二乘法（Moving
Least Squares Method）可确定u(x)的一个近似函数：
其中，a(x)为m维系数向量，p(x)为m维基向量。 为了保证收敛性，p(x)应取自完全多项式，如：一维情况下，
由滑动最小二乘法可将式（1）写成如下形式[8,
9]：
式中，wi(x)为节点的权函数在点x
= (x,y)T的取值，ni(x)为i节点的形函数在点x的取值。 由式（3）可得到形函数关于坐标的偏导数为：
Ak-1= -Ak-1AkA-1
(7) B、权函数 式（2）所建立的近似函数Gu(x)并不能保证通过每个已知点，即不一定满足插值条件Gu(xi)=ui
(i=1,2,…,n)，这与权函数的选取有关。如果将权函数wi(x)取一类特殊形式，则可使Gu(x)满足插值条件，这就是插值型滑动最小二乘法。有如下结论[10、11]：设wi(x)在已知数据点xi
奇异，即：
则
Gu(x)满足插值条件
Gu(xi)=ui(i=1, 2,
..., n)。 权函数的选取非常重要，它直接关系到曲面（地质层面）的计算精度，权函数的选取应遵循如下几个原则：①非负；②某节点的权函数应在自身取最大值，由近及远逐渐衰减，且在某个影响半径之外为0或可忽略不计，可使近似计算局部化，使一任意点的最小二乘计算只限于在该处权不为零的少数点，以减少计算量；③连续可导，以保证近似函数光滑。
式中，ri=║x-xi║为xi与x的距离，rmi为结点i的影响半径，e
为一正的小值，k为正整数。 由式（8）容易得出：wi(x)在riÎ[0,+¥
)内存在关于坐标的k-1阶连续偏导数，则由式（3）得出的形函数ni(x)也存在关于坐标的k-1阶连续偏导数；e值越小，wi在i节点自身取值越大，而在远离i点处取值越小，并在距离rmi以外为0。当e
=0时，权函数具有奇异性，这种情况下，滑动最小二乘拟合符合插值条件。 奇异权虽然在理论上可取，但在数值计算时造成滑动最小二乘法正规函数的病态，所以实际计算中应取e
>0，此时，近似函数并不精确通过每一个已知点，但若e
足够小，权函数接近于奇异，滑动最小二乘法接近于插值型滑动最小二乘法。e
及k的选取具有一定程度的任意性，但选取得好可提高计算精度。一般选用k=4，e
=1.0，计算结果较好。rmi需适当选取，在尽量减少计算量的同时满足式（3）中A的非奇异性，在结点均匀分布时可取为：
式中，n
= 3（线性基），6（二次基），10（三次基）；c为结点分布密度；a为大于1的系数，一般取a=4~6。 2.2
多层连续地层曲面的拟合函数 地层面在地质体中是按一定的层序排列，每一地层有相应的厚度。上面讲述的单值曲面拟合函数也可用于描述连续多层地层界面，但要对其进行相应的修正。计算域中各地层厚度可能差别很大，在对地层面进行函数拟合时，为使拟合函数更好地逼近原始地层面，需对各层面由上到下作层序编号，且对不同地层面指定不同的函数值常量来描述。假定地质体中有L个地层面，则相应有（L－1）个地层，同一地层各处厚度也可能不同，各个地层的参照厚度H1,H2,...HL-1可取各勘测点相应地层厚度的平均值，利用各地层面参照厚度计算出各地层面的函数值常量。令最上部第一层面定义函数值V1=1，则第i层面的函数值常量Vi可表示为[7，12－15]：
式中Hmax为H1,H2,...HL-1中的最大值。 为了拟合各地层面，需要一组N点的实测层面定位数据(xi,
yi, zi)和k，其中(xi, yi, zi)为地层面空间测点坐标，k为该测点i所在地层面的序号，则(xi,
yi, zi)、k与Vk满足关系式：
式中rij2
= (xi－xj)2 + (yi－yj)2+(zi－zj)2。对N个实测数据点解此线性方程组可得待定系数Aj
(j = 1, 2, … , n)，将这些值代回插值函数式，即为通过各实测数据点且处处连续光滑的地层界面。 这种多层连续地层的插值拟合方法的优点是，当在受条件限制、数据分布极不均匀、某些层位的数据较少的情况下，可以通过相邻层位的形态、数据的相似性及相关关系，或者利用相邻层位在空间的延展趋势的协调一致性，来修正或改善这些数据稀疏的地层面的模型。 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命令、命令选项、命令序列等，最后生成扩展名为SCR的AutoCAD脚本文件。在AutoCAD中用SCR命令来执行脚本文件，完成计算机显示。该法的优点是对编程技术要求不高，简单实用，可充分利用AutoCAD软件强大的图形显示功能。 表1列出了一组n=20的地表测点数据。采用2.1.3节的基于平面弹性理论的曲面样条插值法和2.2节多层连续地层曲面的拟合函数，通过编制程序生成AutoCAD脚本文件得到的这组数据的地表曲面和多个地层曲面的计算机图形（图2）。由图2可以看出，该曲面为连续光滑面，且函数通过各实测点值。 表1
地表离散点实测数据
|
表示由(x,y) 点到(xk, yk)点的距离，则插值函数F(x
, y) 可表示为：
,当k=j时，
为调节曲率大小的经验参数，一般取10-2~10-6；n为原始数据点个数。
n+3个待定系数可通过下列方程组求得： 
（1）
；对于二维情况，p(x)可如下选取： 
（2）
（3）
（4）
（5）
(6)
（8）
（9）
(i=2, 3, … ,L)