CN102147934A - 基于tin的地质界面三维形态分析方法 - Google Patents

Abstract

本发明公开了一种基于TIN的地质界面三维形态分析方法，包括在使用三维TIN模型来模拟地质界面的基础上，使用空间几何和计算机图形学原理计算地质界面的一般几何形态参数(坡度、夹角等)和距离场，利用空间插值和趋势剩余分析技术，分级提取出地质界面的形态趋势以及形态起伏，最后，将地质界面三维形态分析结果应用于控矿地质因素场的模拟和隐伏矿体预测。本发明能够精确、高效地分析和提取地质界面的各种形态参数，定量地表达控矿地质因素，在隐伏矿体立体定量预测中具有重要意义。

Description

基于TIN的地质界面三维形态分析方法

技术领域

本发明涉及一种基于TIN的地质界面三维形态分析方法。TIN指Triangulated IrregularNetwork，不规则三角网。

背景技术

传统GIS空间分析方法无法有效的获取和表达控矿地质因素。而结合数学形态学提出的基于三维地质体栅格模型的地质体三维形态分析方法(邓浩面向隐伏矿体预测的三维地质建模与空间分析若干技术研究长沙：中南大学2008)(毛先成，唐艳华，邓浩地质体的三维形态分析与方法中南大学学报(自然科学版)，2011)，存在如下缺点：栅格数据模型空间数据量大、低空间精度，特别是对边界表示时不能准确表达地质体的边界，且生成速度、存储空间和表达精度难以平衡；数学形态学球形结构元素半径的选取在量化方面需要改进；不适合对薄片状地质体或存在狭窄部分的地质体进行趋势形态分析，当结构元素半径较大时，趋势形态会出现断裂的错误，这是因为地质体的狭窄部位会因滤波运算而作为局部起伏剔除掉，从而发生断裂。

发明内容

本发明所要解决的技术问题是提出一种基于TIN的地质界面三维形态分析方法，该基于TIN的地质界面三维形态分析方法能够更精确、更高效的提取地质界面的各种形态参数，有效的定量化表达控矿地质因素。

本发明的技术解决方案如下：

一种基于TIN的地质界面三维形态分析方法，包括以下步骤：在使用三维TIN模型来模拟地质界面的基础上，使用空间几何和计算机图形学原理计算地质界面的一般几何形态参数和距离场；利用空间插值和趋势剩余分析技术分级提取出地质界面的形态趋势以及形态起伏。

所述的一般几何形态参数包括坡度和夹角；

(1)坡度的计算方法为：在三维TIN模型中，任一个三角形ABC(平面方程为：z＝a₀+a₁x+a₂y)的三个顶点的坐标为A(x₁，y₁，z₁)，B(x₂，y₂，z₂)，C(x₃，y₃，z₃)，三角形ABC的坡度为：

$S=\arccos \frac{1}{\sqrt{a_1^2+a_2^2+1}}$

其中， $a_1=\frac{(y_2-y_1)*(z_3-z_1)-(y_3-y_1)*(z_2-z_1)}{(x_3-x_1)*(y_2-y_1)-(x_2-x_1)*(y_3-y_1)},$ $a_2=\frac{(x_3-x_1)*(z_2-z_1)-(x_2-x_1)*(z_3-z_1)}{(x_3-x_1)*(y_2-y_1)-(x_2-x_1)*(y_3-y_1)};$

(2)夹角的计算方法为：

所述的夹角指两个地层界面交线上的某点的夹角，即经过该点的两个三角面片之间的夹角，方法如下：

首先将地质界面栅格化以形成多个块体单元即单元，即采用三维规则网格将研究区地质空间划分为许多小立方体，对两个地质界面TIN模型中穿过相同块体单元的三角面片两两求夹角，最后将求得的夹角的平均值作为地质界面之间的夹角α^*，夹角α^*的求法如下：

设两个三角形的法向量分为V₁(x_v1，y_v1，z_v1)、V₂(x_v2，y_v2，z_v2)，两法向量的夹角为α，即两个三角面片之间的夹角，计算的公式为：

$\cos \mathrm{\α}=\frac{V_1\·V_2}{\left|V_1\right|\left|V_2\right|}=\frac{x_{v1}*x_{v2}+y_{v1}*y_{v2}+z_{v1}*z_{v2}}{\sqrt{{x_{v1}}^2+{y_{v1}}^2+{z_{v1}}^2}\sqrt{{x_{v2}}^2+{y_{v2}}^2+{z_{v2}}^2}}$

