地质体的三维形态分析方法与应用

毛先成¹，唐艳华¹，邓浩²

(1. 中南大学 地球科学与信息物理学院，教育部有色金属成矿预测实验室，湖南 长沙，410083；

2. 北京师范大学 地理学与遥感科学学院，北京，100875)

摘要：针对隐伏矿体立体定量预测中复杂形态地质体的空间分析问题，利用数学形态学和欧氏距离变换理论，提出基于三维栅格模型的地质体三维形态定量分析的方法。利用开闭滤波或闭开滤波的数学形态学变换，构建地质体表面形态趋势的提取算法；将形态滤波、欧氏距离变换、全局集合运算结合起来，建立地质体形态起伏定量分析及分级提取模型；采用基于欧氏距离场的计算模型，提出地质体表面与趋势面夹角的提取方法。在建立地质体三维栅格模型的基础上，将地质体三维形态定量分析方法应用于安徽铜陵凤凰山矿田的新屋里岩体分析，实现地质控矿因素指标的定量提取。研究结果表明：采用该方法对地质体进行形态分析，可计算出地质体的各种形态参数，并提取出地质控矿因素量化指标，能有效地应用于隐伏矿体立体定量预测。

关键词：地质体；三维形态分析；数学形态学；欧氏距离变换

中图分类号：P612；P628          文献标志码：A         文章编号：1672-7207(2012)02-0588-08

Three-dimensional morphological analysis method for geologic bodies and its application

MAO Xian-cheng¹, TANG Yan-hua¹, DENG Hao²

(1. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education,School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;

2. School of Geography, Beijing Normal University, Beijing 100875, China)

Abstract: For the spatial analysis of complex geological bodies in stereoscopic quantitative prediction of concealed ore bodies, a new method of three-dimensional morphological analysis for geologic bodies based on 3D raster model was put forward by combining mathematical morphology with Euclidean distance transformation theory. The algorithm for extracting the surface shape trend of geologic body was proposed by using open-close filter and close-open filter of mathematical morphology. By combination of morphological filtering, Euclidean distance transforming, and global set operation, the models for the quantitative analysis and hierarchical extraction of the shape undulance were established. By means of the calculation model of Euclidean distance field, an approach for extracting the angle parameter between the geologic body surface and its trend surface was acquired. After building the 3D raster models of the geologic bodies, a case study of the three-dimensional morphological analysis method was applied in analyzing quantitatively the Xinwuli magmatic body in Fenghuangshan ore field in Tongling, Anhui Province, and the quantitative extraction of the quantitative indexes of geological ore-controlling factors was achieved. The results show that the morphological analysis method for geologic bodies can be used to calculate various morphological parameters of geological bodies, to extract quantitative indexes of geologic ore-controlling factors successfully, and it plays an important role in stereoscopic quantitative predication of concealed ore bodies.

Key words: geologic body; 3D morphological analysis; mathematical morphology; Euclidean distance transform

随着我国经济的连续高速发展，矿产资源的消耗量和需求量呈现跳跃式增加，对矿产资源的需求日益迫切。而近年来，一大批主力矿山由于各种原因导致资源迅速减少或资源枯竭面临关闭，危机矿山问题成为矿产资源可持续发展的最突出问题。我国大多数老矿山及危机矿山仍具有巨大的资源潜力，围绕矿山基地开展找矿勘查是解决危机矿山接替资源问题的有效途径^[1]，因此，迫切需要采用大比例尺下的隐伏矿体立体定位定量预测技术来扩大和发现新的资源量，即在其深部及边部寻找新的接替资源^[2]。这就需要对矿床中的地质对象进行真三维模拟，通过空间分析提取地质控矿指标，进而发掘有利成矿信息并对未知区域进行定量预测^[3-4]。矿化空间分布常常与地质体形态存在着密切的空间关联关系，采用空间分析方法对地质体三维形态进行定量分析对于成矿信息定量提取具有重要意义。目前，大多三维空间分析方法多来自计算机图形学与图像处理领域，这些方法注重的是几何精度和执行效率，已有的研究包括三维距离量算^[5-6]、三维叠加分析^[7]与三维领域查找^[8]等。数学形态学是图像处理理论的一个重要方面，它被广泛地应用到图像处理的很多领域中。近年来，国内外学者对利用数学形态学处理三维图像和模型进行了许多研究，包括医学图像的交互分割^[9-11]、器官疾病特征检测^[12-13]、空间形态滤波^[14-17]、风暴体的自动识别^[18]、雷达云中的建筑提取^[19]等。这些研究多是基于分层扫描图像或准三维图像的分析处理，而不是针对真三维结构。目前，数学形态学在医学领域应用广泛，而尚未在成矿定量预测方面得到较好的应用。在此，本文作者根据地质体形态分析的实际需要，将地质体三维建模与数学形态学、欧氏距离变换理论相结合，提出基于三维地质体栅格模型的地质体三维形态分析方法，包括地质体表面形态趋势与形态起伏分析、地质体形态参数提取等。该方法可适应于超覆、弯曲(褶皱)、错断等形态复杂的地质体，并满足隐伏矿体立体定量预测的真三维空间要求。此外，还利用该方法对安徽铜陵凤凰山矿田的新屋里岩体形态进行研究分析。

