复杂地质曲面三维插值—逼近拟合构造方法

李明超，缪正建，刘 菲，王 刚

(1．天津大学水利工程仿真与安全国家重点实验室，天津 300072;2．中国水电顾问集团成都勘测设计研究院，成都 610072)

1 前言

地质勘探结果大多反映在一些离散不规则分布的数据点上，为了通过这些离散数据建立起区域性连续的整体模型，需要利用插值和逼近的曲面处理方法［1，2］。曲面插值(surface interpolation)是严格通过给定的数据点来构造曲面，并根据原始数据点值来插补空白区的值;这类方法不改变原始数据点值。而曲面逼近(surface approximation)则是利用相对简单的数学曲面来近似构造复杂的地学曲面，根据一定的数学准则，使所给出的数学曲面最大限度地逼近地质曲面;通过拟合处理的曲面，原始数据点一般有所改变，所以曲面逼近的结果往往会取得平滑的效果。曲面的插值与逼近统称为曲面拟合(surface fitting)［3］。

在地质曲面构造中运用较多的插值和逼近方法包括按近点距离加权平均法、按方位取点加权法、双线性插值法、移动曲面插值法、二元三点插值法、Kriging插值法和三次样条函数拟合法、趋势面拟合法、加权最小二乘拟合法等［4～9］。插值方法能保证构造曲面严格通过原始数据点集合，适用于均匀分布的数据点，但外推能力和唯一性较差［1］;地质工程师一般倾向于使用插值方法，这在二维环境下能够有效地操作，但在三维环境下解决基于分布不均匀原始数据的复杂地质曲面构造时就会遇到较大的困难［6，9］。基于逼近思想的曲面拟合方法对原始数据点的分布没有任何要求，且外推能力和唯一性均较强，并能够保证曲面对原始数据点集合具有最佳逼近效果和很好的光顺性［10］;但是直接采用逼近方法应用于分布不均匀原始数据点的复杂曲面拟合也存在两点困难［3］，一是如何确定拟合曲面的控制网格顶点个数保证其满足精度要求，二是如何保证所构造的逼近曲面的形状能够满足实际需要。

由于自然界地质结构的复杂性，各种地质曲面的拟合构造方法始终是地质学家和计算机科学家关注的热点问题，目前仍然缺乏完善的处理方法［11］。文章以水利水电工程为例，在分析多源地质数据特征的基础上，针对曲面插值和逼近方法各自的优缺点，提出了插值—逼近相结合的复杂地质曲面拟合方法，以满足工程区域三维地质建模精度高且存储量小的要求［12］。

2 多源地质数据

在水利水电工程中，通过多种地质勘测方式如地质测绘、遥感、地质勘探等，得到了各种丰富多样的地质数据，如地质点资料，遥感图像，钻孔、平硐及物探信息等，如图1所示。同时，对这些原始数据加以解译分析，可获得相应的地层界线、断层、褶皱等构造迹线以及平面构造地质图、不同位置的横/纵剖面图和不同高程的平切图等。

图1 多源地质数据示例Fig．1 Examples of multi-source geological data

这些来源不同的数据在精度、分辨率、数量、质量等方面都存在较大的差异，为了使所有有效数据一起成为三维地质曲面拟合可利用的、可靠的、一致的信息，应采用统一的数据结构进行管理。因此，采用面向对象分类和系统的思想，将直接可用的钻孔、平硐信息与解译分析得到的剖面数据进行耦合处理［11］，按不同地层、断层等对象进行自动分层和三维化，最后统一以点集合的数据结构进行存储，供三维地质曲面拟合使用。图2为某水电工程一个地质界面的点集合数据实例，图中数据比较集中且分布均匀的区域是工程主要枢纽建筑物布置区域，而在其相关的周边地区勘探分布较少且离散。

3 地质曲面的插值—逼近拟合构造

图2 地层界面的多源数据集合Fig．2 Data set of a horizon boundary surface

根据上述分析可知，原始的工程地质勘探数据往往分布范围较大且不均匀，单独采用插值方法或者逼近方法都会遇到较大的困难，因此文章提出插值与逼近相结合的方法。在构造复杂地质曲面的过程中，充分考虑原始数据点的分布特点和曲面构造的精度要求，对于集中且均匀分布的原始数据，采用插值方法，使曲面严格通过这些数据点;对于分布离散不均匀的数据，若精度要求很高则采用插值方法进行加密，一般情况下采用逼近拟合方法，使曲面在给定精度下充分逼近原始数据。在水利水电工程坝区，)地质勘探数据集中且分布较均匀，其构造的地质曲面精度要求很高，这部分须采用插值方法;而在其相关的周边地区，对精度要求不高，采用逼近方法是可行的。