$\mathrm{\α}=\arccos \frac{x_{v1}*x_{v2}+y_{v1}*y_{v2}+z_{v1}*z_{v2}}{\sqrt{{x_{v1}}^2+{y_{v1}}^2+{z_{v1}}^2}\sqrt{{x_{v2}}^2+{y_{v2}}^2+{z_{v2}}^2}}$

则地质界面之间的夹角α^*的表达式为：

其中n等于两个地质界面在同一块体单元中的三角面片的夹角个数。

一些成矿因素受到地质界面的夹角和坡度大小的控制，对地质界面夹角和坡度等一般几何形态参数进行提取可以定量表达这种控矿因素，对于发掘有利成矿信息具有重要作用。-这一段后续由我放到说明书中去。

距离场的计算方法为：

控矿地质因素场与空间中某点到相关联的地质界面的距离有关，即控矿地质因素场是空间中某点到地质界面距离的空间分布函数，

单元到地质界面的最近距离的求解，实际上是寻找单元中心点到地质界面TIN距离集合的最小值；TIN模型是由许多三角面片组成，所以点到地质界面的距离实际上是点到所有三角片面的最短距离；

求出预测空间即前述的研究区地质空间中每个单元对应的距离值，建立距离场模型，所述的距离值指某一个单元的中心点到TIN模型中各个三角面片的距离值的最小值【一个单元只有一个距离值】。

分级提取出地质界面的形态趋势以及形态起伏包括：

空间趋势形态反映的是空间物体在空间区域上的主体特征，它忽略了局部形态起伏以揭示空间区域上的主体特征。在矿床地质空间中，相关地质体的主体形态将对周围成矿起到一定因素的影响，起伏则是地质异常作用的侧面反映。因此，对地质体形态(地质体形态即该地质体对应的地质界面的形态)的分析将是控矿因素定量提取的重要步骤。

所述的地质界面的形态趋势是指对地质界面TIN模型中的三角面片的顶点的Z值按照距离平方反比法重新赋值得到的新顶点建立的新的TIN模型，新的TIN模型称为地质界面的趋势面，地质界面的趋势面的形成过程如下：

①确定搜索距离d(d的大小依据勘探线工程间距确定，假如勘探线工程间距为100米，则一级趋势分析的d取100米，二级趋势分析的d取2倍工程间距，则为200米)，只有与估测点距离小于d的样品点(用TIN中三角形的顶点来描述实际的采样数据点，此处样品点即实际的采样数据点，也是三角形的顶点)才能参与计算，其中距离是指两点在XY平面上投影点间的距离，而不是两点的绝对距离；估测点指原来的地质界面TIN模型中的任意一个顶点在XY平面上的投影点；建立新的TIN模型的过程，就是将所有的估测点都计算出新的z值，形成新的TIN模型的顶点，再根据该顶点构建新的TIN模型。

②利用距离平方反比法重新计算样品点的所有z值，公式如下：

$z=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{z_i}{d_i^2}}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{d_i^2}}$

其中，z_i表示该估测点周围、与估测点距离小于d的样品点的垂直坐标分量；d_i表示该估测点与样品点的距离；此处距离都指两点在XY平面上投影点间的距离；

③利用计算的z值，建立新的地质界面TIN模型(新的地质界面TIN模型采用原TIN模型顶点集的x、y坐标，与新计算的相应的z坐标做成新的顶点，边集不变，拓扑结构不变)，形成地质界面的趋势面；

按照上述方法，采用不同的搜索距离d(滤波级别一般为2级，距离可分别取1d和2d)，即可分级提取得到地质界面的形态趋势；

所述的形态起伏是指地质界面相对于形态趋势面的起伏幅度(也称为波幅)，即地质界面中的z值减去形态趋势面的z值，即可得到形态起伏。由于地质界面的形态趋势的提取是分级的，所以，形态起伏也是分级的。

通过前述的方法，得到了地质界面的一般几何形态参数(坡度、夹角等)、距离场、地质界面的分级形态趋势以及分级形态起伏等三大方面的控矿地质因素指标。

有益效果：

本发明公开了一种基于TIN(Triangulated Irregular Network，不规则三角网)的地质界面三维形态分析方法，包括在使用三维TIN模型来模拟地质界面的基础上，使用空间几何和计算机图形学原理计算地质界面的一般几何形态参数(坡度、夹角等)和距离场，利用空间插值和趋势剩余分析技术，分级提取出地质界面的形态趋势以及形态起伏，最后，将地质界面三维形态分析结果应用于控矿地质因素场的模拟和隐伏矿体预测。本发明是基于TIN模型的，TIN模型具有适应任意形态曲面、表达精度可自由调节、适合于稀疏数据、建模速度快等优点，因此，本发明可以发挥TIN模型的优点，能够精确、高效地分析和提取地质界面的各种形态参数，定量地表达控矿地质因素，为隐伏矿体立体定量预测提供成矿信息依据。

