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

毛先成 ，赵莹 ，唐艳华 ，彭正林 ，陈进 ，邓浩

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

界面是指两相间接触的交界部分，两物态间的转换面[1]。地质体是地质空间中不连续的空间实体，包括沉积成因的岩层、侵入成因的岩浆岩体以及受力变形的构造等[2]。宏观的地质界面可以理解为地质体的分界面或间断面，或者称为面状地质体[2-3]。在地质学上，界面两侧的岩石具有不同的岩性类型、不同的矿物组成、不同的化学成分、不同的结构构造。地质界面是地球内部能量集中、汇集、传递、转化和释放的地带，是地球物质变化、活化、迁移和沉淀富集的部位，也是成矿作用发生的场所。自然界中，地质界面成矿现象普遍存在，矿产资源往往在各种形式的地质界面处超常富集。运用地质界面控矿原理进行找矿预测，关键是准确识别并定位地质界面，研究其性质并理解其控矿作用，并融合其他成矿理论进行综合分析，同时注意地质界面控矿在地质异常解释中的运用[1]。对地质体数学特征的研究有助于查明地质体成因，追溯地质体形成过程和机理，并有助于发现一些隐蔽规律性，从而更好地解决认识、预测和评价等理论和实际问题[4]。所以，采用空间分析方法对地质体和地质界面三维几何形态等数学特征进行提取和定量分析，对于发掘有利成矿信息具有重要意义。现有的研究多关注于地质体和地质界面模型的建立和表达[5-6]等，但是，在利用三维空间形态分析提取地质控矿因素方面，见到的研究文献还比较少。针对控矿地质因素提取的三维形态分析问题，本文作者曾提出采用基于栅格模型的地质体三维形态分析方法进行控矿地质因素的分析和提取[7]，该方法是地质三维建模与数学形态学[8-10]相结合形成的一种基于栅格的三维空间分析方法。基于栅格模型的地质体三维形态分析方法使用体素模型(三维栅格模型)表达地质对象空间数据，具有位置隐含、结构简单、操作方便的优点，可描述地质现象(地质场)连续变化的特性，特别适应于超覆、弯曲(褶皱)、错断等复杂形态的地质体。但是，研究结果分析表明，该方法存在如下缺点：栅格数据模型空间数据量大、空间精度相对较低；在进行边界表示时不能准确表达地质体的边界[11]，且生成速度、存储空间和表达精度难以平衡；数学形态学球形结构元素半径的选取在量化方面需要改进；不适合对薄片状地质体或存在狭窄部分的地质体进行趋势形态分析，当结构元素半径较大时，会出现断裂的错误。因此，还需提出一种更合适地质界面表达与三维形态分析的方法。基于TIN模型的表达方法可有效地模拟地质界面[12]。不规则三角网(Triangulated Irregular Network, TIN)模型可以根据表面的复杂程度变化三角形面片的大小和数量，消除可视化时的数据冗余，并且具有三维面表示模型的易于采集与构建、表达精度高、便于显示和更新的优点，适于进行地质体轮廓的几何建模，建模理论成熟，方法易于实现。基于上述分析，本文作者根据地质界面形态分析的实际需要，将利用TIN模型表达和分析地质界面精度高等方面的优势，与三维地质建模、三维空间分析等方法结合起来，提出了一种精度更高的地质界面形态分析方法—基于TIN的地质界面三维形态分析方法，可以提取地质界面的几何形态参数、分析地质界面的多级形态趋势与起伏等，并建立控矿地质因素场模型，实现控矿地质因素的定量表达与可视化。

1 地质界面三维形态分析方法

1.1 TIN模型与地质界面三维形态分析

TIN模型将空间无重复点的散乱点按某种规则(如Delaunay规则)进行三角剖分而形成连续却不重叠的不规则三角面片网，以此来模拟三维物体的表面。它具有固定的结构和简单的处理过程，除了表达各种元素的空间位置以外，同时还可以较好地表示三角形之间的拓扑关系，能够保持测量数据的原始性，且在表达地质界面的不连续方面具有优势，因而目前主要采用TIN模型来模拟地质界面。利用TIN模型来建模时，在模型简单的区域，三角形较少，而在模型复杂的区域，三角形较多。因此，TIN模型能较好地顾及模型特征，逼真表示模型的高低起伏变化，同时能够克服模型平坦区域的数据冗余。TIN的基本单元三角形的几何形状直接决定着TIN的应用质量。TIN模型由数据结构、三角划分准则和三角剖分算法等三部分组成。从结构上讲，TIN是一种矢量数据结构，对三角形的点、线、面的拓扑描述不同而形成不同的数据结构。近年来，借助于飞速发展的计算机软硬件技术，在TIN的快速建立、压缩存储、数据组织等方面均取得了突破性的进展[13-16]。TIN在地形分析[17]、空间统计分析[18]、几何分析[19]等空间分析中发挥着重要作用。