3．1 NURBS 曲面技术

NURBS(Non-uniform rational B-splines，非均匀有理样条曲线)技术是ISO(International Standard Organized，国际标准化组织)1991年颁布的STEP标准中自由型曲线曲面的唯一表示方法，它是在B样条函数基础上发展起来的，利用非均匀节点向量表达式来构造有理B样条，它对标准的解析图形(如圆锥曲线、二次曲面、旋转曲面等)和自由型曲线曲面提供了统一的数学描述。随着计算机辅助几何设计(CAGD)的发展，NURBS技术有了快速的发展，已被广泛应用于包括三维地质建模在内的工程实践中［13］。虽然如此，但由于NURBS技术的复杂性、地质对象的错综复杂和所具有的不规则形态，阻碍了NURBS在三维工程地质中更多的应用。实际上，针对复杂地质体形态的无规律性变化，选择NURBS技术进行水利水电工程三维地质建模与分析，具有节省存储空间、计算机处理简便易行、数据库管理方便、保证空间唯一性和几何不变性等优点，在地质表示方面有很高的应用价值。

给定一组(m+1)(n+1)的网格控制点Pij(0≤i≤m，0≤j≤n)的NURBS曲面可定义为:

式(1)中，wij为相应于控制点Pij的权因子;k、l为阶数;Nki(u)(0≤i≤m)和Nlj(v)(0≤j≤n)分别是定义在u、v向节点矢量U={u0，…，um+k|ui≤ui+1，i=0，…，(m+k － 1)}和 V={v0，…，vn+l|vj≤ vj+1，j=0，…，(n+l－1)}上的 k、l阶 B 样条基函数;Sij(u，v)表示拟合的曲面段，u∈［ui，ui+1］，i= (k － 1)，…，m ，v ∈ ［vj，vj+1］，j=(l－1)，…，n。在实际工程中，一般 k、l取3基本上就可以满足要求。

3．2 基于NURBS技术的插值—逼近拟合实现

以图2中给定的数据点集合为例，设该地层面的原始数据集合为 D={ps，s=1，2，…，m}，其中分布密集均匀的子集为 D1={pi，i=1，2，…，n}，剩余分布不均匀的子集为 D2={pj，j=1，2，…，r1}和 D3={pk，k=1，2，…，r2}，这里 n+r1+r2=m，如图3(a)所示。则基于NURBS技术的地质曲面插值—逼近拟合方法具体实现如下:

1)对于子集 D1，按照蒙皮法(skinning)［14］插值思想构造曲面，若曲面控制点为n，该方法的数学模型为:

式(2)中，X∈Rn，为曲面控制顶点构成的未知矢量;B∈Rn，为曲面集合D1及边界条件构成的已知矢量;A∈Rn×Rn，为系数矩阵。

由式(2)可知，蒙皮法是一种基于插值思想的曲面构造方法，能够保证曲面严格通过集合D1中的所有数据点。该方法使曲面在进行插值时顺序地通过一系列截面曲线，从而很好地生成复杂的自由曲面。在NURBS蒙皮过程中，首先要求各截面曲线的拓扑性质应该一致，即都是开曲线或都是闭曲线;其次在蒙皮插值过程中，插值的方式将影响最终的结果，根据工程地质实际需要，一般选择线性插值。图3(b)为子集D1通过蒙皮插值得到的曲面S1。

2)对于子集D2和D3，则直接利用NURBS曲面技术采用反算法［15］进行拟合。为了满足曲面的边界约束条件(邻接曲面边界和区域边界约束)，这里直接提取已构建曲面S1与子集D2、D3相邻的边界线加入子集D2、D3，同时根据研究区域将原始数据转化为一系列u、v方向的曲线矢量，利用反算法求出各曲线矢量上的控制点。进而，基于设定精度和控制点数据运用分段NURBS函数进行逼近拟合，获得按精度充分逼近子集D2、D3的拟合曲面S2、S3，如图3(c)所示。这里曲面控制点的权值均设置为1(可进行调整)，曲面边界公差为0．01，内部公差为0．1。

3)合并曲面S1和曲面S2、S3，由于构建S2、S3曲面的右边界和左边界数据均来自S1曲面，因此两者能够达到实际精度要求的缝合，最终得到完整的地质拟合曲面S，如图3(d)所示，该图对完整曲面进行了渲染表达，同时显示了相应的原始数据集，以便进行对比调整。