附图说明

图1为基于TIN的地质界面三维形态分析流程图；

图2为不整合面坡度因素(gU)分布栅格模型示意图(-150米至-200米标高范围)；

图3为地层界面与不整合面的夹角因素(aU_S)分布栅格模型示意图(-150米至-200米标高范围)；

图4为不整合面距离因素(dU)分布栅格模型示意图(-150米至-200米标高范围)

图5为不整合面原模型和趋势形态模型示意图，(a)不整合面的三维TIN模型示意图；(b)插值搜索半径为100米的不整合面一级趋势TIN模型示意图；

图6为不整合面形态起伏分级提取示意图，(a)一级外凸内凹部分；(b)二级外凸内凹部分；

图7为不整合面形态起伏因素(waU和wbU)分布栅格模型示意图(-150米至-200米标高范围)，(a)一级起伏waU；(b)二级起伏wbU。

具体实施方式

以下将结合附图和具体实施实例对本发明做进一步详细说明：

一种基于TIN的地质界面三维形态分析方法，包括：在使用三维TIN模型来模拟地质界面的基础上，使用空间几何和计算机图形学原理计算地质界面的一般几何形态参数(坡度、夹角等)和距离场，利用空间插值和趋势剩余分析技术，分级提取出地质界面的形态趋势以及形态起伏，最后以栅格模型来表达控矿地质因素场，实现各种控矿地质因素的定量模拟。1、根据权利要求1所述的基于TIN的地质界面三维形态分析方法，其特征在于，所述使用空间几何和计算机图形学原理计算地质界面的一般几何形态参数(坡度、夹角等)包括：

(1)坡度提取

由于坡度的计算都是基于地质界面的TIN模型，因此对于给定点的坡度计算实质上求该点所对应的三角面片的坡度。具体过程如下：

TIN模型上任意三角形ABC可以用如下平面方程表示：

z＝a₀+a₁x+a₂y    (1)

平面上坡度处处相等，可以用下式计算该三角形的坡度：

$\text{S=arccos}\frac{1}{\sqrt{a_1^2+a_2^2+1}}---\left(2\right)$

因此只要计算出a₁，a₂就能计算出坡度。

设三角形ABC的三个顶点的坐标为A(x₁，y₁，z₁)，B(x₂，y₂，z₂)，C(x₃，y₃，z₃)，则垂直于该三角形平面的法向量为：

设三角形平面上有一个点M(x，y，z)，则

由公式(3)可以得出三角形所在平面的方程：

$\left|\begin{array}{ccc}x-x_1& y-y_1& z-z_1\\ x_2-x_1& y_2-y_1& z_2-z_1\\ x_3-x_1& y_3-y_1& z_3-z_1\end{array}\right|=0---\left(4\right)$

结合公式(1)和公式(4)可以得出：

$a_1=\frac{(y_2-y_1)*(z_3-z_1)-(y_3-y_1)*(z_2-z_1)}{(x_3-x_1)*(y_2-y_1)-(x_2-x_1)*(y_3-y_1)}---\left(5\right)$

$a_2=\frac{(x_3-x_1)*(z_2-z_1)-(x_2-x_1)*(z_3-z_1)}{(x_3-x_1)*(y_2-y_1)-(x_2-x_1)*(y_3-y_1)}---\left(6\right)$

最后代入公式(2)就可以求出S。

求得的S就是对应于每个三角面片的坡度信息。

(2)地质界面间夹角提取

由于地质界面都是采用TIN模型进行建模，所以求两个地层界面交线上的某点的夹角就是求经过该点的两个三角面片之间的夹角。该方法首先将地质界面栅格化，然后再求相交的立方体单元，最后求地质界面之间的夹角。

地质界面栅格化：

地质界面是以TIN模型表达的，其栅格化的实质就是将地质界面TIN穿过的块体单元提取出来，形成地质界面的块体模型。具体步骤如下：

①三角面片投影。首先将三角面片的三个顶点向三个坐标平面：XY平面、XZ平面和YZ平面进行投影，得到三个投影三角形。下面以XY平面为例，设投影后的三角形为ABC。

②获取最小外包矩形。根据投影面的三角形的三个顶点计算三角形的最小外包矩形EFGH，该外包矩形与传统的外包矩形不同，必须与矿床块体模型在相应平面投影的网格线相重合。外包矩形的目的在于减少搜索范围，提高计算效率。

