基于MAPGIS空间分析的长坑幅信息量法成矿远景预测

卓福星

(福建省闽东南地质大队，泉州，362021)

成矿预测涉及地质、物理、化学、遥感等多源多尺度地学数据，为有效管理与充分应用多源、多尺度地学数据，许多研究者将GIS(地理信息系统)应用于成矿预测中，如在地理信息系统软件MAPGIS基础上，研制了矿产资源综合信息评价系统MRAS(中国地质科学院矿产资源研究所)[1]、金属矿产资源评价分析系统MORPAS(中国地质大学(武汉)数学地质研究所、中国地质调查局发展研究中心等)[2]，应用证据加权法进行攀西地区铜镍铂族元素找矿远景预测(2009)[3]、祁连山地区铜矿成矿预测(2009)[4]，通过模糊逻辑法优选秦岭地区铅锌矿找矿远景区(2009)[5]、矿化系数等值线法预测黑龙江金矿(2007)[6]，及采用MORPAS特征分析法的云南个旧西区锡多金属成矿远景预测(2009)[7]等；基于ArcGIS或ArcView等GIS软件，开发ArcSDM、ArcWofe、GeoDAS等成矿预测模块或软件[8，9]，应用证据权法的成矿远景预测(2007)[10]，通过权重叠加模型的成矿预测(李随民等，2008)[11]等。虽然我国地勘工作中也广泛使用了GIS，但是多应用于制图、成果管理等，较少涉及空间分析。笔者以信息量法进行成矿预测，并以1∶5万长坑幅为实验区，探讨通过数据驱动方法在MAPGIS中实现综合信息定量预测的空间分析建模方法。

1∶5万长坑幅实验区位于南武夷晚古生代坳陷区大田—龙岩坳陷带东侧，闽东火山断坳带上，政和—大埔与福安—南靖断裂带之间，闽江口—永安和漳平—仙游断裂带交会部，属于洛阳—剑斗东西向铁多金属成矿带(东段)。区内构造、岩浆活动强烈，具多期次叠加改造特点，铁、铜、铅锌成矿地质条件优越，已发现铁、铜、钼、铅锌多金属矿床(点)17处，其中内生铁、铜、铅锌矿12处。

1 MAPGIS成矿预测的空间数据建立

根据福建漳平洛阳—安溪剑斗地区区域地质矿产调查*福建省闽东南地质大队，福建漳平洛阳-安溪剑斗地区区域地质矿产调查报告 (1∶5万草测)，1986。、福建漳平洛阳—安溪剑斗铁多金属成矿预测*福建省闽东南地质大队，福建漳平洛阳-安溪剑斗铁多金属成矿预测 (1∶5万)，1998。、福建漳平洛阳—安溪潘田地区矿产远景调查*福建省地质调查研究院，福建漳平洛阳—安溪潘田地区地质矿产远景调查报告 (1∶5万)，2014。等成果，对实验区1∶5万长坑幅铁、铜、铅锌矿的成矿背景，典型矿床、成矿条件、控矿因素、找矿标志等进行深入的分析研究。在此基础上，以福建漳平洛阳—安溪潘田地区地质矿产图(1∶5万)为基本数据源，根据试验区综合信息成矿预测目标、MAPGIS 6.6空间数据组织特点，及相关规范、要求、标准，进行数据分类分层采集，建立基于MAPGIS 6.6的空间数据库。此外，对与该区铁、铜、铅锌矿成矿预测具比较重要意义的矽卡岩化、角岩化、硅化等围岩蚀变不能从1∶5万地质矿产图获取(虽图中有各种蚀变符号，但不能进行准确定位相应蚀变的空间分布)，以及考虑到保持整个试验区的信息水平的一致性，因此忽略了蚀变图层的提取，仅是在划分远景区等级时予以充分考虑。