地质界面三维形态分析，以控矿因素定量提取为目的，以地质界面三维TIN模型为基础，利用空间几何原理计算地质界面的一般几何形态参数(坡度、夹角等)和距离场，利用空间插值技术和趋势剩余分析技术，提取出地质界面的局部形态特征和总体趋势，并以栅格模型来表达控矿地质因素场，实现各种控矿地质因素的定量模拟，为进一步建立矿体的立体定量预测模型提供基础。本研究提出并实现的地质界面三维形态分析方法包括：(1) 地质界面的几何形态参数提取；(2) 地质界面的形态趋势与起伏分析。

1.2 地质界面的几何形态参数提取

1.2.1 坡度提取

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

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

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

设三角形ABC的三个顶点的坐标为A(x1，y1，z1)，B(x2，y2，z2)，C(x3，y3，z3)，用向量法可以得出三角形所在平面的方程：

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

最后将式(4)和(5)代入式(2)就可以求出α，即对应于每个三角面片的坡度信息。通过距离场算法能获得地质空间中的立体单元到最近的三角面片的标示ID，该三角面片的坡度信息就是该立体单元所要赋值的坡度。

1.2.2 地质界面间夹角提取

求TIN模型中2个地层界面交线上的某点的夹角就是求经过该点的2个三角面片之间的夹角。但是，在TIN模型中，确定地质界面之间相交的三角面片难度较大，需要对2个界面上的三角面片逐一进行几何求交计算来获得实际相交的三角面片，寻找相交的三角面片的计算量与构成2个界面的三角面片的数量的乘积呈正相关，计算量大、时间效率低，所以本文提出了一种简单高效的方法来求取地质界面之间相交的三角面片，从而计算两地层界面之间的夹角。该方法首先将地质界面栅格化，获取三角面片的穿过的块体单元的网格索引，通过网格索引直接定位有可能实际相交的三角面片，然后采用向量法直接计算这些三角面片之间的夹角，最后通过赋值处理得到地质界面之间的夹角。

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

① 将三角面片的 3个顶点向坐标轴 XY平面投影，得到三角形为ABC(如图1所示)。

② 获取最小外包矩形。

③ 采用计算机图形学中的 Liang_Barsky参数化线段裁剪多边形算法[20]，计算三角形的三条边穿过的平面网格。

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

⑤ 重复这些步骤分别求得在XZ平面和YZ平面上相应的内部网格索引，然后将三个投影面上的网格索引的分量求交，最后求得三角面片通过的块体单元的索引，并且记录该三角形的ID。

⑥ 通过对所有三角面片循环计算，可以求出地质界面通过的块体单元，最后生成地质界面块体模型。

图1 TIN 模型中某个三角面片在XY平面上投影图Fig. 1 XY plane projection of a triangular of TIN model

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

2个三角形的夹角计算经常采用向量法，首先规定三角形的法向量向上，然后求2个三角形的法向量，法向量的夹角即为2个平面的夹角。

设 2个三角形的法向量分为 V1(x1，y1，z1)和V2(x2，y2，z2)，两法向量的夹角为α，则：

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

1.3 地质界面的形态趋势与起伏分析

在矿床地质空间中，相关地质体的主体形态将对周围成矿起到一定因素的影响。因此，对地质体形态的分析将是控矿因素定量提取的重要步骤。

地质体形态可以通过地质界面的波状起伏来描述。借鉴传统的曲面的趋势形态分析方法，结合地质界面的实际情况，本文采用地质界面的原始TIN模型，利用距离平方反比法对地质界面进行趋势形态分析。地质界面TIN模型的趋势形态分析实际上就是利用估测点周围一定距离之内的样品点对估测点的z重新赋值，然后利用计算得到的新的z生成地质界面的趋势面。具体过程如下。

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

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

(3) 利用计算得到的新的 z，建立新的地质界面的TIN模型，形成地质界面的趋势面。

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

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

对于地质曲面上的一个点 z(x0，y0)，可通过R(x，y)绝对值来描述该点所处的局部的形态起伏程度，通过R(x，y)的符号来描述该点所处的局部状态，有：

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