③计算投影面上三角形轨迹网格。采用计算机图形学中的Liang Barsky参数化线段裁剪多边形算法，计算三角形的三条边穿过的平面网格。但是，当三角形投影为一直线段，而且也与矿床块体模型在相应平面投影的网格线重合时，则认为平面上三角形的边穿过的网格为该投影直线段左右网格(或者上下网格)。

④填充三角形。采用扫描线算法计算XY平面上三角形ABC内部的网格索引(i_x，j_y)。

⑤三个投影面的网格索引求交。与XY平面上的处理类似，按照上面的步骤分别求得TIN模型中三角面片在XZ平面和YZ平面上的投影，并计算出相应的内部网格索引，然后将三个投影面上的网格索引的分量求交，最后求得三角面片通过的块体单元的索引，并且将穿过块体单元的三角形的ID记录下来。

⑥通过对地质界面TIN模型中的所有三角面片循环计算，可以求出地质界面通过的块体单元(采用三维规则网格将研究区地质空间划分为许多小立方体，小立方体称为立体单元，或称单元、体元、块体。)，最后生成地质界面块体模型。

夹角的计算：

首先将两地层界面的块体模型进行求交，其实质就是计算块体模型中重合的块体单元的索引值，同时记录穿过重合块体单元的三角面片的标识ID值，然后计算穿过同一块体单元的三角面片之间的夹角。其中，需要对一种情况进行特殊处理：当相交处的块体单元为三角面片的顶点时，地质界面的TIN模型就会有多个三角面片同时穿过相应的块体单元，这种情况下，对两个地质界面TIN模型中穿过相同块体单元的三角面片两两求夹角，最后将这些夹角的平均值作为最后结果。

两个三角形的夹角计算经常采用向量法，首先规定三角形的法向量向上，即法向量的z分量大于等于0，如果z小于0，则所有分量取反；然后求两个三角形的法向量，即三角形的两条边的向量的叉积。法向量的夹角即为两个平面的夹角。

设两个三角形的法向量分为V₁(x_v1，y_v1，z_v1)、V₂(x_v2，y_v2，z_v2)，两法向量的夹角为α，则：

$\cos \mathrm{\α}=\frac{V_1\·V_2}{\left|V_1\right|\left|V_2\right|}=\frac{x_{v1}*x_{v2}+y_{v1}*y_{v2}+z_{v1}*z_{v2}}{\sqrt{{x_{v1}}^2+{y_{v1}}^2+{z_{v1}}^2}\sqrt{{x_{v2}}^2+{y_{v2}}^2+{z_{v2}}^2}}---\left(7\right)$

$\mathrm{\α}=\arccos \frac{x_{v1}*x_{v2}+y_{v1}*y_{v2}+z_{v1}*z_{v2}}{\sqrt{{x_{v1}}^2+{y_{v1}}^2+{z_{v1}}^2}\sqrt{{x_{v2}}^2+{y_{v2}}^2+{z_{v2}}^2}}---\left(8\right)$

则地质界面之间的夹角为α^*(其中n等于两个地质界面在同一块体单元中的三角面片的夹角个数)：

${\mathrm{\α}}^{*}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_i}{n}$

夹角的赋值：

通过上一步骤，获得了两地层界面相交处的块体集合，并且该集合中的块体包含了两者的夹角信息。因此，对于整个地质空间中的立体单元就通过求取与上述块体集合中的块体的最近距离来获取相应的夹角信息。

因为获得的相交块体集合数量有限，所以采用暴力算法就能快速求取地质空间中的立体单元中心点到该块体集合中块体单元中心点的最短距离。求得的块体集合中最短距离所对应的块体的夹角信息就是该立体单元所需要赋值的夹角。

距离场的计算方法为：

控矿地质因素场与空间中某点到相关联的地质界面的距离有关，即控矿地质因素场是到地质界面距离的空间分布函数。在地质空间中，本研究选择欧式距离作为空间距离。地质界面之间的距离或地质空间中某点到地质界面的距离，用以表示和研究地质界面之间的几何接近程度或地质界面对空间中某点的影响程度。在研究中，我们将点到地质界面的距离约定为点到地质界面的最小距离，即预测空间中某单元(体元)到地质界面的最近距离作为控矿地质因素场对单元的影响程度。

按前述，地质界面可以采用三维矢量模型即三维TIN模型来描述，因此，单元到地质界面的最近距离的求解，实际上是寻找单元中心点到地质界面TIN距离集合的最小值。求任意一个空间点与任一实体距离的计算，最简单的方法是对地质空间中所有点(所有点指研究区地质空间中单元的中心点)进行距离量算，找出其中距离最小的点。TIN模型是由许多三角面片组成，所以点到地质界面的距离实际上是点到所有三角片面的最短距离。从数学上讲，实际上就是求点到面的最短距离，但是由于实际情况，需要考虑以下两种情况：①点到三角面的垂足落在三角面内；②点到三角面的垂足在三角面外。因此，点到地质界面的首要问题是要判断垂足是否在三角形内。