1  地质体三维形态分析方法

1.1  数学形态学

数学形态学是以严格的数学理论和集合理论为基础，用具有一定形态的结构元素去量度和提取图像中的对应形状以达到对图像分析和识别的目的。数学形态学的基本运算主要包括膨胀、腐蚀、开、闭等^[20]，其中腐蚀与膨胀、开与闭互为对偶运算。

腐蚀(erosion)运算“Θ”对集合元素采用向量减法，将2个集合合并；膨胀(dilation)运算“”是腐蚀运算的对偶运算，可以通过对补集的腐蚀来定义。开(opening)运算“”就是先对图像进行腐蚀然后膨胀其结果；闭(closing)运算“”是先对图像进行膨胀然后腐蚀其结果。集合A被B腐蚀，表示为；A被B膨胀表示为；和分别表示开启变换和闭合变换。

设A表示一幅二值图像，B为结构元素：

           (1)

             (2)

             (3)

            (4)

式中：为A的补集；为B关于坐标原点旋转180°后的集合。

1.2  欧氏距离变换

在基于栅格数据结构的地质体体素模型的空间分析中，需要大量的距离量算，因此，需要用一种快速量算距离，在运行栅格空间中生成距离场的解决方案，而距离变换(Distance transform)^[21]正是这样一种方法。由于离散化表示的矿床地质空间是一种地质现象复杂的欧氏空间，为了在控矿地质作用空间分析中与实际欧氏距离相吻合，保证基于距离量算的控矿因素定量提取精度，采用欧氏距离作为距离测度，即需要对整体或部分的离散化矿床地质空间进行三维欧氏距离变换。本文采用8个模板进行三维欧氏距离变换的方  法^[22]，对体素模型的栅格空间进行扫描，生成欧氏距离场。

1.3  地质体表面形态趋势与起伏分析

1.3.1  地质体形态趋势提取

在矿床地质空间中，相关地质体的总体形态对周围成矿起到一定程度的影响和控制，所以，地质体总体形态的分析将是控矿因素定量提取的重要步骤。一般分析方法选用趋势面数学方程，而这种趋势面分析的方法在地质体形态复杂、存在有超覆现象、无法依靠地质体曲面重建来构造地质体表面形态时，无法满足要求。在此，本文在地质体体素栅格模型的基础上，提出一种基于数学形态学的通用趋势形态提取方法。

定义1个球形结构元素，因为球形在三维空间中是各向同性的，没有任何一个方向具有特殊性。用这种结构元素对对象进行开运算或闭运算时，若球在整个三维空间中进行滚动，则必将包括球在地质体边界(可视为一个自我封闭曲面)内外的滚动。球在地质体边界内部滚动时为开运算。一次开运算结果的边界将是球状内边界上滚动的轨迹组合成的一个薄壳的外边界，因为球在内边界上滚动将会削平边界上的凸峰(见图1(b)和(c))，因此，该闭合曲面较原地质边界滤除了球无法滚进的凸峰部分。同理，闭运算时球在地质体外边界上滚动，此时，球的滚动将会填补边界上的凹谷(见图1(d)和(e))，从而得到一个被滤除了凹谷的轮廓。因此，用这2种运算的组合对地质体进行运算操作，可分别削平和填补地质体边界上的凸峰和凹谷，得到平滑的趋势形态。

图1  球形结构元素开运算与闭运算效果