其中：R1(x，y)，R2(x，y)，…，Rn(x，y)表示n级局部起伏变化，各级局部变化的筛分可用趋势剩余分析方法。

利用上述方法可实现形态起伏的分级提取。因为插值搜索半径d决定可滤除的波形的大小，所以通过改变d可提取出不同级别的起伏。

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

地质界面TIN模型中任意三角形ABC(A(x1，y1，z1)，B(x2，y2，z2)，C(x3，y3，z3))，三角形 ABC 所在平面中的任意点M(x，y，z)，平面方程同式(3)。

可得：

其中：a=(x2-x1)(y3-y1)-(y2-y1)(x3-x1)，b=(z2-z1)×(y3-y1)(x-x3)+(y-y1)(x2-x1)(z3-z1)-(x-x1)(y2-y1)(z3-z1)-(y-y1)(z2-z1)(x3-x1)。

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

2 地质界面三维形态分析应用实例

实例：福建省尤溪县丁家山铅锌矿床的不整合面。

丁家山铅锌矿床是属于梅仙铅锌矿田的一个大型矿床。其编号为Ⅰ，Ⅱ和Ⅲ的三层矿体分布于龙北溪组上段(Pt2-3l3)绿片岩地层，产状平缓，矿体产状与含矿地层产状协调一致，矿床空间分布特征与该含矿岩性段总体延伸特征保持一致。在侏罗系长林组(J3c)地层中产出铅锌矿体(Ⅳ号矿体)，呈长条状或透镜状，矿体倾向与不整合面产状基本一致[21]。

龙北溪组上段绿片岩系为主要赋矿层位，燕山期石英斑岩、花岗斑岩脉在形成过程中矿区地层和含矿岩系是主要的控矿地质因素。某些铅锌矿体赋存于岩层构造裂隙中，矿体受断裂构造控制，产状与断裂构造产状一致，呈脉状，受裂隙构造控制的矿体的形成和分布可能与后期的热液叠加改造有关[21]。

不整合面在矿区内分布广泛，是控制矿体分布的主要构造。产于龙北溪组上段(Pt2-3l3)绿片岩地层中的矿体同时还受不整合构造控制。当不整合面较陡，且与下伏绿片岩地层层面呈较大角度斜交时，则在不整合面两侧一定范围内均有利于成矿。一方面，在不整合面构造两侧，矿体的产状受不整合面产状的影响，且矿体的厚度和品位有增大变富的趋势；另一方面，在不整合面构造中，还产出有充填型的脉状矿体。

从定性的角度，可以发现丁家山铅锌矿床不整合面的空间形态与成矿存在一定联系，然而对这种关系的认识是有限的，而且存在主观性，于是采用三维形态分析来对其进行定量化研究成为必然。又由于丁家山铅锌矿床不整合面属于正常产状，产状平缓，因此采用本方法进行形态分析。首先，采用Datamine软件建立地质界面三维TIN数据模型，并且建立矿化空间的栅格模型，然后利用研究开发的三维形态提取算法和程序，对不整合面地层界面进行三维形态分析，如对不整合面的形态趋势提取以及形态起伏分级提取,不整合面坡度提取，不整合面的夹角提取，不整合面距离场因素指标提取等，最后利用矿化空间栅格模型来表达各种控矿地质因素场。

图2所示为不整合面形态起伏分级提取。不整合面趋势-起伏因素分析主要是揭示不整合面起伏对不整合面周围地质空间的控矿作用影响。经过反复试验，选择100 m和200 m作为一级滤波和二级滤波的插值搜索范围半径，对不整合面原始TIN模型进行一级形态滤波和二级形态滤波，对应的得到不整合面的一级趋势与一级起伏(图 2(a))、二级趋势与二级起伏(图2(b))。

waU和 wbU分别代表不整合面一级起伏和二级起伏的起伏程度，用单元到不整合面的最近距离处的不整合面趋势部分的欧式距离 Einner或 Eouter来度量。显然，若单元到不整合面的最近距离处属于不整合面相应级别起伏的外凸部分，则waU和wbU大于0；若单元到不整合面的最近距离处属于不整合面相应级别起伏的内凹部分，则waU和wbU小于0。计算结果如图3所示。

陡倾不整合构造是层控型铅锌矿床后期改造最有利的构造，陡倾不整合构造旁侧矿体均强烈显示变厚变富的趋势，因此不整合面坡度因素gU的提取具有重要意义。计算结果的栅格模型表达如图4所示。

