基于GIS与层次分析法的综合成矿预测——以新疆库米什地区为例

doi:10.11872/j.issn.1005-2518.2020.02.139

基于GIS与层次分析法的综合成矿预测——以新疆库米什地区为例

陈超民^,¹^,², 冷成彪^,¹, 司国辉³

1. 东华理工大学核资源与环境国家重点实验室，江西南昌 330013

2. 上海师范大学环境与地理科学学院，上海 200234

3. 西安地质矿产勘查开发院有限公司，陕西西安 710100

Comprehensive Metallogenic Prediction Based on GIS and Analytic Hierarchy Process：A Case Study of Kumishi Region in Xinjiang

CHEN Chaomin^,¹^,², LENG Chengbiao^,¹, SI Guohui³

1. State Key Laboratory of Nuclear Resources and Environment，East China University of Technology，Nanchang 330013，Jiangxi，China

2. School of Environmental and Geographical Sciences，Shanghai Normal University，Shanghai 200234，China

3. Xi’an Institute of Geological and Mineral Exploration Co. ，Ltd. ，Xi’an 710100，Shaanxi，China

通讯作者: 冷成彪（1982-），男，山东临沭人，教授，从事矿床学方面的科研与教学工作。lcb8207@163.com

收稿日期: 2019-07-25 修回日期: 2019-11-25 网络出版日期: 2020-05-07

基金资助:

国家自然科学基金项目“新疆精河县色勒特果勒还原性斑岩—矽卡岩Cu-Mo矿的矿床地球化学研究”. 41872097

Received: 2019-07-25 Revised: 2019-11-25 Online: 2020-05-07

作者简介 About authors

陈超民（1995-），男，江西寻乌人，硕士研究生，从事遥感与GIS方面的研究工作ccm2222@163.com , E-mail：ccm2222@163.com

摘要

新疆库米什地区植被覆盖率低、基岩裸露度高，为遥感地质找矿提供了良好条件。为了提高该地区的找矿效率、实现找矿突破，系统收集了库米什地区的矿产、地质及遥感资料，并采用“掩膜+Crosta主成分分析+阈值分割”方法从ETM+遥感数据中提取了蚀变异常信息。在矿产、地质和遥感多源信息的基础上，总结出矿（化）点、岩体内外接触带、蚀变带、断层、羟基蚀变异常信息和铁染蚀变异常信息6个控矿因子，采用基于知识驱动的层次分析法建立了成矿预测模型，利用数学方法和GIS平台完成了综合成矿预测。最后，以部分未加入模型的矿点与野外实地考察结果，验证了成矿预测效果。结果表明：运用层次分析法在新疆库米什地区初步进行多源信息综合成矿预测，其结果具有一定的准确性，能够为该区进一步地质找矿工作提供参考。

关键词： 遥感 ; 蚀变异常提取 ; 主成分分析 ; 层次分析法 ; GIS ; 成矿预测 ; 新疆库米什

Abstract

The Kumishi region in Xinjiang has low vegetation coverage and high bedrock exposure，which are conducive to ore-prospecting by remote sensing.In order to improve the efficiency of ore-prospecting in this region and achieve a breakthrough in ore-prospecting，we systematically collected data of mine，geology and remote sensing in the Kumishi region.Then we extracted abnormal information of alterations from the ETM+ image by adopting the method of “Mask+Principal component analysis of Crosta+Threshold segmentation.”To check the accuracy of abnormal information of alterations，we compared the extracted alteration information with the field phenomenon of the wall-rock alteration in the Kalatage-Qigebu area and found that it was basically consistent.Based on multi-source information of mine，geology and remote sensing，six ore-controlling factors including mineralized points （ore spots），contact zones of magmatic rocks，alteration zones，faults，anomalies of hydroxyl alterations and anomalies of iron-stained alterations were selected.Since the number of deposits or ore spots in the Kumishi region was not enough to support the data-driven model，the analytic hierarchy process which belongs to a knowledge-driven decision-making method was adopted in this study.Consequently，the hierarchical model for the metallogenic prediction was established by the analytic hierarchy process which can present people’s subjective experience and thinking in digital form and realize the combination of qualitative analysis and quantitative analysis.In this model，the metallogenic prospective prediction was regarded as the target layer.The above six ore-controlling factors were considered as the criterion layer.0~300 m buffers，300~600 m buffers，and 600~1 000 m buffers were taken as the index layer in the lower layer of the first four criterion layers.The last two criterion layers were divided into three grades that were deemed as the index layer.By the combination of some previous research and expertise，some judgment matrices for each criterion layer and each index layer were constructed.In addition，we assigned values to each index layer.On this basis，some weights of each criterion layer and each index layer were calculated by the mathematical method.Meanwhile，the raster calculation and the kernel density analysis were carried out on the GIS platform to complete the comprehensive metallogenic prediction.Finally，some ore spots that were not added to the model were used to evaluate the effect of the metallogenic prediction.To test the accuracy of the predicted results，we drilled two holes in the Kalatage-Qigebu area.The scheelite mineralization in this area was preliminarily found because of the bluish-purple fluorescence of some samples under the ultraviolet light.It was confirmed more precisely by the scanning electron microscope（SEM） and the energy dispersive spectrometer（EDS）.In conclusion，the results show that the application of the analytic hierarchy process in the metallogenic prediction of multi-source information in the Kumishi region has a good effect.Therefore，it can provide some references for further work in this region.

