正则化薄板样条函数拟合地层界面_王宝龙.pdf

第 45 卷 第 5 期 煤田地质与勘探 Vol. 45 No.5 2017 年 10 月 COAL GEOLOGY 2. Key Laboratory of Geo-Inatics of State Bureau of Surveying and Mapping, Chinese Academy of Surveying and Mapping, Beijing 100830, China; 3. College of Geosciences and Surveying Engineering, China University of Mining thin-plate spline function TPS; regularization; equilateral triangle mesh 三维地层界面构建是 3DGM三维地质建模、 3DGIS三维地理信息系统、3DCAD三维计算机辅 助设计的基础。地层界面的几何形态随三维空间位 置而变化，除非被断层切断，在大多数情况下是遵 ChaoXing 24 煤田地质与勘探 第 45 卷 循板状体受力弯曲的形式进行变形，形成复杂的曲 面，这与数学上薄板样条的力学机理十分相似。 由于技术与经济等多种因素的限制，地下采样 一般比较稀疏。由这些离散、稀疏的地下采样点， 尽可能合理地构建地层界面的几何形态，是三维地 质建模软件的基本功能。如果仅用这些原始采样点 进行三角剖分来构建地层界面往往比较粗糙，与真 实形态差别较大。为了尽可能地逼近真实地层界面 形态，提高其可视化效果，需要添加一些虚拟采样 点对地层界面形状进行细化描述，这些加密的虚拟 采样点坐标可以利用曲面拟合加以确定[1]。 通常情况下，地层界面的构建是根据竖直钻井 数据点构建的，数据点分布较为均匀，通过构建不 规则三角网就能构建出粗糙的地层界面。近年来， 在页岩气、煤层气、致密砂岩气的勘探开发中广泛 使用水平井钻井技术，这极大地提高了甲烷的采收 效益，但对常规的三维地质建模方法提出了新的问 题水平井数据在目的层中点的分布不像竖直井那 样相对均匀，往往一口水平井会多次穿过同一地层 界面，如果仍采用传统的不规则三角网剖分方法， 容易产生病态的三角网，因而需要寻找新的方法来 拟合地层界面。 笔者针对某页岩气勘探生产区水平钻井数据点 稀疏且分布不均的特点，采用等边三角形网格加密数 据点的方式来构建地层界面，加密的数据点可以用带 有正则化的薄板样条Thin Plate Spline简称 TPS函数 来确定。由于正则化参数与薄板样条弯曲能量之间负 相关，即正则化参数越大，弯曲能量越小，拟合的曲 面越趋于平缓，通过调节正则化参数的值，使拟合的 地质界面在可视化效果和精度上都能满足要求。 1 常用的曲面拟合方法 曲面的构建是三维建模的重要组成部分。曲面 的构建往往需要进行网格细分，从而得到一些加密 点，这些加密的数据点往往会取一些规则的格网的 顶点，如矩形格网和三角形格网。拟合这些加密点 的方法主要有距离反比权重法IDW、自然邻域法、 克里金法Kriging、 径向基函数法RBF和离散光滑 法DSI等。距离反比权重虽然简单，但是当采样点 的密度小时，无法达到满意的精度，同时也无法对 误差进行理论估计。胡进娟等[2]利用距离反比权重 实现了基于六面体空间数据的插值；周小平等[3]在 应用中采用二维 Voronoi 图进行邻域插值；颜辉武 等[4]研究了基于动态 Voronoi 图的距离反比权重的 改进研究，但仍有比较大的局限性。对于构造复杂 的地层界面且原始采样点稀疏的情况下，这些方法 在理论和应用上效果都不是很好。自然邻域法基于 区域大小按比例对原始样本应用权重来进行插值拟 合，效果上与距离反比权重相似。克里金法是根据 已知点采用预测最小误差方法来估计未知点的值， 同时可求出相应的最小估计方差，即能给出估计值 可能误差的数量指标，但计算量大，插值效果依赖 于克里金理论模型中变差函数的选取[5]，但地质模 型所遵循的变差函数是未知的，因此克里金方法的 人为因素主要体现在变差函数选取的合理与否。J L Mallet[6]提出了离散光滑插值，当数据非常密集时， 这种方法能够保持数值稳定且能够对曲线、曲面及 体的几何形态和属性进行描述等；周晓琴等[7]提出 一种基于离散光滑插值改进的高精度三维地质曲面 建模方法，通过综合利用测量、钻孔和物探等数据 及其相互间的相关性， 对地质体对象进行联合建模。 该方法不仅能拟合三维地表模型，给出拟合的置信 区间，而且能够构建出复杂的地质构造曲面，辅助 工程勘察设计等，但其操作方法复杂。径向基函数 能够拟合出光滑的三维表面，并有控制地生成表面 光滑程度的参数，薄板样条函数法作为径向基函数 的一种，能够根据离散、稀疏的控制点生成光滑的 三维表面，实现简单，对于数据贫乏区的陡坡问题， 有很多校正方法，如施加一个正则化参数[8]。Chen Chuanfa 等[9-11]提出了一种健壮的基于改进的正交 最小二乘薄板样条拟合方法，避免了普通 TPS 当两 个采样点距离过近时的病态问题，并将其应用于构 建大尺度的数字高程模型DEM，发现其效果好于 距离反比权重IDW和克里金Kriging等经典的拟 合方法；之后在此基础上提出了正则化的正交最小 二乘薄板样条法，避免了因控制点存在噪声而产生 过度拟合的问题，这种方法拟合的 DEM，在均方误 差和阴影渲染上都达到了令人满意的效果；为了能 更好地顾及地形特征，提出了基于地形特征点的 TPS 方法，运用于 LiDAR 点云数据构建 DEM，在 数据点尽可能少的情况下，很好地保持了 DEM 的 地形特征。 近年来， 有限元法已被广泛应用于拟合三 维表面， 其基本思想是将一个连续的“场”分隔成有限 个分区或单元， 通常是三角形或是矩形网格， 对每一 个单元给出一个近似解， 再将所有单元按照一定的方 式进行组合，来模拟或者逼近原始“场”，从而将一个 连续的无限问题转化成离散的有限问题求解。 陈娟等 [12-13]研究了三次样条插值构造四边形薄板单元及基 于面积坐标与 B 网方法的四边形样条单元，为有限 元法拟合复杂“场”，提供了有利的支持。 ChaoXing 第 5 期 王宝龙等 正则化薄板样条函数拟合地层界面 25 2 薄板样条函数基本原理 径向基函数RBF插值最早是由 Krige 于 1951 年提出的广泛应用于矿藏分析的 Kriging 方法。 1977 年，Duchon 从样条弯曲能量最小的理论出发，得 到了研究多元问题的薄板样条[14]。如今薄板样条函 数法被广泛应用于拟合高程、地下水位高度或污染 程度等三维表面。 TPS 生成一个通过控制点的表面，并使所有点 连接形成的所有坡面的斜度变化最小，即基于生成 最小曲率的面来拟合控制点。薄板样条函数的近似 表达式为式1。 2 1 , ln1,2,, p xyiii i f x yaa xa yw ddip    （）1 式中 x、 y 是待插值点 x、 y 坐标， di2x–xi2y–yi2， xi、yi是控制点 i 的 x、y 坐标。薄板样条函数包括 2 个部分aaxxayy 表示局部趋势函数，其表达式类 似于一阶线性趋势面；di2lndi表示基本函数，目的 是生成最小曲率的表面。参数系数 wi、a、ax和 ay 来源于线性方程组2线性方程组4。 2 1 ln p xyiiii i aa xa yw ddf    2 1 0 p i i w    3 11 0 pp i iii ii w xw y    4 式中 p 表示控制点的个数，fi表示控制点 i 的已知 值，系数的估算需要 p3 个方程，组成线性方程组， 然后解方程组而得。 其矩阵表达形式为 T          wv KP PO ao 5 2 ln ijijij ddK 6 11 22 3 1 1 1 p pp xy xy xy           P 7 3 3 000 000 000      O 8 11 22 13 113 1 0 0 0 ppx y pp vw a vw a a vw                     ，，，vowa9 薄板样条弯曲能量计算公式为 T f Iw Kw 10 3 正则化方法 正则化方法Regularization是求解不适定问题 的有效方法，起源于研究数学物理的反问题，以20 世纪60年代苏联科学院院士A N Tikhonov等[15]提 出处理不适定问题的正则化方法为标志。求解不适 定问题Kxy 的Tikhonov正则化方法的基本思想是 把正则化泛函数 Jax||Kx-y||2a||x||2的极小元 xa 作为Kxy 的正则化近似解，其中 a 称为正则化参 数。正则化方法的基本问题是，决定一种选择正则 化参数 a 的策略使得极小元收敛到原问题 Kxy 的 精确解。对式6增加一个正则化参数 a 得到式11。 22 ln ijijijij KddIa 11 其中 2 11 1 pp ij ij d p     12 这种情况下TPS得到的曲面不一定通过所有的 控制点，避免了控制点的噪声影响。 如果 a0就是普通的TPS，如果 a 为无穷大， TPS就退化为最小二乘法，此时弯曲量为0。当某 个区域内控制点数据贫乏时，若使TPS得到的曲面 经过控制点，就会产生陡坡现象，通过调整正则化 参数，可以避免陡坡现象的发生。 4 地层界面的拟合 在三维地质建模中，需要对地质体、地层、褶 皱、断层及矿体等三维自然地质现象进行描述。对 于地层体的构建离不开地层界面的构建，地层界面 是划分相邻地层的依据。对地层界面的描述，一般 采用不规则三角网TIN表达。但是，实际情况下， 原始采样点数据稀疏且分布不均匀，直接使用不规 则三角网， 会使一些三角形的边特别长， 形成“病态” 三角网。针对这种情况，本文采用加密的规则等边 三角形网格来生成虚拟的控制点数据。相比不规则 三角网，规则等边三角形网是最理想的Delaunay三 角网，三角形等角而且可以根据需要进行加密。 4.1 实验平台 本文实验是在“地质工程师三维助手”G3DA 上进行的，该软件是中国测绘科学研究院地理信息 工程重点实验室开发的面向三维地质建模的3DGIS 系统。G3DA软件是采用Visual Studio 2010、OpenGL 编程进行软件开发，功能上基本实现了三维地质模 型构建、编辑、可视化及分析。可用于地质矿山、 ChaoXing 26 煤田地质与勘探 第45卷 地下工程等相关领域的地下三维地质建模。其软件 主要是基于钻孔数据、等高线数据、二维勘探剖面 图等数据构建三维地质模型[16]。在G3DA中设计的 三维地质模型主要包括DataSource数据源点、线 要素等、DrillHole钻孔、Unconity不整合面、 Stratum Interface地 层 界 面、Fault断 层 面和 Stratum地层等。这些模型从整体上描述了地质体 的各种属性及地质体之间的拓扑关系。 4.2 拟合步骤及结果讨论 以某页岩气勘探区水平钻井资料为例，用图1 所示步骤对地层界面进行拟合。 图 1 地层界面拟合流程图 Fig.1 Flowchart of fitting the stratum interface a. 首先，从研究区原始钻井资料中提取几何 点、线等数据源要素，在G3DA中用这些数据源构 建该区域初始三维点数据源模型图2。从图中来 看，水平钻井的水平段，在目的层中延伸开来，相 比传统只有竖直井的情况下，大大提高了页岩气采 收率，节约了成本。但是，在目的层地层界面上， 相比竖直井，这些水平钻井数据点稀疏且分布不均 图3。 图 2 水平钻井数据空间分布 Fig.2 Spatial distribution of horizontal drilling data b. 以某一地层界面为例，在无边界线约束的条 件下，首先按照该层界面上的原始采样点的空间范 围生成规则的等边三角形网格加密点，三角形的边 长不同会对拟合结果产生一定的影响，在进行拟合 时，要选择适当的边长。如图3所示，生成的等边 三角形网格加密点，加密点之间的距离即三角边长 取1 000 m。 图 3 原始采样点空间分布及等边三角形加密点 Fig.3 Spatial distribution of original sampling points and equilateral triangle mesh vertices c. 加密点高程值由正则化薄板样条函数法拟 合确定。具体算法①根据式5式9和式11 式12，将原始采样点数据xi,yii1,2,⋯,p代入， 并设定正则化参数值，得到一个线性方程组；②解 线性方程组式5；③将加密点的x, y坐标代入式 1，求解加密点值。 d. 最后，根据拟合的数据点生成地层界面模 型，如图4所示。从上到下，从左到右，依次是不 同正则化参数下地层界面拟合的可视化效果。当正 则化参数 a0.0时，结果跟普通的薄板样条函数法 相同，曲面严格过原始采样点，由于原始采样点过 于集中在某一区域导致在没有采样点的区域形成非 常陡的坡，即“过冲”现象。当正则化参数 a0.001、 a0.01和 a0.1时，随着正则化参数的增大，拟合 的地层界面部分经过或是逼近原始采样点，一定程 度上避免了原始采样点噪声，且曲面越趋于平缓。 表1中列出了图4中不同正则化参数下相应的弯曲 能量及均方根误差的对比关系。从可视化的效果和 表中的对比来看，正则化参数与薄板样条函数的弯 曲能量负相关，即正则化参数越大，弯曲能量越小， 相应的与原始采样点的偏差越大，拟合的曲面越趋 于平缓。因此，可以选择适当的正则化参数来拟合 地层界面，以达到更好的拟合效果。 通过对比，综合考虑精度及地层界面平滑度， 确定正则化参数 a0.001，图5为最终G3DA构建 的目的层中所有地层界面模型。当不加地层边界线 ChaoXing 第5期 王宝龙等 正则化薄板样条函数拟合地层界面 27 图 4 不同正则化参数拟合结果对比 Fig.4 Fitting result of different regularization parameters 表 1 不同正则化参数拟合结果弯曲能量及均方误差 Table 1 Bending energy and mean square error of fitting results of different regularization parameters 正则化参数 a 弯曲能量 均方根误差 RMS 0.00 74.44 0.00 0.001 1.06 2.82 0.01 0.17 13.31 0.1 0.05 35.60 时图5a，则把采样点数据的最小外界矩形范围作 为地层界面的边界，增加地层边界线时则把矿区的 实际边界范围作为地层界面的边界，相比不附加地 层边界线的结果，增加地层边界线则忽略了大部分 不相关的区域。从结果上来看，对于稀疏数据点拟 合地质界面时，此方法可以达到较好的效果。 图 5 地层界面的拟合结果正则化参数 a0.001 Fig.5 Fitting result of stratum interfaces 5 结 语 针对原始采样点数据稀疏且分布不均问题，采 用正则化薄板样条函数法，通过编程实现地层界面 的拟合，通过对比不同正则化参数拟合地层界面的 可视化效果和误差， 得到一个较好的拟合地层界面。 正则化薄板样条函数不但可以拟合地层界面，而且 可以通过调整正则化参数使原始采样点与拟合的地 层界面接近，避免了数据点的噪声，保持地层界面 整体的光滑性，同时也避免了无数据点处地层界面 的“过冲”现象，对于构建三维地质曲面具有重要实 用意义。 参考文献 [1] 明镜. 三维地质建模技术研究[J]. 地理与地理信息科学， 2011，27415–18. MING Jing. A study on three-dimensional geological model- ing[J]. Geography and Geo-Ination Science，2011，274 15–18. [2] 胡进娟，盛业华. 基于规则六面体的三维空间数据插值方法[J]. 现代测绘，2005，28510–12. HU Jinjuan，SHENG Yehua. Three dimensional interpolation based on regular[J]. Modern Surveying and Mapping，2005， 28510–12. [3] 周小平，周瑞忠. 基于 Voronoi 图的新型几何插值及其与传统 代数插值方法的比较[J]. 岩石力学与工程学报，2005，241 133–138. ZHOU Xiaoping， ZHOU Ruizhong. New geometric interpolation based on voronoi diagram and the comparison with tradi- tional algebric [J]. Chinese Journal of Rock Mechanics and Engineering，2005，241133–138. [4] 颜辉武，祝国瑞，徐智勇. 基于动态 Voronoi 图的距离倒数加 权法的改进研究[J]. 武汉大学学报信息科学版，2004， 29111017–1020. YAN Huiwu， ZHU Guorui， XU Zhiyong. Improvement research on inverse distance weighting based on dynamic voronoi diagram[J]. Geomatics and Ination Science of Wuhan Uni- versity，2004，29111017–1020. [5] 武强，徐华. 虚拟地质建模与可视化[M]. 北京科学出版社， 2011. [6] MALLET J L. Geomodeling[M]. OxfordOxford University Press，2002. [7] 周晓琴， 王长海， 陈碧宇. 基于离散光滑插值的三维地质体构 造网格模型[J]. 深圳大学学报理工版， 2014， 314 396–400. ZHOU Xiaoqin，WANG Changhai，CHEN Biyu. 3D geological grid model based on discrete smooth interpolation[J]. Journal of Shenzhen University Science and Engineering，2014，314 396–400. 下转第 32 页 ChaoXing 32 煤田地质与勘探 第 45 卷 瓦斯突出模拟实验的一条可行途径。 参考文献 [1] 张玉贵，张子敏，曹运兴. 构造煤结构与瓦斯突出[J]. 煤炭学 报，2007，323281–284. ZHANG Yugui，ZHANG Zimin，CAO Yunxing. Deed coal structure and control to coal-gas outburst[J]. Journal of China Coal Society，2007，323281–284. [2] 刘明举， 龙威成， 刘彦伟. 构造煤对突出的控制作用及其临界 值的探讨[J]. 煤矿安全，2006，481045–50. LIU Mingju，LONG Weicheng，LIU Yanwei. Study on con- trolling action of tectonic coal to outburst and its critical value[J]. Safety in Coal Mines，2006，481045–50. [3] 李会军. 煤的变形行为及其产气机理研究以沁水盆地南部 为例[D]. 北京中国科学院，2011. [4] 张春华，刘泽功. 实验室煤与瓦斯突出模拟试验回顾及展望[J]. 中国安全科学学报，2011，21348–53. ZHANG Chunhua，LIU Zegong. Review and prospect on coal and gas outburst simulation experiment in laboratory[J]. China Safety Science Journal，2011，21348–53. [5] 尹光志， 赵洪宝， 许江， 等. 煤与瓦斯突出模拟实验研究[J]. 岩 石力学与工程学报，2009，2881674–1680. YIN Guangzhi， ZHAO Hongbao， XU Jiang， et al. Experimental study of simulation of coal and gas outburst[J]. Chinese Journal of Rock Mechanics and Engineering，2009，2881674–1680. [6] 颜爱华，徐涛. 煤与瓦斯突出的物理模拟和数值模拟研究[J]. 中国安全科学学报，2008，18937–42. YAN Aihua，XU Tao. Laboratory experiments and numerical simulation on coal and gas outbursts[J]. China Safety Science Journal，2008，18937–42. [7] 赵志刚，胡千庭， 耿延辉. 煤与瓦斯突出模拟实验系统的设计[J]. 矿业安全与环保，2009，3659–12. ZHAO Zhigang，HU Qianting，GENG Yanhui. Design of simulation test system of coal and gas outburst[J]. Mining Safety Environmental Protection，2009，3659–12. [8] 李涛，高凌薇，刘天林. 煤岩流变分析[C]//煤炭科学院地质勘 探分院文集第一集. 西安陕西人民出版社，1987. [9] 周建勋. 煤的变形与光性组构的高温高压变形实验研究及煤 田构造中石英的显微构造与组构[D]. 徐州中国矿业大学， 1991. [10] 姜波，秦勇，金法礼. 煤变形的高温高压实验研究[J]. 煤炭学 报，1997，22180–84. JIANG Bo，QIN Yong，JIN Fali. Coal deation test under high temperature and confining pressures[J]. Journal of China Coal Society，1997，22180–84. [11] HUANG He. Methane ation and retention[D]. Salt Lake City，Utah，USAUniversity of Utah，1992. 责任编辑 范章群 上接第 27 页 [8] DONATO G，BELONGIE S. Approximate thin plate spline mappings[J]. Lecture Notes in Computer Science，2002，2352 21–31. [9] CHEN Chuanfa，LI Yanyan. A robust of thin-plate spline and its application to DEM construction[J]. Computers and Geo- sciences，2012，4899–16. [10] CHEN Chuanfa， LI Yanyan， CAO Xuewei， et al. Smooth surface modeling of DEMs based on a regularized least squares of thin plate spline[J]. Mathematical Geosciences，2014，468 909–929. [11] CHEN Chuanfa，LI Yanyan. A thin plate spline-based fea- ture-preserving for reducing elevation points de- rived from LiDAR[J]. Remote Sensing， 2015， 79 11344– 11371. [12] 陈娟， 李崇君. 三次样条Hermite 插值基构造四边形薄板单元[J]. 中国科学数学，2015，4591523–1536. CHEN Juan，LI Chongjun. Area coordinates and B-net for quadrilateral spline elements[J]. Scientia Sinica Mathe- matica，2015，4591523–1536. [13] 陈娟，李崇君，陈万吉. 基于面积坐标与 B 网方法的四边形 样条单元[J]. 力学学报，2010，42186–90. CHEN Juan，LI Chongjun，CHEN Wanji. The cubic spline hermite interpolation bases for thin plate bending quadrilateral elements[J]. Chinese Journal of Theoretical and Applied Me- chanics，2010，42186–90. [14] DUCHON J. Splines minimizing rotation-invariant semi-norms in Sobolev spaces[M]. BerlinSpringer，197785–100. [15] TIKHONOV A N，ARSENIN V Y. Solutions of ill-posed prob- lems[M]. Washington，D CV H Winston Sons，1977. [16] 董前林，李青元，曹代勇，等. 多源地质数据综合三维建模及 地质学分析以青海木里三露天勘探区为例[J]. 煤田地质 与勘探，2017，45337–44. DONG Qianlin， LI Qingyuan， CAO Daiyong， et al. 3D modeling with integration of geological data of multi-sources and geologi- cal analysisTake the No.3 exploration area of Muli, Qinghai[J]. Coal Geology Exploration，2017，45337–44. 责任编辑 范章群 ChaoXing