基本单元.wp是根据该区工作程度、矿床/点空间格局、控矿构造特征及1∶5万成矿预测要求等[12]，在MAPGIS 6.6输入编辑模块中，通过造平行线方法，按东西向、北南向以1 250 m×1 250 m划分网格线，拓扑处理后，经条件检索提取面积≥1 041 666.7 m2(基本网格面积的2/3)的多边形为成矿预测基本单元，共280个单元。

2 MAPGIS空间分析的找矿信息量提取

信息量法、证据权法等综合信息成矿远景定量预测[13-15]属综合评价问题，包括预测/评价指标(标志/变量)体系建设、指标赋值与规范化、指标赋权、多指标综合等环节。因此，在典型矿床研究及区域成矿条件、控矿因素、找矿标志等综合分析基础，根据福建漳平洛阳—安溪潘田地区地质矿产图为基本数据源构建的空间数据库，利用数据驱动方法通过MAPGIS空间分析提取找矿信息。

2.1 地层侵入岩的找矿信息量

2.1.1 围岩标志的初选

在空间分析模块，分别对侵入岩.wp、地层.wp进行点对面相交分析，并对派生的点文件进行双属性分类统计(图1)。南园组、童子岩组、林地组、文笔山组及燕山早期第二、三次与晚期第一次侵入等地层、侵入岩与区内铁、铜、铅锌等矿床/点具空间关联性，可能为有利的围岩，其有利程度或重要性可通过矿化度/单位面积见矿率等测度，笔者采用找矿信息量确认。

图1 围岩见矿(铁、铜、铅锌)频数直方图Fig.1 Freguency histogram of wall vocks beering iron,copper,lead and zinc ore A1—南园组；A2—童子岩组；A3—文笔山组；A4—栖霞组；A5—林地组；A6—大岭组；A10—燕山晚期第一次侵入；A12—燕山早期第三次侵入；A13—燕山早期第二次侵入

2.1.2 单元的地层、侵入岩存在状态统计