Fig.1  Open and close operation of spherical structure element

设A为地质体对象，B_ball为球形结构元素，将这2种运算进行组合定义成形态滤波变换。

(1) 开闭滤波：

          (5)

(2) 闭开滤波：

          (6)

通过这2种变换均能得到地质体的趋势形态。变换后的地质体轮廓光滑度取决于球形结构元素的半径。半径越大，能削平和填补的凸峰与凹谷部分越多，轮廓(趋势形态)越光滑；半径越小，能削平和填补的凸峰与凹谷部分越少，轮廓(趋势形态)将相对趋近地质体变换前轮廓(形态)。

1.3.2  地质体表面形态起伏提取

对于复杂地质体，可能无法通过一般方法^[23-24]来对地质界面进行重建，进而描述界面的形态起伏。借助欧氏距离变换技术，本文作者提出一种新的形态起伏分析方法，将边界上所有体元的凹凸属性及其起伏程度量化。这种方法只需进行1次形态滤波、1次欧氏距离变换、2次全局的集合运算就可以得到边界上所有体元的凹凸属性及其起伏程度。具体步骤如下。

Step 1：初始化，建立要进行形态起伏分析的地质体的三维二值图像，定义一定半径r的球形结构元素B_ball(r)。r决定可滤除波形的幅度，可滤除波形幅度则决定可得到的起伏程度。

Step 2：以B_ball(r)为结构元素，采用开闭或闭开滤波算子对地质体A进行形态滤波，得到集合。

Step 3：获得外凸部分集合和内凹部分集合，有：

   (7)

Step 4：通过欧氏距离变换，建立集合对象外的欧氏距离场E_outer(i，j，k)和内的欧氏距离场E_inner(i,j,k)。

Step 5：对于地质体边界上的每个体元v，有。

若，则v处于外凸处，且有，转Step 6；

若，则v处于内凹处，有，转Step 7；

否则，说明v处于平坦处。

Step 6：求在体元v坐标处的距离场E_outer，该值说明了v所在的局部外凸程度，转Step 5。

Step 7：求在体元v坐标处的距离场E_inner；该值说明了v所在的局部内凹程度，转Step 5。

1.3.3  形态起伏特征的分级提取

同样，对以上步骤进行改进便可对形态起伏特征进行分级提取。因为球形结构元素半径r决定可滤除的波形幅度，所以，通过改变r可提取出不同程度的起伏。下面给出分级提取算法的具体步骤。

Step 1：初始化，建立要进行形态起伏分析的地质体的三维二值图像，定义需要提前的级别数、各个级别的r。设初始地质体为A₀，建立表达三维实体(地质体或趋势形态体)外部的欧氏距离场E_outer(i，j，k)和内部的欧氏距离场E_inner(i，j，k)。

Step 2：取第i级r构造球形结构元素，对上一级滤波结果A_i-1进行形态滤波，得到。

Step 3：获得第i级下外凸和内凹部分的集合，有：

   (8)

Step 4：对于地质体边界上的每个体元v_0,k，在距离场E_outer(i，j，k)或E_inner(i，j，k)中获取其在A_i-1上的最近体元v_i-1,k；存储每个v_0,k?v_i-1,k的映射。

Step 5：进行欧氏距离变换，更新E_outer(i，j，k)和E_inner(i，j，k)：在E_outer(i，j，k)中建立A_i外的欧氏距离场，在E_inner(i，j，k)中建立A_i内的欧氏距离场。

Step 6：对每个v_i-1,k，进行如下判断：

(1) 若，则v_i-1,k处于外凸处，且有，转Step 7。

(2) 若，则v_i-1,k处于内凹处，且有，转Step 8。

(3) 否则，说明v_i-1,k处于平坦处。

Step 7：求在体元v_i-1,k坐标处的距离场E_outer，该值反映了v_i-1,k所在的局部外凸程度。通过v_0,k?v_i-1,k映射得到地质体边界上体元v_0,k，该值也反映了v_0,k在第i级的起伏程度，将该值赋为正值并存储，转回Step 6。

Step 8：求在体元v_i-1,k坐标处的距离场E_inner，该值反映了v_i-1,k所在的局部内凹程度。通过v_0,k?v_i-1,k映射得到地质体边界上体元v_0,k，该值也反映了v_0,k在第i级起伏程度。将该值赋为负值并存储，转回Step 6。