设空间中有一个点P和一个三角形ABC，设P点的坐标为(x₀，y₀，z₀)，A点坐标为(x₁，y₁，z₁)，B点坐标为(x₂，y₂，z₂)，C点的坐标(x₃，y₃，z₃)。三角形ABC所在平面的法向量(l，m，n)可以由下式求出：

$l=\left|\begin{array}{ccc}1& y_1& z_1\\ 1& y_2& z_2\\ 1& y_3& z_3\end{array}\right|=\left|\begin{array}{cc}{y}_2-y_1& z_2-z_1\\ y_3-y_1& z_3-z_1\end{array}\right|=(y_2-y_1)(z_3-z_1)-(y_3-y_1)(z_2-z_1)---\left(9\right)$

$m=\left|\begin{array}{ccc}{x}_1& 1& z_1\\ x_2& 1& z_2\\ x_3& 1& z_3\end{array}\right|=\left|\begin{array}{cc}{z}_2-z_1& x_2-x_1\\ z_3-z_1& x_3-x_1\end{array}\right|=(z_2-z_1)(x_3-x_1)-(z_3-z_1)(x_2-x_1)---\left(10\right)$

$n=\left|\begin{array}{ccc}{x}_1& y_1& 1\\ x_2& y_2& 1\\ x_3& y_3& 1\end{array}\right|=\left|\begin{array}{cc}{x}_2-x_1& y_2-y_1\\ x_3-x_1& y_3-y_1\end{array}\right|=(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)---\left(11\right)$

三角形ABC所在平面的方程为

lx+my+nz＝lx₁+my₁+nz₁    (12)

通过P(x₀，y₀，z₀)点垂直于三角形ABC所在平面的直线的参数形式方程为