在空间分析模块：① 通过条件检索，由地层.wp的组属性分别检出南园组、童子岩组、林地组、文笔山组、栖霞组等地层，由侵入岩.wp的期、次属性分别检出燕山早期第二、三次及晚期第一次侵入的侵入岩，从而派生地层、侵入岩数据，再通过文件合并或地层.wp 的系、统、侵入岩.wp的期分别派生A7、A8、A9及A13。②通过已知有矿单元点.wt、基本单元点.wt分别与前述派生数据(A1～A13)依次进行点对面相交分析，派生Bi、Di(下标i=1,2…13，与A下标对应，为单元中心点落入Ai对象/标志的单元集合。③以单元中心是否位于Ai对象中作为识别单元是否存在该标志Ai的准则，则通过Bi、Di的属性统计分别得已知有矿单元、基本单元中存在Ai的单元数Ni,j、Si,j，即存在Ai时的有矿单元数Ni,j、单元总数Si,j。找矿信息图层及其信息量统计(表1)。

表1 找矿信息图层(文件)及其信息量统计

续表1

注：B表示已知有矿单元点.wt，D表示基本单元点.wt，j表示“存在”；Bi=B·AND·Ai，Di=D·AND·Ai，Di-=D·NOT·Ai，AND、NOT分别表示相交、相减叠置。

2.1.3 信息量算法与地层侵入岩的找矿信息

某种标志/变量对研究对象的作用由信息量测度，找矿标志/变量的信息量表示对找矿的贡献大小，通过条件概率计算：

IAi,j→B=lg(P(B/Ai，j)/P(B))

(1)

式中：IAi，j→B为Ai标志j状态提供事件B发生的信息量(Ai的找矿贡献)，P(B/Ai,j)为Ai标志j状态存在条件下事件B实现的概率(存在Ai时的见矿概率)，P(B) 为事件B发生的概率(研究区总的有矿概率)。

由于P(B)在工作初期不易估计，及概率一般是通过样本频率估计的，因此，根据概率乘法原理，式(1)进一步等价变换：

IAi,j→B=lg(P(Ai，j/B)/P(Ai，j))=lg((Ni,j/N)/(Si,j/S))

(2)

式中：P(Ai,j/B)为已知事件B发生条件下Ai存在的概率，有矿单元Ai存在的概率，由已知有矿单元中存在Ai的频率Ni,j/N估计。根据式(2)提取的地层、侵入岩的找矿信息量。

2.2 物化探异常的找矿信息

基于MAPGIS空间分析的物化探异常找矿信息量统计方法、步骤与前述地层侵入岩的找矿信息量统计过程基本一致，因此，仅以土壤金属量异常以例，进行简要说明。在空间分析模块：① 探索哪些元素土壤异常与已知矿床/点有关，经矿点.wt与土壤.wp点对面相交分析或面对点相交分析，对派生的点或面文件进行属性分析使土壤Cu、Zn异常与已知矿床/点具关联性；②信息图层提取，通过土壤.wp元素组合的条件检索使土壤Cu异常A21、土壤Zn异常A20区对区合并分析(或在输入编辑模块通过“工作区/添加区文件……”实现)派生A22；③Ni,j、Si,j统计及相应派生Ai、Di(i=20，21，22)，已知有矿单元点.wt·AND·Ai至Bi、基本单元点.wt·AND·Ai至Di，分别查看Bi、Di属性或属性统计趋于 Ni,j、Si,j；④ 找矿信息量计算，根据表1，由式(2)生成I20、I21、I22。采用类似方法，提取分散流异常、重砂异常及物探航磁异常的找矿信息量。

2.3 构造的找矿信息

试验区中，与铁、铜、铅锌多金属等成矿关系密切的构造以断裂构造为主，因此，通过信息量法探测其找矿意义。在MAPGIS空间分析模块：①矿点.wt对断裂构造.wl的点对线叠加分析到派生点文件的双属性分类统计(图2-a)；②根据断裂的产状特征，将断裂走向归纳为北东、北西、东西、北南向，通过走向从断裂构造.wl依次检出，分别与矿点.wt进行点对线叠加分析，派生的点文件依次进行双属性分类统计，其中，矿床/点与北南向断裂点线最近距离(图2-b)，分别有一定量的矿床/点与北东、北西、东西、北南向断裂的距离小于500 m，确定500 m为断裂缓冲(控/找矿)半径；③对断裂构造.wl派生的东西、北东、东西、北南向断裂数据分别进行缓冲区分析形成A15、A16、A17、A18，A16·AND·A17生成A19(北东、北西向断裂交会部位)；④Ai(i=1, …5)分别与已知有矿单元点.wt、基本单元点.wt分别进行点对面相关分析，查看派生的Bi、Di属性或进行属性统计形成Ni,j、Si,j。由式(2) 生成Ii(i=15, …19)。

图2 矿床/点与断裂最短距离频数直方图Fig.2 The shortest distance histogram of deposit / point and fracture

3 基于MAPGIS空间分析的成矿远景综合评价

3.1 基于信息量法的地质变量筛选

基于前述的找矿信息分析，主要基于Nij、Sij、Ii进一步筛选：①A3、A4、A5区内小规模零星分布，仅各识别有1单元存在，找信息量大(均为1.368 0)，若单独分别设置将显著影响信息量法变量筛选，因考虑到下石炭统-下二叠统是我省铅锌的重要贮矿层位而设置了A7。同时，大岭组或建瓯群也具一定找矿意义，且试验区已知的铁、铜、铅锌矿化位置与基底“天窗”空间关联密切，因此，宜构置为前震旦-下二叠统A9。②A20、A21与A22(A22=A20∪A21)重复，存在A20、A21的单元数相差较多，A21的少且同属关系密切的土壤异常类，采用合并的A22；③A23、A24与A25(A25=A23∪A24)重复，存在A23、A24的单元数相当，信息量不很大且分别为分散流、重砂异常，宜分别构置(A23、A24)；④A16、A17与A19(A19=A20∩A21)虽有重叠，但A19是基于断裂交会的找/控矿意义而构置的；⑤ 根据信息量法，直接删除与负信息量对应的变量A1、A11、A13。从而初选地质变量13项，在此基础，大于0的信息量Ii按从大到小排序后，分别按式③计算正信息的总和I+、式④计算有用信息△I+：

(3)

(4)

式中：i+为正信息从大至小重新排序的序号，序号越小对应的地质变量/标志对找矿的重要性越大；m+为正信息的地质变量/标志；k为有用/保留信息水平，一般取k=0.75，或根据经验给定。

根据k=0.75，I+=6.608 5、△I+=4.956 4，对比从大至小重新排序后的正信息，当i+=9时，其正信息累加为5.283 8>△I+，从而筛选A23、A24、A9、A17、A13、A26、A19、A18、A22等9个地质变量(未入选变量/标志仍具找矿意义)，构成成矿远景预测的评价指标/变量体系，对应的信息量 依次记为γ1, γ2, …γ9。

3.2 基于MAPGIS的地质变量赋值

3.3 综合信息成矿远景定量预测

3.3.1 信息量法成矿远景综合评价模型

则信息量法采用线性加权综合多变量：

(5)

式中：s为单元序号，s=1, 2, … S表示按成矿预测划分的全部单元；i=1, 2, … M为信息量法筛选的M个变量，XS×M为所有单元M个变量观测值矩阵，RM×1为变量信息量/权重矩阵/向量；xs,i为XS×M的s行i列元素，表示s单元的x(i)观测值；γi为RM×1的第i个元素，表示x(i)的信息量/权重。H则为综合信息量矩阵/向量，其中Hs为s单元的综合找矿信息量，0≤Hs≤△I+，越大越具找矿远景，H所有大于某临界值HO的单元为远景单元。

3.3.2 模型的实现

在MAPGIS 6.7属性库管理模块，将前述已赋A23、A24、A9、A17、A13、A26、A19、A18、A22值的基本单元点.wt的 x(1), x(2), …x(9)及“单元编号”等属性导出(导为DBF格式)得X280×9。根据导出的X280×9、式(4)(以Excel方式打开，将R9×1输入其中的一列，利用MMULT函数)得H280×1，并将H280×1连同x(1), x(2), …x(9)等导入基本单元.wp中(另存为DBF格式)。对远景单元“临界值HO”的确定，该文根据H280×1，采用以下方法：①异常法，以“平均值0.483 9+1.645×标准差0.687 6”确定临界值HO=1.615 0，则全区有远景单元22个(已知有矿单元9个)；②曲线法，H280×1按组距0.323 4划分为11组进行频数统计(图3)，频数曲线拐点大致位于信息量第5组(1.617 0)上(1.293 6)限间，若以其组中值为临界值，取HO≈1.455 3，则全区有远景单元25个(已知有矿单元10个)；③类比法，将12个已知有矿单元的综合信息量按从大到小排序，以漏矿可能性小于10%，即不漏矿的可靠系数(91.7%)确定临界值HO=1.289 6，则有远景单元34个(已知有矿单元11个)。因为异常法、曲线法的远景单元数量偏少，类比法的远景单元数量相对较合适，及作图法临界值HO为0.94～1.25(图4)，已知矿床/点除各有1处分别位于信息量小于0.94和1.25区域外，其余均位于信息量大于1.25的区域。因此，选择HO=1.20为临界值，则全区分别有成矿远景单元35个(沿燕山期侵入岩东、北侧周边基底窗口分布)、非成矿远景单元235个，35个远景单元中有11个已知有矿单元(图5)，与全区已知的12个有矿单元重复率93.7%，远景单元见矿率大于31.4%，而245个非远景单元中仅1个已知有矿(占0.4%)，即非远景单元的见矿率约0.4%，所以，定HO=1.20为临界值是可行的。

图3 信息量的单元频数曲线图Fig.3 Unit frequency curve of information content

图4 等信息量线示意图Fig.4 Schematic diagram of equal information

图5 长坑幅铁、铜、铅锌成矿远景预测示意图Fig.5 Iron, copper, lead and zinc mineralization amplitude prospect prediction diagram of Changkeng

3.3.3 预测结果

以HO=1.20为临界值，全区共有35个找矿远景单元。 因为信息量法仍是基于类比找矿的，是类比法的一种量化，而前述的找矿信息分析/提取是基于已知铁、铜、铅锌矿床/点的，所以，预测的35个远景单元为铁、铜、铅锌矿的找矿远景单元。分别以2.0977、1.615 0将成矿远景单元划分为A、B、C类，其中A类11个(6个有矿、有矿概率54.5%)、B类11个(3个有矿、有矿概率27.3%)、C类13个(2个有矿、有矿概率15.4%)，并圈定为12处成矿远景区。根据A、B、C类远景单元、已知有矿单元的分布，有利成矿条件、围岩蚀变及其它找矿标志的发育程度等，将圈定的12处远景区划分为A、B、C三类，其中，A类3处、B类3处、C类6处，以信息量法为切点，着重探讨基于MAPGIS空间分析的成矿远景定量预测方法、程序，圈定的铁、铜、铅锌成矿远景区12处及各种标志、变量的找矿信息量等成果，对该区进一步的地质勘查工作具指导借鉴意义。

本文是福建省地质调查研究院2011～2014年承担的“1∶5万福建漳平洛阳—安溪潘田地区矿产远景调查”项目等的进一步成果，使用了项目组建设的漳平洛阳—安溪潘田铁多金属成矿带多源地学数据库的数据，在此表示感谢。

1 肖克炎,张晓华,宋国耀.等.应用GIS技术研制矿产资源评价系统.地球科学:中国地质大学学报,1999,24(5).

2 胡光道,陈建国.金属矿产资源评价分析系统设计.地质科技情报,1998,17(1).

3 张小静,李佑国.基于证据加权法的攀西地区铜镍铂族元素找矿远景预测.物化探计算技术,2009,31(2).

4 赵江南,陈守余,商伟.基于GIS证据权重法的祁连山地区铜矿床成矿预测.地质找矿论丛,2009,24(1).

5 商伟,陈守余,赵江南.秦岭地区铅锌矿床模糊逻辑法找矿远景区优选.地质找矿论丛,2009,24(1).

6 徐文喜,张晖,张德会.矿化系数等值线法及其在矿产资源预测中的应用.地质通报,2007,26(12).

7 李堃,胡光道,段其发. 基于MORPAS平台特征分析法的成矿远景区预测——以个旧西区锡多金属矿为例.地质科技情报,2009,28(4).

8 Bonhanm-Carter G F, Agterberg F P, Wright D F. Weights of evidence modeling: A new approach to mapping mineral potential//Agterberg F P, Bonham-Carter G F. Statical App lications in the Earth Sciences.Geological Survey of Canada,1989.

9 Cheng Q. GeoData analysis system(GeoDAS) for mineral exploration: User's guide and exercise manual. Material for the Training Workshop on GeoDAS held at York University, Toronto, Canada,2000.

10 李随民,姚书振,周宗桂，等. 基于ARCVIEW证据权重法的成矿远景区预测——以陕西旬北铅锌矿富集区为例.地质找矿论丛,2007,22(3).

11 李随民,姚书振,韩玉丑，等.权重叠加模型在冀西北铅锌矿富集区资源预测中的应用.金属矿山,2008,(4).

12 李裕伟,赵精满,李晨阳.基于GMS、DSS和GIS的潜在矿产资源评价方法.北京:地质出版社,2007.

13 王永军,李名松,全旭东，等. 基于GIS的层次分析法在张家口北东地区金矿成矿预测中的应用.地质科技情报,2007,26(4).

14 何海洲,杨志强.找矿信息量法在广西大厂矿田新一轮成矿预测中的应用.矿产与地质,2007, 21(5).

15 赵鹏大,胡旺亮,李紫金，等.矿床统计预测(第二版).北京:地质出版社,1994.