Step 9：若下一级r未超过最大阈值，则令i=i+1，转Step 2，否则，A_i为最终趋势部分(主体形态)，结束。

通过以上步骤可得到地质体边界上所有体元的各级别起伏程度。由这种方法提取出来的起伏程度随着级别的增高而增大，而非级别越高，起伏越小。

1.4  地质体表面与趋势面夹角提取

一些成矿因素受到地质体接触面(如岩浆岩体与围岩的接触带)夹角的控制。为了定量表达这种控矿因素，需要对接触面走向与其趋势面走向方向的夹角进行提取，也需要对接触面走向与区域应力场方向夹角进行提取。地质体以体素栅格模型离散表达，这种夹角提取问题具体为求2个空间实体上对应两体元切平面之间的角度问题。如图2所示，地质体表面(实线) 上有1点P，它在地质体的趋势形态表面(虚线)上的最近点为P′，求过这2点的切平面间夹角θ。针对这个夹角提取问题，本文作者提出基于欧氏距离场的夹角提取方法。

图2  接触面间夹角示意图

Fig.2  Angle between interfaces

对于切平面，可由体元v中心坐标(x₀，y₀，z₀)过切点P的法向量N(A，B，C)来定义。设曲面由以下方程给出：

                       (9)

≤  是曲面上的1点，并设函数的偏导数在该点连续且不同时为0。在曲面上，通过点M任意引1条曲线Γ，假定曲线Γ的参数方程为：

   (10)

令t=t₀对应于点，且，和不全为0。该曲线的切线方程为：

         (11)

可以证明：在曲面上通过点P且在点P处具有切线的任何曲线，它们在点P处的切线在同一个平面上。曲面在点M的切平面方程为：

                  (12)

其法向量为：

(13)

对于F(x，y，z)，本文通过欧氏距离场的方法来近似模拟。对地质体内外体元进行欧氏距离变换，可得到整个空间的欧氏距离场。设体元v中心坐标(x₀，y₀，z₀)满足：当体元在地质体内部、距离小于0时，F(x₀，y₀，z₀)＜0；体元在地质体表面、距离为0时，此时，F(x₀，y₀，z₀)=0；体元在地质体外、距离大于0时，F(x₀，y₀，z₀)＞0。而对于F_x，F_y和F_z，可通过梯度算子来模拟，本文采用Sobel算子^[25]来近似计算，设算子模板大小为3×3×3。求出x，y和z 3个方向梯度G_x，G_y和G_z后，进而求得点P的法向量(G_x，G_y，G_z)。

设原始岩体A接触面S(A)上有体元v，其在趋势形态体A′表面S(A′)上对应体元为v′，过v和v′的法向量分别为n和n′，则求原始接触面上某体元与其趋势形态体上对应体元切平面间夹角α的具体步骤如下。

Step 1：求岩体接触面，用数学形态学方法，有(其中，B为方形结构元素)。利用前面的方法求地质体趋势形态。同样，按数学形态学方法，有：。

Step 2：求S(A)外和内的欧氏距离场，分别为E_outer(i，j，k)和E_inner(i，j，k)；求S(A′)外和内的欧氏距离场，分别为E′_outer(i，j，k)和E′_inner(i，j，k)，以E_outer(i，j，k)和E_inner(i，j，k)所表示的距离模拟F(x，y，z)；E_outer(i，j，k)的距离为正值，E_inner(i，j，k)中的距离为负值，S(A)上距离为0；同理，以E′_outer(i，j，k)和E′_inner(i，j，k)距离模拟F_趋势(x，y，z)(地质体趋势部分的曲面函数，下同)。

Step 3： 对S(A)中需要求夹角的体元v，在 E′_outer(i，j，k)或E′_inner(i，j，k)中找出其最近在S(A′)上的最近体元作为对应体元v′。

Step 4：用模板对F(x，y，z)进行卷积运算，求得在v处x，y和z 3个方向的梯度G_x，G_y和G_z，n=(G_x,G_y,G_z)；同样，对F_趋势(x,y,z)进行卷积运算，求得n′。

Step 5：求两切平面间夹角α。

               (14)

Step 6：判断S(A)中有没有尚未求的夹角体元，若有，则回到Step 3。

2  地质体形态分析应用实例