不整合面与下伏绿片岩地层层面呈较大角度斜交的这种对成矿有利的情况可以利用不整合面夹角因素aU_S来表达。不整合面与多个地层界面有相交，定义单元的不整合面夹角指标值取值为单元到不整合面的最近距离处对应的最近的夹角。计算结果如图5所示。陡倾不整合构造旁侧矿体均强烈显示变厚变富的趋势，并且产于火山岩盖层及其他层位的“充填型”矿体，受裂隙构造控制，这种构造因素可以用不整合面的距离场模型描述。从丁家山铅锌矿床不整合面三维TIN模型获取不整合面的线框模型，以单元到不整合面的最小距离作为距离场值来表达，即可得到不整合面距离场因素指标dU，如图6所示。为区分不整合面上下单元的距离场值的差异，将位于不整合面分布范围之上的单元场值置为正(正距离场)，位于不整合面分布范围之下的单元场值为负(负距离场)。

图2 不整合面形态起伏分级提取Fig. 2 Hierarchical extraction of undulate shapes of unconformity surface

图3 不整合面形态起伏因素(waU和wbU)分布栅格模型(-150～-200 m标高范围)Fig. 3 Raster model of shape factors of unconformity surface (waU and wbU)

图4 不整合面坡度因素(gU)分布栅格模型(-150～-200 m标高范围)Fig. 4 Raster model of slope factors of unconformity surface (gU)

图5 地层界面与不整合面的夹角因素(aU_S)分布栅格模型(-150～-200 m标高范围)Fig. 5 Raster model of angle factors of stratum interface and unconformity surface (aU_S)

图6 不整合面距离因素(dU)分布栅格模型(-150～-200 m标高范围)Fig. 6 Raster model of distance factors of unconformity surface (dU)

3 结论

(1) 结合三维地质建模、空间插值、地质场和空间分析等理论，提出了基于TIN的地质界面三维形态分析方法，通过该方法，可以快速有效的进行地质界面坡度和夹角等几何形态参数提取、地质界面距离场分析、地质界面形态趋势与起伏分析等，克服了栅格模型和数学形态学用于地质界面形态分析的缺点，是一种精度更高、用途更广泛的新方法。

(2) 以福建省尤溪县丁家山铅锌矿床为例，利用提出的方法对地层界面和不整合面等地质界面进行形态分析，提取了与铅锌成矿的各种地质形态因素，可以合理地定量化描述地质现象，进而在成矿预测中发挥作用，如陡倾不整合构造是成矿有利因素，而通过形态分析提取的不整合面坡度因素gU，可以便于定位到坡度陡倾的位置，或者便于建立坡度陡倾与成矿的非线性函数关系。

(3) 矿产资源的三维立体定量预测及三维评价系统是危机矿山可接替资源评价与找矿创新体系建立过程中亟待探索和实践的新理论、新方法。采用地质界面三维形态分析方法，以地质界面三维TIN模型的建立和可视化为基础，提取出地质界面空间形态指标，从而得到控矿地质因素场模型，进而用于控矿地质因素的定量分析，在隐伏矿体三维可视化立体定量预测中具有重要意义。

[1] 张善明, 吕新彪, 邓国祥, 等. 地质界面控矿原理及其运用要点[J]. 地质科技情报, 2009, 28(6): 51-56.ZHANG Shanming, LÜ Xinbiao, DENG Guoxiang, et al.Ore-controlling mechanism of geological interface and the key of application[J]. Geological Science and Technology Information, 2009, 28(6): 51-56.

[2] 朱志澄, 宋鸿林. 构造地质学[M]. 武汉: 中国地质大学出版社, 1990: 5-14.ZHU Zhicheng, SONG Honglin. Structural geology[M]. Wuhan:China University of Geosciences Press, 1990: 5-14.

[3] 邓浩. 面向隐伏矿体预测的三维地质建模与空间分析若干技术研究[D]. 长沙: 中南大学地球科学与信息物理学院, 2008:1-17.DENG Hao. Study on 3D geological modeling and spatial analysis for prediction of buried ore bodies[D]. Changsha:Central South University. School of Geosciences and Info-Physics, 2008: 1-17.

[4] 赵鹏大. 试论地质体的数学特征[J]. 地球科学, 1982(1):145-155.ZHAO Pengda. On the mathematical characteristics of geological bodies[J]. Earth Science-Journal of China University of Geoscience, 1982(1): 145-155.

[5] 周访滨, 刘学军. 基于栅格DEM自动划分微观地貌形态的研究[J]. 武汉理工大学学报: 信息与管理工程版, 2008, 30(2):172-175.ZHOU Fangbin, LIU Xuejun. Research on the automated classification of micro landform based on grid DEM[J]. Journal of Wuhan University of Technology: Information &Management Engineering, 2008, 30(2): 172-175.

