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

2
R2
ln
rij
2
R2
+1−
rij
2
R2
)
式 中 rij2=(xi － xj)2+(yi － yj)2 。 对 N 个 实 测 数 据 点 求 解 此 线 性 方 程 组 可 得 待 定 系 数
A j (j=1,2, … ,n)，将这些值代回插值函数式，即为通过各实测数据点且处处连续光滑。
2.1.4 曲面样条插值法 该方法的优点是原始数据点可以任意排列，拟合出的光滑曲面通过已知数据点。其数学 表达式[3,21]为：
(x - x k ) 2 + (y - y k ) 2 表示由(x,y) 点到(xk, yk)点的距离， 则插值函数 F(x , y)
F ( x, y ) =
∑ [d
k =1 n
n
k
∑ [d
k =1
n ( x, y ) ] = ∑ f k ⋅ Wk ( x , y ) 1 k =1 2
fk
k
( x, y ) ]
改写成矩阵形式: AX = B 其中，
⎡ 0 r12 2 ln(r12 2 + ε ) ⎢ 2 2 0 ⎢ r12 ln(r12 + ε ) ⎢ # # ⎢ 2 2 2 2 r ln( r + ε ) r ln( r 1, n −1 2, n −1 2, n −1 + ε ) A = ⎢ 1,n −12 ⎢ r ln(r 2 + ε ) r2, n 2 ln(r2, n 2 + ε ) 1, n ⎢ 1, n 1 1 ⎢ ⎢ x1 x2 ⎢ ⎢ y1 y2 ⎣
工程地质计算机应用
※
2006 年 第 3 期 ※
总 43 期
1
专题论述
问题探讨
三维地质建模的数学模型与显示方法
曾钱帮 1 何小萍 2
（1 中国科学院地质与地球物理研究所 北京 100029） (2 北京软通动力科技有限公司 北京 100027)
【摘要】基于离散数据集的曲面插值拟合方法，精确通过工程勘察数据点，获得光滑连续的地质 界面的数学模型，可以用于表达地形、地下水位面、岩层面、构造面等各种地质界面和岩土体物 理力学参数的空间分布。单值界面的数学模型中的插值型滑动最小二乘法是局部插值方法，可避 免全局插值方法的缺憾。编制程序生成 AutoCAD 脚本文件的地质界面计算机显示方法的优点是对 编程技术要求不高，简单实用，可充分利用 AutoCAD 软件强大的图形显示功能。 【关键词】工程地质 三维建模 勘测数据 单值曲面 拟合函数 地层曲面 计算机显示
z A (r ) = A1 (
r12 r12 r12 ln + 1 − ) R2 R2 R2
式中 r12 = (x - x1)2 + ( y - y1)2，R 为平面影响半径，其影响函数 zA(r)随 r1 的增大而逐渐减 小；它在 r =0 和 r =R 时 dzA(r)/dr = 0。利用上式对 N 个实测数据点进行叠加可以得到离散点 数据的曲面拟合函数[7, 20]：
0
图 1 加权法绘制的地层曲面
Shepard 方法的插值结果只能是 C 连续的，当增加、删除或改变一个点时，权函数
Wk (x, y) 均需重新计算，因而该方法是一个全局插值算法。图1是根据离散数据利用与距离
成反比的加权方法绘制的地层曲面。 2.1.2 径向基函数（Radial Basis Function）插值法（Multiquadric 方法） 该方法采用的插值函数为[6]：
2
工程地质计算机应用
式中权函数
2006 年 第 3 期
总 43 期
3
Wk ( x, y ) =
∏ [d
n j≠k n n k =1 j ≠ k
j
( x, y )
j
]
2
∑∏ [d
( x, y )
]
2
有如下性质：① Wk (x, y) 非负；② Wk (x, y) 是C0连续 的；③ Wk (x, y) = δ kj ,当k=j时， δ kj =1，否则 δ kj =0； ④Σ Wk (x, y) =1，具有加权性质。
T
0 1 xn yn
1 0 0 0
xn 0 0 0
y1 ⎤ ⎥ y2 ⎥ # ⎥ ⎥， yn −1 ⎥ yn ⎥ ⎥ 0 ⎥ 0 ⎥ ⎥ ⎥ 0 ⎦
［3 ］
） ，用于地质界面的模拟也是不合适的。离散数据集的单值曲面拟合插值法是较为
成熟、使用频率高的方法。与距离成反比的加权方法、径向基函数插值方法、样条曲面方法 等，拟合过程大都是全局性的，当增加、删除和修改数据点时都要重新计算，有些方法还必 须要求解维数很大的线性方程组，这就使得这些方法的效率和稳定性大大降低，寻求局部插 值拟合方法就成为今后地质界面数学模拟的发展目标。对于用来模拟复杂褶皱的多值曲面拟
1 前言
工程地质三维建模与可视化是应用计算机图形学和图像处理技术，是将工程地质勘测数 据和工程地质岩土体力学数值模拟分析的计算结果转换为图形图像在计算机屏幕上显示出 来并进行交互处理的理论、方法和技术。复杂地质体中的各类地质信息都可以被看作是三维 空间的函数，利用各种野外勘测数据分别建立相应的曲面拟合函数，进而利用计算机建立三 维地质模型，逼真反映地质结构全貌，达到直观地表达地质信息的分布规律、提高对于地质 规律的认识、指导地质工程项目的勘测施工的目的。因此，工程地质三维建模与可视化研究 有其重要的理论和现实意义。 工程地质三维建模与可视化研究中，地质界面的数学模拟是基础。由于地质界面必须精 确通过控制点（工程勘测数据点） ，通过离散的工程勘测数据点建立三维地质模型及其计算 机图形可视化显示属于重构问题，与计算机图形学中对于机械设计特别有效的形体的构造问 题是有根本区别的， 所以针对机械设计发展起来的自由曲线曲面造型技术
N
F ( x, y ) =
∑
i =1
Ai (
ri r r ln i 2 + 1 − i 2 ) 2 R R R
2
2
2
为使 (x 1 , y1 ) 点处 F(x,y)等于实测数据点值 f 1 ，需建立联立方程组：
4
工程地质计算机应用
2006 年 第 3 期
总 43 期
fi =
∑A
j =1
n
j
(
rij
W ( x, y ) = a0 + a1 x + a2 y + ∑ Fi ri 2 ln(ri 2 + ε )
i =1
பைடு நூலகம்
n
式中， a0，a1，a2，f i (i = 1, 2,..., n) 为待定系数；ri = (x-xi)2+(y-yi)2； ε 为调节曲率大小的经
2
验参数，一般取 10-2~10-6；n 为原始数据点个数。 n+3 个待定系数可通过下列方程组求得：
F ( x, y ) = ∑ a j ( x − x j ) 2 + ( y − y j ) 2 + c 2
j =1
n
[
]
1
1
2
, j = 1,2,..., n
其中c为常数，一般取1。将 n 个点 (x i , y i ) 的实测值 f i 代入上式建立联立方程：
f i = ∑ a j [( xi − x j ) 2 + ( y i − y j ) 2 + c 2 ] 2 , i = 1,2,..., n
n ⎧ W a a x a y Fij rij 2 ln(rij 2 + ε ) = + + + ∑ 0 1 j 2 j ⎪ j i =1 ⎪ n n n ⎪ Fij = ∑ xi Fi = ∑ yi Fi ⎨ ∑ i =1 i =1 i =1 ⎪ 2 2 ⎪ rij = ( x j − xi ) + ( y j − yi ) 2 ⎪ j = 1, 2,..., n ⎩
2
工程地质计算机应用
2006 年 第 3 期
总 43 期
合插值函数也是一个公认的难题。 地质体通常是不规则形体， 由多个各种成因类型形状各异的结构面
［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) 的距离成反 令：d k = 比[6]。 可表示为：
［1］
在地质层面模拟
［2 ］
中的应用受到极大限制。另外，现今热门的通过 Voronoi 图和 Dulauny 三角剖分
在空间构
造不规则三角网（Triangular Irregular Net, TIN）方法构造层面，一方面由于原始数据点一般 相距较远，常需要进行三角网的插值加密，另一方面，插值曲面不光滑，无法求出层面上某 点沿坐标轴的坡向、坡度、曲率和产状（TIN 方法无数学曲面的解析表达式，就无法求得对 于 x 或 y 的一阶和二阶偏导，而地质层面的走向、倾向和倾角又与曲面方程的偏导数有一定 的关系
j =1
n
求解此线性方程组可获得待定系数 aj (j=1,2,…,n)，将这些值代回插值函数式，即为通过 各实测数据点且处处连续光滑的曲面方程。 在数据点数量不太大的情况下，Multiquadric 法计算不太麻烦。这一方法提出的近 20 年 间，在水文测量、大地测量、地质及采矿、地球物理等领域得到广泛应用，效果良好。 2.1.3 基于平面弹性理论的曲面样条插值法 对地表地形或地下水位等离散数据进行函数拟合可借助平面弹性理论，把计算域视为无 限延伸平面，在平面 (x 1 , y1 ) 点上给定一垂直位移 A 1 时，平面各处将均会随着产生位移，若 则由 A 1 引起距离 (x 1 , y1 ) 点为 r(r<R) 处的点 限制距离 (x 1 , y1 ) 点半径为 R 以外位移为零时， (x,y)的位移 z A (r) 为：
T
" r1,n −12 ln(r1, n −12 + ε ) r1, n 2 ln(r1, n 2 + ε ) 1 x1 2 2 " r2,n −1 ln(r2, n −1 + ε ) r2, n 2 ln(r2, n 2 + ε ) 1 x2 " # # # # " 0 rn −1, n 2 ln(rn −1, n 2 + ε ) 1 xn −1 " rn −1,n 2 ln(rn −1,n 2 + ε ) " 1 " xn −1 " yn −1