以安徽铜陵凤凰山矿田的新屋里岩体为例。由于安徽铜陵凤凰山矿田新屋里岩体存在有广泛的超覆现象，因此，采用本方法进行形态起伏分析。地质体三维形态分析和控矿地质因素场建模方法都以三维栅格模型为数据模型，将已建立的各种地质体块体模型(栅格模型)作为其原始输入数据。根据凤凰山矿田矿体定位预测概念模型，控制矿体分布的主要地质体与地质因素包括新屋里岩体及其形态起伏、接触带、地层及褶皱构造、断层等，需要描述这些抽象控矿地质因素三维分布的场变量。针对凤凰山矿田中的各类地质体的三维块体模型，采用前述的地质体三维形态分析技术和控矿地质因素场建模技术，研究开发了三维形态学处理算法与程序，对新屋里岩体等地质体进行三维形态分析，如岩体形态趋势提取、岩体表面形态起伏提取、形态起伏特征分级提取、地质体表面和趋势面夹角提取等，实现各种控矿地质因素场的定量模拟。

新屋里岩体原始模型和趋势形态模型见图3。由图3可以看出：当选择的结构元素半径足够大时，获得的形态趋势面变得平滑，由此可以分析岩体的形态轮廓。

根据前述地质体表面形态起伏提取方法对新屋里岩体进行形态起伏分析。图4所示为按这种方法用半径为100 m和1 000 m的球形结构元素对新屋里岩体进行二级形态起伏分析所获得的1级形态起伏提取图与2级形态起伏提取图。岩体形态因素w_r1G和w_r2G分别代表岩体表面一级起伏和二级起伏的起伏程度，用单元距离最近的岩体体元的岩体趋势部分的欧氏距离E_outer或E_inner来度量，计算结果如图5所示。从图5可见：若岩体体元属于岩体相应级别起伏的外凸部分，则w_r1G和w_r2G均大于0；若岩体体元属于岩体相应级别起伏的内凹部分，则w_r1G和w_r2G均小于0。

控制凤凰山矿田铜矿体分布的成矿断裂具有明显的偏北西专属性，这种专属性主要与区域构造应力场有关。为了表达区域挤压远应力场对矿化分布的控制作用，构造挤压应力方向面与新屋里岩体原始接触面的夹角a_IP作为该因素的指标。计算挤压应力方向面与原始接触面的夹角；对于地质空间中的任意单元，该结果取其最近接触面单元的夹角，如图6所示。

图3  新屋里岩体原始模型和趋势形态模型

Fig.3  Original model of rock mass and its model of trend

图4  新屋里岩体形态起伏分级提取外凸、内凹部分

Fig.4  Hierarchically extracting concave and convex parts of rock mass’s undulate shape

图5  岩体形态因素(w_r1G和w_r2G)分布栅格模型(-200 ~ -700 m标高范围)

Fig.5  Raster model (from -200 m to -700 m elevation) of rock mass factors (w_r1G and w_r2G)

图6  区域挤压远应力场因素(a_IP)分布栅格模型(-200 ~ -700 m标高范围)

Fig.6  Raster model (from -200 m to -700 m elevation)of regional extruding far crustal stress field factor a_IP

3  结论

(1) 结合数学形态学和欧氏距离变换理论，提出了基于三维栅格模型的地质体三维形态分析数学模型和算法，包括地质体表面形态趋势分析、分级形态起伏提取与定量分析、地质体表面与趋势面夹角提取等，为隐伏矿体立体定量预测中复杂地质体形态分析提供了一种新的方法与实现途径。

(2) 将提出的方法应用于安徽铜陵凤凰山矿田新屋里岩体等复杂地质体的形态分析，解决了趋势分析等传统方法无法处理复杂地质体的超覆现象等问题。趋势形态定量提取避免了复杂地质体表面无法进行趋势面重建的矛盾，只需要对地质体进行多次开闭滤波运算即可实现。地质体表面形态起伏定量提取，是通过1次形态滤波、1次欧氏距离变换、2次全局的集合运算实现的，可以得到边界上所有体元的凹凸属性及其起伏程度；当改变球形结构元素的半径时，即可进行形态起伏的分级提取。通过地质体表面与趋势面的欧氏距离场的计算，可获取到地质体表面及其趋势面的夹角、区域主应力方向面与岩体原始接触面的夹角。