Keywords： remote sensing ; extraction of alteration anomalies ; principal component analysis ; analytic hierarchy process ; GIS ; metallogenic prediction ; Kumishi region in Xinjiang

PDF (10443KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

陈超民, 冷成彪, 司国辉. 基于GIS与层次分析法的综合成矿预测——以新疆库米什地区为例[J]. 黄金科学技术, 2020, 28(2): 213-227 doi:10.11872/j.issn.1005-2518.2020.02.139

CHEN Chaomin, LENG Chengbiao, SI Guohui. Comprehensive Metallogenic Prediction Based on GIS and Analytic Hierarchy Process：A Case Study of Kumishi Region in Xinjiang[J]. Gold Science and Technology, 2020, 28(2): 213-227 doi:10.11872/j.issn.1005-2518.2020.02.139

库米什地区位于新疆吐鲁番托克逊县境内，目前该地区已发现彩虹铜多金属矿床^[1,2]、彩华沟铜多金属矿床^[3,4]、硫磺山多金属矿床^[4]、铜花山钴多金属矿床^[4]、天彩金矿床^[5]、桑树园子钨矿床^[6]和忠宝钨矿床^[7]等具有一定规模的矿床，显示出良好的找矿潜力。该地区海拔在700~2 300 m之间，地势总体上呈现出西北高、东南低的特征。相比传统找矿方法，遥感找矿技术具有宏观观测能力强、综合性强、信息量大和找矿效率高等优势。本研究区基岩裸露度高、植被覆盖率较低，为遥感地质找矿提供了良好的前提条件。

地质、矿产和遥感信息在成矿预测方面均有重要的指示意义。随着大数据时代的到来，基于GIS的成矿预测与评价逐渐完善^[8,9]。目前成矿预测模型主要有数据驱动和知识驱动2类，其中数据驱动模型主要研究矿床（点）与多源成矿信息之间的关系，通过定量统计分析建立数学模型进行评价，其结果取决于样本数据的详细程度；知识驱动模型主要是通过分析已有区域地质与成矿规律方面的资料，系统总结控矿因素并根据已有经验对其赋值，进而完成评价与预测^[10]。常用的数据驱动模型主要有证据权模型和逻辑回归模型^[11]。由于研究区尺度（比例尺为1/20万）内所掌握的矿床（点）数量及成矿信息不足以支撑数据驱动模型，故本研究采用一种知识驱动的决策评价方法——层次分析法，试图以简单直观的方式进行初步评价。该方法可以将人的经验和思维以数字化的形式呈现，实现定性与定量相结合。利用层次分析法可以将成矿预测问题分解成若干个影响成矿概率高低的控矿因素，根据前人总结的成矿规律和经验对各层次、各因素进行相对重要性的比较并赋予一定数值，通过数学方法计算各层次、各因素的权重值，最后在GIS平台上进行空间分析，以直观显示成矿预测结果。前人通过对地质、矿产、遥感、地球化学和地球物理等数据的研究，运用层次分析法在不同地区开展成矿预测，取得了良好的找矿效果^{[12,13,14,15,16,17,18,19]}。

本文在前人研究工作的基础上，采用“掩膜+Crosta主成分分析+阈值分割”方法从Landsat 7 ETM+遥感图像中提取铁染、羟基蚀变异常信息，以矿产、地质和遥感信息为依托，以GIS为平台，采用层次分析法在库米什地区进行多金属矿床的综合成矿预测与评价，结合预留的一些矿点及野外实地考察结果进行验证，旨在为该区后续找矿工作提供指导。

1 研究区地质概况

新疆库米什地区位于东天山吐—哈盆地的南缘（图1），该区横跨南天山、中天山和北天山3个性质不同的构造单元，被2条NWW向的逆冲断裂（中天山北缘断裂、中天山南缘断裂）分开^[20]。区内出露地层主要有元古宇中天山群星星峡组（Ptx）黑云母斜长片岩、眼球状片麻岩夹大理岩；志留系阿哈布拉克群（Sah）泥岩、粉砂岩、石英片岩、大理岩和凝灰岩；下泥盆统阿尔彼什麦布拉克组（D₁a）石英片岩和钙质片岩夹大理岩；中泥盆统阿拉塔格组（D₂a）石英片岩、钙质千枚岩、灰岩和大理岩；上泥盆统哈孜尔布拉克组（D₃h）粉砂岩、霏细岩和灰岩；下石炭统小热泉子组（C₁x）安山岩、斜长玢岩、凝灰砾岩和火山角砾岩；下石炭统牙曼苏组（C₁y）斜长玢岩、安山玢岩、闪长玢岩、凝灰砂岩和霏细斑岩；下石炭统马鞍桥组（C₁m）灰岩和砂岩；中石炭统桑树园组（C₂s）钙质砂岩、灰岩夹泥灰岩；中石炭统底坎尔组（C₂d）长石砂岩、凝灰岩、凝灰砾岩和凝灰砂岩；下二叠统阿其克布拉克组（P₁a）砾岩、砂岩、粉砂岩夹灰岩；下—中侏罗统三工河组（J_1-2s）砂岩、粉砂岩、泥质砂岩和黏土岩；第三系鄯善群（Esh）砾岩、砂岩和黏土岩；第三系桃树园组（N₁t）砾岩和砂岩；第三系葡萄沟组（N₂p）粗砂岩、砂砾岩和黏土质泥岩，以及第四系（Q）洪积堆积、风成堆积、湖相沉积和化学沉积。区内断裂发育，主要呈NE和NW向展布，区内中部发育有一处规模较大的背斜构造，其南部发育有一处规模较小的向斜构造，枢纽走向均为NWW向。区内存在以加里东期和华力西期花岗岩为主的多期次岩浆侵入活动，花岗岩呈岩基、岩株或岩枝状沿着研究区背斜构造的轴部和两翼侵位，构成NWW向花岗岩带^[20]。围岩蚀变主要有云英岩化、硅化、碳酸盐化、矽卡岩化、滑石化、黄钾铁矾化、绿帘石化、高岭土化、绿泥石化、钨矿化、褐铁矿化和绢云母化等。区内矿化现象较好，已发现一些矿化点^[21]和不同规模的矿床^[20]（如彩华沟铜多金属矿、桑树园子钨矿和忠宝钨矿等），具有良好的找矿前景。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 库米什地区地质简图（据文献[22]修改）

1.第四系；2.葡萄沟组；3.桃树园组；4.鄯善群；5.三工河组；6.阿其克布拉克组；7.底坎尔组；8.桑树园组；9.马鞍桥组；10.牙曼苏组；11.小热泉子组；12.哈孜尔布拉克组；13.阿拉塔格组；14.阿尔彼什麦布拉克组；15.阿哈布拉克群；16.中天山群星星峡组；17.加里东期变形花岗岩；18.华力西期花岗岩；19.辉橄岩；20.橄榄岩；21.断层；22.石膏矿点；23.石灰岩矿点；24.多金属矿点；25.褐铁矿化点；26.铁矿化点；27.铜矿（化）点；28.钨矿点

Fig.1 Brief geological map of the Kumishi region （modified after reference [22]）

2 遥感数据预处理

本研究采用Landsat 7 ETM+ SLC-on卫星数字产品（1999~2003）数据，数据下载自中国科学院计算机网络信息中心地理空间数据云平台（http://www.gscloud.cn），来源于美国地质调查局（United States Geological Survey，USGS）。ETM+数据共有8个波段，可见光波段（1、2、3）、近红外波段4和短波红外波段（5、7）的空间分辨率均为30 m，热红外波段6的空间分辨率为60 m，全色波段8的空间分辨率为15 m。本研究选取两景2002年10月18日无云的Landsat 7 ETM+遥感数据（轨道号分别为P142R30和P142R31）。

下载的产品类型为L1T级（Level 1T标准地形校正），采用UTM-WGS84南极洲极地投影和3次卷积取样。所下载的数据产品已经过系统的辐射校正、精确的地面控制点几何校正以及采用高精度数字高程模型（Digital Elevation Model，DEM）数据的地形校正。

在ENVI 5.3.1软件上对所下载的两景遥感图像进行辐射定标、FLAASH大气校正、多波段合成、图像拼接和研究区范围（88°00′~89°00′ E，42°00′~42°40′ N）的裁剪，以获得符合要求的遥感图像。

3 遥感蚀变异常信息提取

3.1 干扰信息的去除

遥感蚀变异常信息是一种弱信息^[23]，因此遥感蚀变信息的提取对遥感图像的要求较高，一般要求在冰雪覆盖较少、植被发育较弱的季节获取遥感数据，同时该时相的云量也要较少。但是受地形地貌的影响，有些因素仅靠数据时相的选择难以克服，比如阴影、水体、高山上的冰雪和第四系等，因此一般可以采用掩膜的方法来解决，使得这些较强的干扰信息不参与统计分析。由于研究区范围内图像无水体、云和冰雪，所以只需去除植被、阴影和第四系信息的干扰。

（1）植被的去除。由于归一化差值植被指数（Normalized Difference Vegetation Index，NDVI）是植被生长状态及植被覆盖度的最佳指示因子，且增强了对植被的响应能力^[24]，所以可通过NDVI来提取植被信息，从而完成剔除。计算公式如下：

N D V I = \frac{ρ_{N I R} - ρ_{R}}{ρ_{N I R} + ρ_{R}}

(1)

式中： $ρ_{N I R}$ 为近红外波段的反射率； $ρ_{R}$ 为红光波段的反射率。

根据标准假彩色合成图像中植被的特征，查询计算后的像元值，设定阈值为0.25，将像元值大于0.25的像元单独提取并进行掩膜。

（2）阴影的去除。本研究采用阴影植被指数（Shaded Vegetation Index，SVI）进行阴影的去除^[25]。SVI增强了明亮区域、阴影区域和水体之间的差异。计算公式如下：

S V I = N D V I \times ρ_{N I R}

(2)

根据真彩色合成图像中阴影的特征，查询计算后的像元值，设定阈值为0.002，将像元值小于0.002的像元单独提取并进行掩膜。

（3）第四系的去除。第四系为松散沉积和堆积物，存在较强烈的泥化现象，因此会对结果产生干扰，导致出现假异常^[26]。本研究根据1/20万新疆库米什地区地质简图（图1）对第四系进行掩膜。

3.2 4个波段特征向量主成分分析法（Crosta法）

Hunt等^[27]于20世纪70年代开始系统地测定各种岩石、矿物的可见光—近红外光谱，总结出一系列岩石、矿物在不同波段所表现的光谱特性和变化规律，这为矿化蚀变信息的提取奠定了基础。前人对于遥感蚀变信息的提取相继提出了众多行之有效的方法，如波段比值变换法、主成分分析法、光谱角填图法、对应分析法、Gram-Schmidt投影方法、混合像元分解法和MPH（Mask PCA and HIS）等方法^[23]。Crosta等^[28]于1989年采用PCA（TM1、TM3、TM4、TM5）和PCA（TM1、TM4、TM5、TM7）对巴西米纳斯吉拉斯州西部的亚热带地区残积土壤进行铁染、羟基蚀变信息填图。随后该技术被广泛应用于TM/ETM+遥感图像的铁染、羟基异常蚀变信息的提取，被称为Crosta技术。

（1）铁染蚀变异常信息的提取。根据含铁离子矿物的光谱特征（含铁离子矿物在ETM+图像中的波段1和4有较强的吸收带，在波段3存在反射峰）^[29]，在ENVI 5.3.1软件对掩膜后图像的波段1、3、4、5进行主成分分析，得到主成分分析特征向量和特征值表（表1）。能够代表铁染蚀变异常主分量图像的判断标准一般为：构成此主分量的特征向量，其ETM+3系数应该与ETM+1、ETM+4的系数符号相反，ETM+3一般与ETM+5系数符号相同^[30]。当铁染蚀变信息主分量中ETM+3、ETM+5对应的特征向量载荷因子符号为正，ETM+1、ETM+4对应的特征向量载荷因子符号为负时，铁染蚀变异常位于主分量的亮值区，否则位于暗值区^[31]。

表1 研究区掩膜后图像1、3、4、5波段主成分分析的特征向量及特征值

Table 1 Eigenvectors and eigenvalues of principal component analysis in band 1,3,4 and 5 of the masked image in the study area

主成分	ETM+1	ETM+3	ETM+4	ETM+5	特征值	信息量/%
PC1	0.271432	0.458413	0.523877	0.664632	4 531 404.269010	99.13
PC2	0.373363	0.434416	0.360416	-0.736195	30 152.496484	0.66
PC3	0.804067	0.004837	-0.580966	0.126218	8 510.120933	0.19
PC4	-0.374701	0.775318	-0.508066	0.018739	908.433548	0.02

新窗口打开| 下载CSV

由表1可知，PC4中波段3与波段1、波段4反差较大，符合判断标准，也与含铁矿物在ETM+1和ETM+4上具有较强吸收特征、在ETM+3上具有高反射特征正好相符。因此，PC4可作为铁染蚀变异常的指示主分量。

（2）羟基蚀变异常信息的提取。根据黏土类矿物（含羟基矿物）光谱特征（含羟基矿物在ETM+图像中的波段7表现为特征吸收带，在波段5呈现相对高反射）^[29]，在ENVI 5.3.1软件上对掩膜后图像的波段1、4、5、7进行主成分分析，得到主成分分析特征向量和特征值表（表2）。能够代表羟基蚀变异常主分量图像的判断标准一般为：构成此主分量的特征向量，其ETM+5系数应该与ETM+4、ETM+7的系数符号相反，ETM+5一般与ETM+1系数符号相同^[30]。当羟基蚀变信息主分量中ETM+1、ETM+5对应的特征向量载荷因子符号为正，ETM+4、ETM+7对应的特征向量载荷因子符号为负时，羟基蚀变异常位于主分量的亮值区，否则位于暗值区^[31]。

表2 研究区掩膜后图像1、4、5、7波段主成分分析的特征向量及特征值

Table 2 Eigenvectors and eigenvalues of principal component analysis in band 1,4,5 and 7 of the masked image in the study area

主成分	ETM+1	ETM+4	ETM+5	ETM+7	特征值	信息量/%
PC1	0.246298	0.475324	0.607979	0.586316	5 441 535.328711	98.90
PC2	0.359951	0.686251	-0.077387	-0.627301	41 000.148234	0.75
PC3	0.679599	-0.060805	-0.612567	0.399010	13 197.824440	0.24
PC4	-0.589846	0.547202	-0.499134	0.321742	6 266.013310	0.11

新窗口打开| 下载CSV

由表2可知，PC4中波段5与波段7反差较大，符合判断标准，也与黏土类矿物（含羟基矿物）在ETM+7具有较强吸收特征、在ETM+5上具有高反射特征正好相反。因此，PC4可作为羟基蚀变异常的指示主分量。

3.3 阈值分割

在阈值分割之前，对图像分别进行5×5均值滤波处理，排除孤立的散点，从而有效地反映蚀变异常信息的分布情况。由于ETM+数据及其线性处理结果近似服从正态分布，因此可以采用门限阈值进行分级处理，即以“均值+N×标准偏差”来划分异常强度的等级^[30,32]。铁染蚀变异常和羟基蚀变异常中的N值均设定为1.5，2.0，2.5，依此将异常强度等级划分为三级、二级和一级（图2和图3）。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 铁染蚀变异常分布（由ETM+遥感图像真彩色合成叠加）

1.铁染三级蚀变异常；2.铁染二级蚀变异常；3.铁染一级蚀变异常

Fig.2 Distribution of anomalies of iron-stained alterations（superimposed with the true color synthesis for ETM+ image）

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 羟基蚀变异常分布（由ETM+遥感图像真彩色合成叠加）

1.羟基三级蚀变异常；2.羟基二级蚀变异常；3.羟基一级蚀变异常

Fig.3 Distribution of anomalies of hydroxyl alterations(superimposed with the true color synthesis for ETM+ image）

3.4 蚀变异常信息验证

将遥感蚀变信息应用于成矿预测的地质依据是热液型金属矿床的围岩蚀变。当围岩受到成矿热液作用时围岩组分会发生变化，围岩蚀变便是成矿物质在逐步富集过程中留下的印记，因此可将其作为热液型金属矿床的找矿标志^[33]。一般来说，热液型金属矿床的出现会伴随着不同程度的围岩蚀变，但并非所有的围岩蚀变现象均与矿床的出现有关，因此通过遥感手段提取的蚀变异常信息一般起着指示作用。此外，异物同谱现象的存在也会引起蚀变异常信息的偏差。为了解释蚀变异常信息与矿床勘查之间的关系，以彩华沟铜多金属矿床^[3,4]和喀拉塔格—其格布查证区为例进行说明。

（1）彩华沟铜多金属矿床。该矿床位于研究区的东部（图1东部铜矿点所处位置），赋存于下泥盆统阿尔彼什麦布拉克组（D₁a），该矿床的围岩蚀变较强，主要蚀变类型有硅化、绢云母化、含铁碳酸盐化和绿泥石化，其次为石膏化和重晶石化。然而，由图2和图3可见，该矿床所处的区域并未呈现出应有的强烈蚀变。

（2）喀拉塔格—其格布查证区。本研究选取了ETM+遥感图像的7，4，1波段进行假彩色合成［图4（a）］，该波段组合具有兼并可见光、近红外和短波红外波段信息的优势，图像具有较丰富的色彩和地质信息。本次查证地段位于喀拉塔格—其格布地区，查证地段地层属元古宇中天山群星星峡组（Ptx），岩性主要为黑云母斜长片岩和眼球状片麻岩，在图像上的特征表现为红褐色、NW向带状展布。查证地段的花岗岩在图像上的特征表现为浅绿色、NW向带状展布。由图4可知，在该地段铁染蚀变信息比羟基蚀变信息少得多，这与野外实际情况相吻合。经查证，该地段发现有大量绢英岩化，而褐铁矿化却很少。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 喀拉塔格—其格布区蚀变异常信息野外验证结果

（a）ETM+7，4，1波段合成图像；（b）铁染蚀变异常信息；（c）羟基蚀变异常信息；（d）地貌景观；（e）面状分布的褐铁矿化与绢英岩化；（f）条带状分布的褐铁矿化与绢英岩化

Fig.4 Field validations for abnormal information of alterations in the Kalatage-Qigebu area

由上述可知：利用Crosta主成分分析法从ETM+遥感图像中提取的蚀变异常信息基本上是具备指示意义的，但考虑到其在部分地区存在偏差，因此在后续层次分析法中并没有赋予蚀变信息太高的权重。未来可以引入高分辨率和高光谱等遥感数据，这将有助于提高蚀变异常信息的准确性。

4 基于层次分析法的成矿预测

层次分析法（Analytic Hierarchy Process，AHP）是由Saaty教授^[34]于20世纪70年代提出的一种简便、灵活且实用的多准则决策分析方法，该方法不断地被完善和推广^[35,36]，随后被国内引进并作了系统的论述^[37,38]。通过层次分析法，可以将成矿预测分析问题分解成若干个影响成矿可能性高低的因素（本研究根据所掌握的资料选取了与成矿关系密切的矿产、地质和遥感信息作为预测因子），将其按照支配关系形成层次结构，两两比较，构建判断矩阵，求解各控矿因子的权重，建立成矿预测数学模型进行预测。

4.1 构建层次结构模型

结合前人研究和所掌握的资料数据，以科学性、可比较性、可度量性和可实施性筛选预测因子，建立层次结构模型（图5），其中A为目标层，B1~B6为准则层（B1为矿产因子，B2~B4为地质因子，B5和B6为遥感因子），C1~C18为指标层。值得注意的是，地质因子中的蚀变带（B3）和遥感因子中的羟基、铁染蚀变异常信息（B5、B6）之间并不是完全独立的，二者存在部分重叠，但考虑到地质因子中的蚀变带是已被证实的，且其数量和范围远远小于遥感技术提取的蚀变异常信息，因此蚀变带引入所造成的冗余度很小，有助于增强预测结果的准确性。此外，指标层中的缓冲区阈值的设定，对于预测结果具有十分重要的意义，阈值偏大将很难体现预测结果的价值，而阈值偏小则很难开展初步预测。一般而言，阈值是基于经验而设定的，本文采用何珊等^[19]的缓冲区阈值进行划分。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 成矿预测层次结构图

Fig.5 Hierarchical chart for the metallogenic prediction

倘若矿床（点）数据足够多，可以统计矿床（点）与各控矿因素缓冲区距离之间的关系，从而找出合适的缓冲区阈值。许多学者在使用证据权模型时统计了矿床点与控矿因素之间的距离关系，其结果随地区与样本数据而异，但基本都处在4 km范围内^{[39,40,41,42]}。

4.2 各控矿指标的分级与赋值

各控矿因子的量纲不一，若要进行成矿预测分析，则需按照一定的标准对各控矿因子进行分级赋值。B1~B4中将0~300 m缓冲区赋值为6，300~600 m缓冲区赋值为4，600~1 000 m缓冲区赋值为2，大于1 000 m区域赋值为0。B5和B6中将一级蚀变异常区赋值为6，二级蚀变异常区赋值为4，三级蚀变异常区赋值为2，无蚀变异常区赋值为0。值得注意的是，当所赋具体数值在同一个标准等级内时，具体的数值大小并不会对结果产生影响，因为对最终计算所得的结果而言，其数值间的差异才是关注点。

4.3 构建判断矩阵

判断矩阵是表示本层所有的因子对上一层某个因子相对重要性的比较，这是层次分析法的关键^[37]。通常采用9标度作为判断标准，标度代表着要素k相对于要素l对上层要素的重要程度，其中1表示k与l同等重要，3表示k比l稍微重要，5表示k比l比较重要，7表示k比l强烈重要，9表示k比l极端重要^[38]。针对目标层A，准则层B中各因素相互比较，并赋予标度，从而构建准则层判断矩阵A-B。针对准则层B_j，指标层C中各因素相互比较，并赋予标度，从而构建6个指标层判断矩阵B_j-C（j=1,2,…,6）。因此结合前人研究资料与分析者的判断构建出7个比较判断矩阵：

A - B = [\begin{matrix} 1 & 2 & 4 & 6 & 8 & 9 \\ 1 / 2 & 1 & 2 & 4 & 6 & 8 \\ 1 / 4 & 1 / 2 & 1 & 2 & 4 & 5 \\ 1 / 6 & 1 / 4 & 1 / 2 & 1 & 2 & 3 \\ 1 / 8 & 1 / 6 & 1 / 4 & 1 / 2 & 1 & 2 \\ 1 / 9 & 1 / 8 & 1 / 5 & 1 / 3 & 1 / 2 & 1 \end{matrix}]

（3）

B_{j} - C = [\begin{matrix} 1 & 3 & 5 \\ 1 / 3 & 1 & 2 \\ 1 / 5 & 1 / 2 & 1 \end{matrix}] (j = 1,2, \dots, 6)

（4）

4.4 层次单排序及一致性检验

层次单排序的过程是计算判断矩阵特征值与特征向量的过程，即对于判断矩阵A，计算满足 $A W = λ_{m a x} W$ 的特征值与特征向量。其中， $λ_{m a x}$ 为判断矩阵A的最大特征值，W为对应于 $λ_{m a x}$ 的归一化特征向量，W的分量便是各对应元素层次单排序的权重值^[37]。

由于客观事物的复杂性与人类认识的多样性，构造出的判断矩阵难免包含一定的非一致性，若非一致性较高时，继续把排序的权重向量结果作为决策依据，就会出现问题，因此在得到 $λ_{m a x}$ 后需要进行一致性检验。步骤如下：

（1）计算一致性指标CI，即

C I = (λ_{m a x} - n) / (n - 1)

(5)

式中： $λ_{m a x}$ 为判断矩阵最大特征值；n为判断矩阵的阶数。

（2）查找相应的平均随机一致性指标RI。Saaty^[38]最初通过求取500个随机判断矩阵的平均特征值给出了1到10阶矩阵的RI值；许树柏^[38]对其进行改进，将随机样本增加到1 000个，给出了1到15阶矩阵的RI值；Saaty^[36]将随机样本个数增加到50 000个，给出了1到10阶矩阵的RI值。本文采用50 000次随机重复试验所得的RI值，通过查表获得符合本文判断矩阵阶数的RI值：当矩阵阶数n=6时，RI为1.25；当矩阵阶数n=3时，RI为0.52。

（3）计算一致性比例CR，即

C R = C I / R I

(6)

当CR<0.10时，认为判断矩阵具有满意的一致性，否则就需要对判断矩阵进行调整。

对于上述各比较判断矩阵，用MATLAB数学软件进行计算，结果如表3所示。

表3 各判断矩阵的计算结果

Table 3 Results of the calculation for each comparison matrix

判断矩阵	层次单排序的权重向量W	$λ_{m a x}$	CI	RI	CR
A-B	（0.4361，0.2639，0.1446，0.0779，0.0463，0.0313）^T	6.1023	0.0205	1.2500	0.0164
B₁-C	（0.6483，0.2297，0.1220，0，0，0，0，0，0，0，0，0，0，0，0，0，0，0）^T	3.0037	0.0018	0.5200	0.0035
B₂-C	（0，0，0，0.6483，0.2297，0.1220，0，0，0，0，0，0，0，0，0，0，0，0）^T	3.0037	0.0018	0.5200	0.0035
B₃-C	（0，0，0，0，0，0，0.6483，0.2297，0.1220，0，0，0，0，0，0，0，0，0）^T	3.0037	0.0018	0.5200	0.0035
B₄-C	（0，0，0，0，0，0，0，0，0，0.6483，0.2297，0.1220，0，0，0，0，0，0）^T	3.0037	0.0018	0.5200	0.0035
B₅-C	（0，0，0，0，0，0，0，0，0，0，0，0，0.6483，0.2297，0.1220，0，0，0）^T	3.0037	0.0018	0.5200	0.0035
B₆-C	（0，0，0，0，0，0，0，0，0，0，0，0，0，0，0，0.6483，0.2297，0.1220）^T	3.0037	0.0018	0.5200	0.0035

新窗口打开| 下载CSV

4.5 层次总排序及一致性检验

层次总排序的过程是计算针对上一层次而言的本层次所有元素重要性的权重值，这一过程是从上到下逐层顺序进行的，同时需要利用同一层次中所有层次单排序的结果^[37]。设B层中有m个元素（B₁，B₂，…，B_m），其关于上一层A的层次总排序权重分别是b₁，b₂，…，b_m，C层中有n个元素（C₁，C₂，…，C_n），这n个元素关于B层中元素B_j（j=1，2，…，m）的层次单排序权重分别是c₁_j，c₂_j，…，c_nj，C层中各元素的层次总排序权重值（c₁，c₂，…，c_n）计算公式如下：

c_{i} = \overset{m}{\sum_{j = 1}} c_{i j} b_{j}, i = 1,2, \dots, n

(7)

为了评价层次总排序结果的一致性，同样需要进行一致性检验，因为各层次的非一致性可能累积并影响到最终的结果^[37]。设C层中元素（C₁，C₂，…，C_n）关于B层中元素B_j（j=1，2，…，m）的层次单排序一致性指标为CI_j（j=1，2，…，m），相应的平均随机一致性指标为RI_j（j=1，2，…，m），则C层的层次总排序一致性比例CR_final为

C R_{f i n a l} = \overset{m}{\sum_{j = 1}} C I_{j} b_{j} / \overset{m}{\sum_{j = 1}} R I_{j} b_{j}

(8)

根据表3结果，在MATLAB数学软件上进行计算得到层次总排序结果（表4）和层次总排序的一致性比例CR_final=0.00346<0.10，符合一致性要求。

表4 层次总排序结果

Table 4 Results of the general hierarchical ranking

准则层	准则层权值b_j	指标层	指标层权值c_i
矿化点（B1）	0.4361	0~300 m缓冲区（C1）	0.2827
		300~600 m缓冲区（C2）	0.1002
		600~1 000 m缓冲区（C3）	0.0532
岩体内外接触带（B2）	0.2639	0~300 m缓冲区（C4）	0.1711
		300~600 m缓冲区（C5）	0.0606
		600~1 000 m缓冲区（C6）	0.0322
蚀变带（B3）	0.1446	0~300 m缓冲区（C7）	0.0937
		300~600 m缓冲区（C8）	0.0332
		600~1 000 m缓冲区（C9）	0.0176
断层（B4）	0.0779	0~300 m缓冲区（C10）	0.0505
		300~600 m缓冲区（C11）	0.0179
		600~1 000 m缓冲区（C12）	0.0095
羟基蚀变异常信息（B5）	0.0463	一级蚀变异常（C13）	0.0300
		二级蚀变异常（C14）	0.0106
		三级蚀变异常（C15）	0.0056
铁染蚀变异常信息（B6）	0.0313	一级蚀变异常（C16）	0.0203
		二级蚀变异常（C17）	0.0072
		三级蚀变异常（C18）	0.0038

新窗口打开| 下载CSV

4.6 构建成矿预测数学模型

在构建成矿预测模型中引入成矿有利度Y进行评价，模型如下：

Y = \overset{n}{\sum_{i = 1}} c_{i} F_{i}

(9)

式中：i为控矿指标的个数；F_i为第i个控矿指标的具体赋值；c_i为第i个控矿指标对应的排序权重值（表4）。

4.7 基于ArcGIS平台的空间分析

将层次结构模型中的各类控矿因子集中在ArcGIS 10.3.1软件（图6）。随后对矿（化）点、岩体内外接触带、蚀变带和断层以图5的阈值进行缓冲区处理，将C1~C18各控矿指标按照成矿预测数学模型进行栅格计算得到成矿有利度图（图7）。为了探讨成矿有利点的集聚程度，将成矿有利度大于0的栅格均转换成点进行核密度分析（图8），ArcGIS中的核密度分析以四次核函数为基础，采用基于Silverman经验规则的带宽搜索算法计算默认带宽。虽然该处理忽略了成矿有利点有利度大小的差异，但能够体现成矿有利点的分布规律。将成矿有利度图与成矿有利点核密度分布图结合起来考虑效果较好。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 各类控矿因子叠加图

1.总矿（化）点；2.断层；3.总蚀变带；4.总岩体；5.羟基三级蚀变异常；6.羟基二级蚀变异常；7.羟基一级蚀变异常；8.铁染三级蚀变异常；9.铁染二级蚀变异常；10.铁染一级蚀变异常

Fig.6 Stacking map of various ore-controlling factors

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 成矿有利度图

Fig.7 Map of the ore-forming favorability

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 成矿有利点核密度分布图

1.模型中使用的矿（化）点；2.预留的矿点；3.验证的矿化点；4.喀拉塔格—其格布项目区

Fig.8 Kernel density distribution map of favorable metallogenic points

5 成矿预测效果评价

图7展示了成矿有利度的分布情况，值越高（颜色越深）代表着该区域越有利于成矿。图8反映了成矿有利点的集聚程度与分布规律，核密度越高代表该区域集聚了较多的成矿有利点。通过对比图7和图8可知，成矿有利度与成矿有利点核密度之间存在一定的关系，即成矿有利度高，其成矿有利点核密度一般也高；但成矿有利点核密度高并不意味着成矿有利度高，这是合乎逻辑的，因为核密度处理过程忽略了各点的成矿有利度。此外，本研究的模型中共使用了18个矿（化）点数据，预留了4个矿点数据作为检验，并在喀拉塔格—其格布地区考察了2个点作为验证。由图7和图8可知，作为预留检验的4个矿点基本处在成矿有利度和成矿有利点核密度较高的位置。但即便是在已使用的矿（化）点数据中也存在一些点并未处于成矿有利度和成矿有利点核密度较高的位置，可能有2个方面的原因：一是其本身与模型中控矿因素的外在关系较弱；二是模型中缺乏其他控矿因素信息。

为了探讨矿（化）点以外的区域中成矿预测的准确性，本研究选取了喀拉塔格—其格布地区进行考察。喀拉塔格—其格布地区位于新疆托克逊县库米什镇北约20 km处（图8）。在该地区实施了2个钻孔（位置如图8中验证的矿化点所示），经查证发现该区存在白钨矿化现象，主要分布于花岗岩的节理裂隙面、花岗岩与石英脉的接触面、花岗岩表面的石英薄膜上、片麻岩的片麻理表面上以及片麻岩与石英脉的接触面。白钨矿主要呈星点状、薄膜状和浸染状产出。白钨矿具有荧光性，其在紫外灯的照射下呈现蓝紫色（图9）。此外，某些样品的扫描电子显微镜—能谱仪（Scanning Electron Microscope-Energy Dispersive Spectrometer，SEM-EDS）分析结果也证明了白钨矿的存在（图10）。在成矿预测效果评价中仅实地考察了喀拉塔格—其格布地区，考察结果表现出一定的准确性，但成矿预测结果在其他未经验证区域的准确性仍需进一步的检验。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 喀拉塔格—其格布地区的白钨矿化

（a）片麻岩与石英脉接触处的白钨矿化；（b）片麻岩与石英脉接触处的白钨矿化（在紫外灯照射下）；（c）、（d）钻孔岩芯中花岗岩的白钨矿化（在紫外灯照射下）

Fig.9 Scheelite mineralization in the Kalatage-Qigebu area

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 片麻状黑云母花岗岩样品SEM-EDS结果（25.0 kV，×100）

Fig.10 SEM-EDS result of the gneissic biotite granite sample （25.0 kV，×100）

6 结论

本研究以前人地质矿产资料和Landsat 7 ETM+数据为基础，总结出矿（化）点、岩体内外接触带、蚀变带和断层4个控矿因子。通过“掩膜+Crosta主成分分析+阈值分割”法从ETM+遥感数据中提取了羟基、铁染蚀变异常信息，并在喀拉塔格—其格布地区进行了考察，发现考察区的围岩蚀变现象与通过遥感手段提取的蚀变异常信息基本吻合。因此，提取的蚀变异常信息除了部分区域存在偏差外基本具备指示意义。将矿（化）点（矿产信息）、岩体内外接触带、蚀变带、断层（地质信息）和羟基、铁染蚀变异常信息（遥感信息）作为控矿判断因子，以GIS为平台，采用基于知识驱动的层次分析法建立数学模型进行成矿预测。在成矿预测效果评价的过程中，除了利用模型中已使用的矿（化）点和预留矿点外，还在喀拉塔格—其格布地区进行了钻孔设计，发现考察区存在钨矿化，表明本研究具有一定的参考价值，能够对该区进一步找矿工作起到指示作用。

层次分析法是一种将主观经验与知识定量化的决策评价分析方法，可将人的主观思维以数字化的形式体现，借助GIS平台能够实现可视化。而数据驱动方法一般基于统计学方法，其结果取决于样本数据的详细程度，随着大数据时代的到来，基于数据驱动方法的优势愈加凸显。倘若能够整合同一尺度下的矿产、地质、遥感和物化探数据，便可以定量揭示隐藏在数据之中的一般规律，从而将传统的定性的成矿预测问题朝着定量方向转变，定量化的预测结果在一定程度上也能反作用于定性结果。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

李爱民，高保明，路枫，等 .

新疆托克逊县彩虹铜多金属矿床地质特征

[J].资源环境与工程，2011，25（4）：293-298.