图3 基于NURBS技术的地质曲面插值—逼近拟合构造Fig．3 Interpolation-approximation fitting construction of geological surface based on NURBS technique

3．3 曲面分析与调整

在完成上述地质曲面的构造后，需要对其进行检查分析与调整，以满足实际精度要求和建模需要。这主要从以下两个方面进行:

1)地质结构合理性和曲面几何性检查。即检查所拟合的地质结构面整体变化趋势是否合理、在几何结构连续性及拓扑关系上是否正确。若发现不合理或错误之处，可快速方便地对其进行局部调整(NURBS曲面局部调整不会影响其他部位)或重新构造。

2)原始数据的精度分析。将所有原始数据与构造曲面进行偏差分析，以原始数据点与曲面间距离为变量，求其标准偏差σ(m):

式(3)中，n为数据点数;di为数据点i与构造曲面间的距离。

若σ小于控制误差ε，则曲面满足要求，否则需要对曲面偏差较大处按原始数据进行局部调整，直到达到误差要求。对于上述曲面，其原始数据点数n=597，控制误差 ε =0．5 m，根据式(3)求得标准偏差值σ=0．17 m，这主要是由逼近拟合曲面引起的误差，小于控制误差，精度已足以满足实际需要。

4 结语

文章针对工程地质多源数据的分布和地质曲面的特征，充分考虑了地质精度要求、曲面连续性和数据存储量等方面的均衡，提出并实现了基于NURBS技术的复杂地质曲面插值—逼近拟合构造方法。该方法对于工程关键区域集中且均匀分布的原始数据，采用NURBS蒙皮插值方法，使曲面严格通过这些数据点;对于周边区域分布离散的数据，采用逼近拟合方法，使曲面在给定精度下充分逼近原始数据;并可通过地质结构合理性检查、曲面几何性和原始数据精度分析对构造曲面进行检查分析和调整。实例研究表明，该方法所构造的地质曲面能满足地质工程师的实际需要，并能为进一步的三维地质建模提供基础，效果良好。

［1］黄健全，罗明高，胡雪涛．实用计算机地质制图［M］．北京:地质出版社，1998．

［2］Horsman J，Bethel W．Methods of constructing a 3D geological model from scatter data［C］//Proceedings of AVS 95，http://www － vis．lbl．gov/Publications/1995/Horsman － AVS95．ps，1995．

［3］彭芳瑜，周云飞，周 济．基于插值与逼近的复杂曲面拟合［J］．工程图学学报，2002，23(4):87－96．

［4］Fernández O．Obtaining a best fitting plane through 3D georeferenced data［J］．Journal of Structural Geology，2005，27(5):855－858

［5］吴春发，李 星．地质模拟中数据插值方法的应用［J］．地球信息科学，2004，6(2):50－52．

［6］周小文，付 晖，吴昌瑜．地层特性随机场插值方法应用研究［J］．岩土力学，2005，26(2):221－224．

［7］李晓军，王长虹，朱合华．Kriging插值方法在地层模型生成中的应用［J］．岩土力学，2009，30(1):157－162．

［8］Wang D，Wang H T，Xi J H，et al．Three-dimensional geological modeling with discrete smooth interpolation algorithm［C］//2010 International Conference on Remote Sensing，IEEE，2010:125－129．

［9］Sirakov N M，Granado I，Muge F H．Interpolation approach for 3D smooth reconstruction of subsurface objects［J］．Computers ＆Geosciences，2002，28(8):877 －885．

［10］秦绪佳，王 青，鲍虎军．基于散乱点的增量式曲面逼近［J］．计算机辅助设计与图形学学报，2006，18(9):1408－1413．

［11］徐 华，武 强．复杂地质体中多值面的网格生成算法［J］．计算机辅助设计与图形学学报，2002，14(7):609－612．

［12］钟登华，李明超．水利水电工程地质三维建模与分析理论及实践［M］．北京:中国水利水电出版社，2006．

［13］Zhong D H，Li M C，Song L G，et al．Enhanced NURBS modeling and visualization for large 3D geoengineering applications:an example from the Jinping first-level hydropower engineering project，China［J］．Computers＆ Geosciences，2006，32(9):1270－1282．

［14］Piegl L，Tiller W．Algorithm for approximate NURBS skinning［J］．Computer－Aided Design，1996，28(9):699－706．

［15］施法中．计算机辅助几何设计与非均匀有理 B样条(CAGD＆NURBS)［M］．北京:高等教育出版社，2001．