(3) 采用地质体三维形态分析方法，建立可视化地质体形态模型，提取出地质体形态相关的几何参数，构建控矿地质因素场，定量地描述和分析了矿化分布的控制作用，实现了控矿地质因素分析的定量化，对危机矿山深边部隐伏矿体立体定量预测及找矿具有重要意义。

参考文献：

[1] LIU Liang-ming, PENG Sheng-lin. Key strategies for predictive exploration in mature environment: Model innovation, exploration technology optimization and information integration[J]. Journal of Central South University of Technology, 2005, 12(2): 186-191.

[2] 毛先成, 戴塔根, 吴湘滨, 等. 危机矿山深边部隐伏矿体立体定量预测研究: 以广西大厂锡多金属矿床为例[J]. 中国地质, 2009, 36(2): 424-435.
MAO Xian-cheng, DAI Ta-gen, WU Xiang-bin, et al. The stereoscopic quantitative prediction of concealed ore bodies in the deep and marginal parts of crisis mines: A case study of the Dachang tin polymetallic ore deposit in Guangxi[J]. Geology in China, 2009, 36(2): 424-435.

[3] MAO Xian-cheng, ZOU Yan-hong, LU Xiao-qin, et al. Quantitative analysis of geological ore-controlling factors and stereoscopic quantitative prediction of concealed ore bodies[J]. Journal of Central South University of Technology, 2009, 16(6): 987-993.

[4] 毛先成, 邹艳红, 陈进, 等. 危机矿山深部、边部隐伏矿体的三维可视化预测: 以安徽铜陵凤凰山矿田为例[J]. 地质通报, 2010, 29(2/3): 10-22.
MAO Xian-cheng, ZOU Yan-hong, CHEN Jin, et al. Three-dimensional visual prediction of concealed ore bodies in the deep and marginal parts of crisis mines: A case study of the Fenghuangshan ore field in Tongling, Anhui, China[J]. Geological Bulletin of China, 2010, 29(2/3): 10-22.

[5] Sud A, Govindaraju N, Gayle R, et al. Interactive 3D distance field computation using linear factorization[C]//Proceedings ACM Symposium on Interactive 3D Graphics and Games. Chapel Hill, USA: Association for Computing Machinery, 2006: 117-124.

[6] Jones M W, Baerentzen J A, Sramek M. 3D distance fields: A survey of techniques and applications[J]. Transactions on Visualization and Computer Graphics, 2006, 12(4): 581-599.

[7] Sakuragi F, Nakamoto M, Sasama T, et al. System simulator for structural description and error analysis of multimodal 3D data integration systems[J]. Electronics and Communications in Japan, 2007, 90(8): 325-338.

[8] Zuliansyah M, Supangkat S H, Priyana Y, et al. 3D topological relations for 3D spatial analysis[C]//2008 IEEE International Conference on Cybernetics and Intelligent Systems, CIS 2008. Piscataway, NJ, USA: Institute of Electrical and Electronics Engineers, 2008: 585-590.

[9] Fetita C, Brillet P Y, Prêteux F J. Morpho-geometrical approach for 3D segmentation of pulmonary vascular tree in multi-slice CT[C]//Pluim J P W, Dawant B M. Medical Imaging 2009: Image Processing. Lake Buena Vista, FL, USA: Proceeding of SPIE, 2009: F1-F12.

[10] Felix J H S, Cortez P C, Holanda M A, et al. Lung and chest wall structures segmentation in CT images[C]//Proceedings of VIPIMAGE 2007: 1st ECCOMAS Thematic Conference on Computational Vision and Medical Image Processing. Netherlands: Taylor and Francis/Balkema, 2008: 291-294.

[11] Naegel B. Using mathematical morphology for the anatomical labeling of vertebrae from 3D CT-scan images[J]. Computerized Medical Imaging and Graphics, 2007, 31(3): 141-156.

[12] Chien K C C, Fetita C, Brillet P Y, et al. Detection and classification of interstitial lung diseases and emphysema using a joint morphological-fuzzy approach[C]//Karssemeijer N, Giger M L. Medical Imaging 2009: Computer-Aided Diagnosis. Lake Buena Vista, FL, USA: Proc. SPIE, 2009: 0.1-0.8.