[6] 姜栋, 赵文吉, 朱红春, 等. DEM地形信息提取对比研究: 以坡度为例[J]. 测绘科学, 2008, 33(5): 177-179.JIANG Dong, ZHAO Wenji, ZHU Hongchun, et al. Comparison of landform information extracted from DEMs: A case study of slope[J]. Science of Surveying and Mapping, 2008, 33(5):177-179.

[7] 毛先成, 唐艳华, 邓浩. 地质体的三维形态分析方法与应用[J]. 中南大学学报: 自然科学版, 2012, 43(2): 588-595.MAO Xiancheng, TANG Yanhua, DENG Hao.Three-dimensional morphological analysis method for geologic bodies and its application[J]. Journal of Central South University:Science and Technology, 2012, 43(2): 588-595.

[8] Serra J. Introduction to mathematical morphology[J]. Computer Vision, Graphics, and Image Processing, 1986, 35: 283-305.

[9] Serra J. Image analysis and mathematical morphology[M]. San Diego: Academic, 1982: 15-58.

[10] Jahjah M, Ulivieri C, Santilli G. Automatic archaeological feature extraction from satellite VHR images by mathematical morphology functions[C]//International Astronautical Federation-59th International Astronautical Congress 2008. Paris:IAF, 2008: 2751-2755.

[11] 吴立新, 史文中, Christopher G. 3DGIS与3DGMS中的空间构模技术[J]. 地理与地理信息科学, 2003, 19(1): 5-11.WU Lixin, SHI Wenzhong, Christopher G. Spatial modeling technologies for 3D GIS and 3D GMS[J]. Geography and Geo-Information Science, 2003, 19(1): 5-11.

[12] 彭仪普. 地形三维可视化及其实时绘制技术研究[D]. 成都:西南交通大学土木工程学院, 2002: 1-19.PENG Yipu. Studying on three-dimension visualization and real-time rendering of terrain[D]. Chengdu: Southwest Jiaotong University. School of Civil Engineering, 2002: 1-19.

[13] WU Hangbin, LIU Chun. Compression of discrete data based on the triangulate irregular network[C]//Proceedings of the SPIE-The International Society for Optical Engineering.Bellingham: SPIE, 2006: 6420.

[14] LU Linlin, GUO Huadong. Visualization of a digital elevation model[J]. Data Science Journal, 2007, 6: 481-484.

[15] Kothuri R, Ravada S, Fisher E. Triangulated irregular network:US, 12069089[P]. 2010-08-10.

[16] ZHU Qing, ZHANG Yeting, LI Fengchun. Three-dimensional TIN algorithm for digital terrain modeling[J]. Geo-spatial Information Science, 2008, 6(11): 79-85.

[17] YU Ming, CHEN Xiaoyu, AI Tinghua. Application research of terrain based on DEM and data mining[C]//CHEN Jingming, PU Yingxia. Geoinformatics 2007: Geosputial information Science.Bellingham, Wash: SPIE, 2007: 67531w.

[18] 钟登华, 郭享, 李明超, 等. 基于三维地质模型的地下洞室参数化设计与方案优选[J]. 天津大学学报, 2007, 40(5): 519-524.ZHONG Denghua, GUO Xiang, LI Mingchao, et al. Parametric design and schemes optimization for underground structure based on 3D geological model[J]. Journal of Tianjin University,2007, 40(5): 519-524.

[19] 许丽敏, 薛安. 基于Delaunay三角网与Voronoi图联合提取等高线骨架的地形重建算法研究[J]. 北京大学学报自然科学版,2009, 45(4): 647-652.XU Limin, XUE An. Terrain reconstruction from contours by skeleton extraction using Delaunay triangulation and Voronoi diagram[J]. Acta Scientiarum Naturalium Universitatis Pekinensis, 2009, 45(4): 647-652.

[20] 孙家广.计算机图形学[M].北京:清华大学出版社,1998:166-215.SUN Jiaguang. Computer graphics[M]. Beijing: Tsinghua University Press, 1998: 166-215.

[21] 吴建设, 黄仁生. 福建省尤溪峰岩铅锌银资源潜力及找矿方向探讨[J]. 中国地质, 2001, 28(12): 13-18.WU Jianshe, HUANG Rensheng. Potentials of lead-zinc and silver resources and their prospecting direction at Fengyan,Youxi, Fujian[J]. Chinese Geology, 2001, 28(12): 13-18.