$\left\{\begin{array}{c}x=x_0+\mathrm{tl}\\ y=y_0+\mathrm{tm}\\ z=z_0+\mathrm{tn}\end{array}\right.---\left(13\right)$

用上式带入ABC所在平面的方程，可以解得

$t=\frac{l(x_1-x_0)+m(y_1-y_0)+n(z_1-z_0)}{l^2+m^2+n^2}---\left(14\right)$

再将这个t代入从P点向三角形ABC所在平面所作的垂线方程，就可以求得垂足Q的坐标(X，Y，Z)如下：

$\left\{\begin{array}{c}X=x_0+\mathrm{tl}\\ Y=y_0+\mathrm{tm}\\ Z=z_0+\mathrm{tn}\end{array}\right.---\left(15\right)$

下面就是判断垂足Q是否落在三角形ABC的三边围成的闭区域内。l，m，n中至少有一个不等于0，不妨设n≠0，依次计算

$\left|\begin{array}{ccc}X& Y& 1\\ x_2& y_2& 1\\ x_3& y_3& 1\end{array}\right|,\left|\begin{array}{ccc}{x}_1& y_1& 1\\ X& Y& 1\\ x_3& y_3& 1\end{array}\right|,\left|\begin{array}{ccc}{x}_1& y_1& 1\\ x_2& y_2& 1\\ X& Y& 1\end{array}\right|$

如果这三个值的绝对值加起来正好等于n的绝对值，就说明垂足Q落在三角形ABC的三边围成的闭区域内。这时，P点到三角形ABC的最短距离就是P点到三角形ABC所在平面的垂直距离，也就是线段PQ的长度

$\mathrm{PQ}=\sqrt{{(X-x_0)}^2+{(Y-y_0)}^2+{(Z-z_0)}^2}---\left(16\right)$

如果

的绝对值之和不等于n的绝对值，说明垂足Q不落在三角形ABC三边围成的闭区域内。这时，分别求P点到线段AB、BC和CA的最短距离，然后取其中最短者。以求P点线段AB的最短距离为例，具体过程如下：

①判断角PAB是否为钝角或为直角，方法为判断向量AP与向量AB的点积是否小于等于0，即若满足此条件，则P到线段AB最短距离为线段PA。

$\mathrm{PA}=\sqrt{{(x_1-x_0)}^2+{(y_1-y_0)}^2+{(z_1-z_0)}^2}---\left(17\right)$

如果不满足则转到下面第②步。

②判断角PAB是否为钝角或者直角，过程①相同。

③当①和②条件不满足，则P点到线段AB的垂足位于线段AB内部，设垂足Q(X，Y，Z)，直线AB的参数方程可写为：

$\left\{\begin{array}{c}x=t(x_2-x_1)+x_1\\ y=t(y_2-y_1)+y_1\\ z=t(z_2-z_1)+z_1\end{array}\right.---\left(18\right)$

则向量PQ可表示为：(x₁+t(x₂-x₁)-x₀，y₁+t(y₂-y₁)-y₀，z₁+y(z₂-z₁)-z₀)

向量AB为：(x₂-x₁，y₂-y₁，z₂-z₁)

由于向量AB与向量PQ垂直，则因此可以求得

$t=\frac{(x_1-x_2)*(x_1-x_0)+(y_1-y_2)*(y_1-y_0)+(z_1-z_2)*(z_1-z_0)}{(x_1-x_2)*(x_1-x_2)+(y_1-y_2)*(y_1-y_2)+(z_1-z_2)*(z_1-z_2)}---\left(19\right)$

将参数t代入到公式(18)，就可以求出垂足点Q的坐标，然后求出PQ的长度。

通过这种点到TIN的距离算法编程可计算出地质空间中每一个单元，即立体单元中心点到地质界面的最近距离，利用该距离在地质空间中的分布实现了控矿地质因素场的离散化建模。

所述利用空间插值和趋势剩余分析技术，分级提取出地质界面的形态趋势以及形态起伏包括：

借鉴传统的曲面的趋势形态分析方法，结合地质界面的实际情况，采用地质界面的原始TIN模型，利用距离平方反比法对地质界面进行趋势形态分析。

地质界面TIN模型的趋势形态分析实际上就是利用估测点周围一定距离之内的样品点对估测点的z值重新赋值，然后利用计算的新z值生成地质界面的趋势面。具体过程如下：

①确定搜索距离d，只有与估测点距离小于d的样品点才能参与计算，其中距离是指两点在XY平面上的相对距离，而不是两点的绝对距离。

②利用距离平方反比法重新计算样品点的所有z值，公式如下：

$z=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{z_i}{d_i^2}}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{d_i^2}}---\left(20\right)$

③利用计算的z值，建立新的地质界面的TIN模型，形成地质界面的趋势面。

地质体起伏形态可以通过地质界面的波状起伏来描述。这种波状起伏可以通过地质界面相对于其趋势平面的起伏幅度来表达。这种起伏幅度又叫波幅，反映地质界面的弯曲程度和在空间中的波状起伏幅度变化。

这种波状起伏的建模可以通过前文所述的趋势形态分析来实现，即将地质界面标高z的观测值的变化分解为二部分：①观测值总体变化重建的一个曲面，即趋势面T(x，y)；②观测值的局部变化，称谓为信号，即剩余趋势面R(x，y)；即

z(x，y)＝T(x，y)+R(x，y)    (21)

则剩余趋势面R(x，y)：

R(x，y)＝z(x，y)-T(x，y)    (22)

对于地质曲面上的的一个点z(x₀，y₀)，可通过R(x，y)绝对值来描述该点所处的局部的形态起伏程度，通过R(x，y)的符号来描述该点所处的局部是外凸或是内凹，还是平坦，有：

对于R(x，y)所表示的局部起伏，可视为不同级别起伏的叠加，有：

R(x，y)＝R₁(x，y)+R₂(x，y)+...+R_n(x，y)    (24)

其中，R₁(x，y)，R₂(x，y).，..，R_n(x，y)表示了n级局部起伏变化，各级局部变化的筛分可用趋势剩余分析方法。

利用上述方法可实现形态起伏的分级提取。因为插值搜索半径d决定可滤除的波形的大小，所以通过改变d值可提取出不同级别的起伏。算法首先用d值较小的插值搜索半径对地质体进行形态滤波，提取出第一级趋势形态和较小的起伏；接着增加d值，对第一级趋势形态再次形态滤波，进一步得出第二级趋势形态和较上一级别起伏更大的新一级起伏；再次增加d值，对上一级别趋势形态进行新的一次形态滤波，得出新一级趋势形态和起伏；这样一直迭代下去，直到d值达到指定最大阈值为止。d值改变，趋势就改变；趋势改变，相应的起伏就改变了

地质界面的剩余趋势面通过波幅很好地反映了地质界面在空间中的起伏变化程度，为了定量描述地质界面上任何一点的波幅，需要计算地质界面与趋势面上具有相同x，y值的点的z坐标的差值，实质上是通过x，y计算地质界面和趋势面的z值。计算过程如下：

地质界面TIN模型中任意三角形ABC，已知A(x₁，y₁，z₁)，B(x₂，y₂，z₂)，C(x₃，y₃，z₃)，三角形ABC所在平面中的任意点M(x，y，z)，平面方程为：

$\left|\begin{array}{ccc}x-x_1& y-y_1& z-z_1\\ x_2-x_1& y_2-y_1& z_2-z_1\\ x_3-x_1& y_3-y_1& z_3-z_1\end{array}\right|=0---\left(25\right)$

解得： $z=\frac{(z_2-z_1)(y_3-y_1)(x-x_1)+(y-y_1)(x_2-x_1)(z_3-z_1)-(x-x_1)(y_2-y_1)(z_3-z_1)-(y-y_1)(z_2-z_1)(x_3-x_1)}{(x_2-x_1)(y_3-y_1)-(y_2-y_1)(x_3-x_1)}+z_1---\left(26\right)$

但是，当三角形ABC垂直于XY平面时，相同的x，y值会有无数个z值对应，对于此种情况，将z确定为最小值与最大值和的一半，即：

$z=\frac{z_{\max }+z_{\min }}{2}---\left(27\right)$

实施例1：

以下六个步骤描述基于TIN的地质界面形态分析方法实施，以福建省尤溪县丁家山铅锌矿床的不整合面为实例。

(1)建立不整合面TIN模型和矿化空间栅格模型：利用Datamine软件把地质剖面上不整合面(当上下两套不同时代地层之间出现过沉积间断或地层缺失时，这两套地层之间的接触界面称为不整合面)的已知位置点(已知位置点是指利用地表填图、勘探工程、井下巷道等技术手段获取的不整合面位置的控制点)用线连接起来，形成一系列线段，然后利用线框模型建模功能建立不整合面的TIN模型。TIN模型把这些线段以三角面片的形式拼接起来形成一个三角形网格来模拟地质界面，并导出到Access数据库。采用精度为10m×10m×10m的立体单元格，划分矿化空间，形成矿化空间的栅格模型。栅格模型存储在Access数据库中，同时，三维形态分析结果和控矿因素指标提取结果也存储在该数据库中。

(2)地质界面三维形态分析软件TriangelCal开发：采用C#语言，在VisualStudio 2005平台下，利用本发明内容中所述的具体实施方式(公式(1)至公式(27))设计并开发地质界面三维形态分析软件TriangelCal，包括五个功能模块：①模型数据读取模块；②地质界面距离场计算模块；③地质界面一般几何形态参数提取模块；④地质界面趋势-起伏分析计算模块；⑤结果保存模块。以下各步骤即步骤(3)至步骤(5)都是采用该软件进行计算的。

(3)地质界面的距离场计算：以单元到不整合面的最小距离作为距离场值，即不整合面距离场因素指标dU。公式(14)和公式(15)求得立体单元中心点到地质界面的垂足坐标，继而判断垂足是否在三角形内，然后根据公式(16)和公式(17)分别求得不同情况下的距离。为区分不整合面上下单元的距离场值的差异，将位于不整合面分布范围之上的单元场值置为正(正距离场)，位于不整合面分布范围之下的单元场值为负(负距离场)。

(4)地质界面的一般几何形态参数(坡度、夹角等)提取：利用公式(5)和公式(6)得到三角形平面方程参数，代入到公式(2)得到不整合面坡度因素gU；求不整合面夹角因素aU_S时，先对不整合面以及其他地质界面栅格化，并求得交集，然后利用公式(7)和公式(8)求得地质界面相交处的平均夹角，由于不整合面与多个地层界面有相交，定义单元的不整合面夹角指标值取值为单元到不整合面的最近距离处对应的最近的夹角。

(5)分级提取出地质界面的形态趋势以及形态起伏：选择100m和200m的插值搜索范围半径，按照距离平方反比法即公式(20)，对不整合面原始TIN模型进行一级形态滤波和二级形态滤波，分别得到不整合面的一级趋势与一级起伏waU、二级趋势与二级起伏wbU。

(6)地质界面三维形态分析结果的应用：将存储在栅格模型Access数据库中的三维形态分析结果，导入Datamine软件，形成控矿地质因素的块体模型，利用Datamine的块体可视化功能，完成控矿地质因素的可视化表达；对Access数据库中的已知矿化指标和控矿地质因素指标进行回归分析，建立矿体预测模型(建立反映控矿指标到矿化分布的非线性映射关系的隐伏矿体立体定量预测模型，包括用于预测单元中品位、金属量的多元回归模型和用于预测含矿性的逻辑斯蒂回归模型，将未知单元的控矿地质因素指标代入模型，即可达到对未知单元中矿的品位、金属量、含矿性指标进行预测的目的)，预测得到铅的远景金属量为123917.37吨，锌的远景金属量为520334.04吨。预测结果可以通过实施找矿工程进行验证。

Claims (2)

1.一种基于TIN的地质界面三维形态分析方法，其特征在于，包括以下步骤：在使用三维TIN模型来模拟地质界面的基础上，使用空间几何和计算机图形学原理计算地质界面的一般几何形态参数和距离场；利用空间插值和趋势剩余分析技术分级提取出地质界面的形态趋势以及形态起伏。

2.根据权利要求1所述的基于TIN的地质界面三维形态分析方法，其特征在于，所述的一般几何形态参数包括坡度和夹角；

(1)坡度的计算方法为：在三维TIN模型中，任一个三角形ABC的三个顶点的坐标为A(x₁，y₁，z₁)，B(x₂，y₂，z₂)，C(x₃，y₃，z₃)，三角形ABC的坡度为：

$S=\arccos \frac{1}{\sqrt{a_1^2+a_2^2+1}}$

其中， $a_1=\frac{(y_2-y_1)*(z_3-z_1)-(y_3-y_1)*(z_2-z_1)}{(x_3-x_1)*(y_2-y_1)-(x_2-x_1)*(y_3-y_1)},$ $a_2=\frac{(x_3-x_1)*(z_2-z_1)-(x_2-x_1)*(z_3-z_1)}{(x_3-x_1)*(y_2-y_1)-(x_2-x_1)*(y_3-y_1)};$

(2)夹角的计算方法为：

所述的夹角指两个地层界面交线上的某点的夹角，即经过该点的两个三角面片之间的夹角，方法如下：

首先将地质界面栅格化以形成多个块体单元即单元，即采用三维规则网格将研究区地质空间划分为许多小立方体，对两个地质界面TIN模型中穿过相同块体单元的三角面片两两求夹角，最后将求得的夹角的平均值作为地质界面之间的夹角α^*，夹角α^*的求法如下：

设两个三角形的法向量分为V₁(x_v1，y_v1，z_v1)、V₂(x_v2，y_v2，z_v2)，两法向量的夹角为α，即两个三角面片之间的夹角，计算的公式为：

$\mathrm{\α}=\arccos \frac{x_{v1}*x_{v2}+y_{v1}*y_{v2}+z_{v1}*z_{v2}}{\sqrt{{x_{v1}}^2+{y_{v1}}^2+{z_{v1}}^2}\sqrt{{x_{v2}}^2+{y_{v2}}^2+{z_{v2}}^2}}$

则地质界面之间的夹角α^*的表达式为：

其中n等于两个地质界面在同一块体单元中的三角面片的夹角个数。

根据权利要求1所述的基于TIN的地质界面三维形态分析方法，其特征在于，距离场的计算方法为：

控矿地质因素场与空间中某点到相关联的地质界面的距离有关，即控矿地质因素场是空间中某点到地质界面距离的空间分布函数，

求出预测空间即前述的研究区地质空间中每个单元对应的距离值，建立距离场模型，所述的距离值指某一个单元的中心点到TIN模型中各个三角面片的距离值的最小值。

分级提取出地质界面的形态趋势以及形态起伏包括：

空间趋势形态反映的是空间物体在空间区域上的主体特征，它忽略了局部形态起伏以揭示空间区域上的主体特征。在矿床地质空间中，相关地质体的主体形态将对周围成矿起到一定因素的影响，起伏则是地质异常作用的侧面反映。因此，对地质体形态的分析将是控矿因素定量提取的重要步骤。

所述的地质界面的形态趋势是指对地质界面TIN模型中的三角面片的顶点的Z值按照距离平方反比法重新赋值得到的新顶点建立的新的TIN模型，新的TIN模型称为地质界面的趋势面，地质界面的趋势面的形成过程如下：

①确定搜索距离d，只有与估测点距离小于d的样品点才能参与计算，其中距离是指两点在XY平面上投影点间的距离，而不是两点的绝对距离；估测点指原来的地质界面TIN模型中的任意一个顶点在XY平面上的投影点；

②利用距离平方反比法重新计算样品点的所有z值，公式如下：

$z=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{z_i}{d_i^2}}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{d_i^2}}$

其中，z_i表示该估测点周围、与估测点距离小于d的样品点的垂直坐标分量；d_i表示该估测点与样品点的距离；此处距离都指两点在XY平面上投影点间的距离；

③利用计算的z值，建立新的地质界面TIN模型，形成地质界面的趋势面；

按照上述方法，采用不同的搜索距离d，即可分级提取得到地质界面的形态趋势；

所述的形态起伏是指地质界面相对于形态趋势面的起伏幅度，即地质界面中的z值减去形态趋势面的z值，即可得到形态起伏。

Images