[13] Jiang Y, Meng J, Jaffer N. A novel segmentation and navigation method for polyps detection using mathematical morphology and active contour models[C]//Proceedings of the 6th IEEE International Conference on Cognitive Informatics, ICCI 2007. Piscataway, NJ, USA: Institute of Electrical and Electronics Engineers, 2007: 357-363.

[14] Kiwanuka F N, Ouzounis G K, Wilkinson M H F. Surface-area-based attribute filtering in 3D[C]//Mathematical Morphology and Its Application to Signal and Image Processing. Netherlands: Springer Berlin / Heidelberg, 2009: 70-81.

[15] Tankyevych O, Talbot H, Dokladál P, et al. Spatially-variant morpho-hessian filter: Efficient implementation and application[C]//Mathematical Morphology and Its Application to Signal and Image Processing. Netherlands: Springer Berlin/ Heidelberg, 2009: 137-148.

[16] 李勇, 吴华意. 基于形态学梯度的机载激光扫描数据滤波方法[J]. 遥感学报, 2008, 12(4): 633-639.
LI Yong, WU Hua-yi. Filtering airborne LIDAR data based on morphological gradient[J]. Journal of Remote Sensing, 2008, 12(4): 633-639.

[17] 刘彩虹. 基于立体结构元素的三维形态滤波算法的研究[J]. 长春工程学院学报: 自然科学版, 2009, 10(2): 85-88.
LIU Cai-hong. Study of three-dimensional morphological filtration algorithm on structuring element[J]. Journal of Changchun Institute of Technology: Natural Sciences Edition, 2009, 10(2): 85-88.

[18] 韩雷, 郑永光, 王洪庆, 等. 基于数学形态学的三维风暴体自动识别方法研究[J]. 气象学报, 2007, 65(5): 806-814.
HAN Lei, ZHENG Yong-guang, WANG Hong-qing. 3D storm automatic identification based on mathematical morphology[J]. Acta Meteorologica Sinica, 2007, 65(5): 806-814.

[19] TAO Jin-hua, SU Lin, LI Shu-kai. Method of extracting building model from lidar point cloud[J]. Infrared and Laser Engineering, 2009, 38(2): 340-345.

[20] 唐常青, 吕宏伯, 黄铮, 等. 数学形态学方法及其应用[M]. 北京: 科学出版社, 1990: 10-22.
TANG Chang-qing, L? Hong-bo, HUANG Zheng, et al. Mathematical morphological method and its application[M]. Beijing: Science Press, 1990: 10-22.

[21] Fabbri R, Costa L D F, Torelli J T, et al. 2D Euclidean distance transform algorithms: A comparative survey[J]. ACM Computing Surveys, 2008, 4(1): 2: 1-44.

[22] 蔺宏伟, 王国瑾. 三维带符号的欧氏距离变换及其应用[J]. 计算机学报, 2003, 26(12): 1645-1651.
LIN Hong-wei, WANG Guo-jin. Three dimensional signed euclidean distance transform and its applications[J]. Chinese Journal of Computers, 2003, 26(12): 1645-1651.

[23] 郭仁忠. 空间分析 [M]. 2版. 北京: 高等教育出版社, 2001: 77-109.
GUO Ren-zhong. Spatial analysis[M]. 2nd ed. Beijing: Higher Education Press, 2001: 77-109.

[24] 於崇文. 数学地质的方法与应用[M]. 北京: 冶金工业出版社, 1980: 246-289.
YU Chong-wen. Mathematical geology method and its application[M]. Beijing: Metallurgical Industry Press, 1980: 246-289.

[25] Gonzalez R G, Woods R E. 数字图像处理[M]. 2版. 阮秋琦, 阮宇智, 等译. 北京: 电子工业出版社, 2003: 463-474.
Gonzalez R G, Woods R E. The digital image processing[M]. 2nd ed. YUAN Qiu-qi, YUAN Yu-zhi, et al, transl. Beijing: China Electronic Industry Press, 2003: 463-474.

(编辑 陈灿华)

收稿日期：2011-01-10；修回日期：2011-03-25

基金项目：国家自然科学基金资助项目(41172297)；国家“十一五”科技支持计划项目(2006BAB01B07)

通信作者：毛先成(1963-)，男，湖南澧县人，教授，博士生导师，从事隐伏矿床预测、地学信息技术研究；电话：0731-88877571；E-mail：xcmao@126.com