王巧,黄顺佳,魏春慧,江启煜,陈文成
(1.泉州师范学院 资环学院,福建 泉州 362000;2.福建省闽东南地质大队,福建 泉州 362000)
德化-尤溪矿集区位于德化、尤溪、永泰3县交界处,处于华夏地块之南武夷晚古生代拗陷区的大田-龙岩拗陷带东侧,闽东火山断拗带中段的云山巨型环状火山构造南西部,北东向政和-大埔及福安-南靖深大断裂所夹的寿宁-华安断隆带与北北东向浦城-尤溪断裂带交汇部位,区内构造活动频繁、岩浆作用发育,成矿条件优越,是福建省主要的金矿矿集区之一,已发现东洋[1]、下坂[2]、仙公寨[3]、双旗山[4]等几十处金矿(点).本区以往数十年的地质找矿工作,既有找矿突破,如发现火山-次火山岩型东洋金矿[1],也积累了丰富的找矿资料,但找矿难度也不断增大了,因此开展成矿预测评价,圈定金矿找矿靶区(远景区),不但条件充分,而且是迫切的需要,具重要意义.成矿预测评价长期是国内外地矿研究的重点、热点,形成了三部式、相似类比、地质异常等[5-9]理论,提出证据权法、特征分析法、逻辑信息法等[10-17]定量方法,尤其引入GIS[10-11、17-20],有利于充分挖掘地、物、化、遥等找矿信息,提高预测效率、水平.本文基于GIS空间分析,主要通过数据驱动建模[10-20],进行德化-尤溪矿集区金矿成矿远景信息量法[16-18]预测评价,包括空间数据准备、变量图层构置(含筛选、重构、赋值)、多变量(图层)综合等[10、17-19]内容,核心是基于GIS空间分析的数据驱动的找矿信息挖掘.
地质变量是地物化遥中各种成矿控矿因素、直接或间接矿化/找矿标志的泛称,因此多源信息融合的成矿远景定量预测评价必须准备的空间数据,包括(表1):(1)地物化遥空间数据,以1∶5万德化-尤溪矿集区(中仙、街面、水口、上涌、赤水、雷峰幅等)区域地质矿产图等为主要数据源采集的.(2)成矿预测评价单元划分及其相关数据,从地物化遥中定量提取地质变量的找矿信息,首先必须单元划分,而且事关定量预测评价的成败.地质变量源于点、线、面要素,其中,点、线作用于一定邻近区域,一般通过缓冲区分析转换为面状,曲面数据(如断层交点数、断层线密度、地层熵等)可通过DEM根据给定阈值进行“二值化”而转为面状或直接按单元统计,宜采用规则的单元[10、17、19].根据本区的控制程度、矿床/点空间格局、控矿构造特征及1∶5万预测要求,按平行/垂直主控矿构造的NE、NW向1000 m×1000 m划分网格,检出面积≥1000×1000×4/5=800000 m2多边形为基本单元G,取G的中心点派生g(点状),与金矿(表1)叠置派生已知有矿单元E(面状)、e(点状),并赋相同的ID、编号及“已知有矿[有已知金矿(化)点取1,无已知金矿(化)点取0]”等属性,总单元数S=1127、已知有矿单元数a=36.
表1 空间数据简表
注:属性仅指本次定量预测可能用到的.
在定量预测评价中,找矿信息挖掘是从分析各种地质变量在找矿预测中的作用、意义、贡献、重要性切入的,本文主要通过数据驱动的方式,应用信息量[16-18]模型进行定量分析.设Ai(i=1,2…n,下同)表示某一地质变量,是能提供找矿信息的地、物、化、遥等控矿因素、成矿条件、找矿标志等的统称.根据信息量法,地质变量Ai的找矿信息计量模型,如下:
Ii=lg(P(B|Ai,1)/P(B)),i=1,2…n
(1)
Ii=lg((P(Ai,1|B)/P(Ai,1))=lg((ai/a)/(si/s)),i=1,2…n
(2)
式中:①Ai表示第i项地质条件、因素、找矿标志,B表示有矿事件,1表示存在B.②s表示研究区划分的总单元数,即G的单元数,s=1127;a表示研究区中存在B的单元总数,即G中已知有矿的单元数,a=36;③si、ai分别表示研究区、已知有矿单元中存在Ai的单元数,即G或g、B或b中分别存在Ai的单元数,其关键是如何辨识单元是否存在Ai,一般采用[10、16、18]中心点法即若单元中心落入Ai中则存在、面积法即若单元内Ai的累积面积达到规定的临界值则存在,否则不存在,其中“临界值”为“0”时称重要性法,本文通过MapGIS 6.x空间分析进行识别、统计.④P(·)、P(·|·)为概率、条件概率符号.⑤Ii为Ai的找矿信息量,反映Ai的找矿意义、作用:若Ii=0则Ai不提供找矿信息;Ii>0,Ai存在有利找矿,且越大越有利;Ii<0,Ai存在不利找矿,一般作对立变换,即构置not·Ai,则not·Ai是有利的.
表1,面状要素包括地层、侵入岩、化探异常等,存在与否一般通过面积法或重要性法识别(表1X1~X8),根据成矿预测评价的“相似类比”原则,其找矿信息挖掘过程主要包括三步:第一步,获取已知金矿(或有矿单元)在哪些面要素中、已知有矿单元中有哪些面要素分布,MapGIS 6.x通过金矿(或b)、B分别与面要素相交叠置、属性分析等实现,如图1,这些要素Ai与金矿具空间关联性,可能提供找矿信息;第二步,利用式(2)分别计算这些要素的找矿信息量,通过条件检索依次得图层Ai即Ai.wp,Ai与G相交叠置(网格化得Ai·and·G.wp)、双属性(编号-面积累积)分类统计、面积约束检索(得到并存为Ai·and·G.wb)、属性分析得si、ai,代入式(2)得Ii,见表2.第三步,构置地质变量Xi,如表1所示(忽略了未入选图层模型的地质变量),与之相应,为Ai·and·G.wb新添属性Xi并统赋值1.
图1 点(已知金矿、单元)、面(地层、侵入岩、分散流异常)叠置派生数据的属性分析图表
表1,线要素有断裂构造及各种界面,点要素有断裂交汇点、破火口构造等,一般在邻域分析基础上,通过中心点法识别存在与否(表2X9-X14),进行找矿信息挖掘,主要环节包括:(A)初选点、线要素并获取相应的要素图层,如通过断层.wl与金矿.wt或有已知金矿单元b.wt叠置分析后的双属性统计(金矿与点线距离)得图2,可见NW向、NE向及环状断裂可能有利找矿,因此分别检索后得NW向.wl、NE向.wl、环状.wl,其中NW向.wl、NE.wl向分别进行10 m(图面0.2 mm)的缓冲区分析后相交叠置,提取相交多边形的中心点作为NW向、NE向断层交点(NW·and·NE.wt).此外,不整合、破火口也可能是有利因素,分别检索得不整合.wl、破火口.wl.(B)通过R-I曲线[17、19]探索点、线要素的合适的缓冲区半径,即求点、线要素的邻域、作用范围,就是在MAPGIS中通过循环Rt=Rt-1+△R(其中t=2,3…、R1=50 m、△R=50 m)的缓冲区分析得相应的缓冲区Ai-Rt.wp,再与g.wt相交叠置(网格化)得g·and·Ai-Rt.wt(中心点法识别)、属性分析得si,t、ai,t,代入式(2)得Ii,t并作R-I曲线(如图2-b所示),根据R-I曲线并结合si、ai等选择合适的半径,见图2、表2.(C)与合适的缓冲区半径对应的g·and·Ai-Rt.wt就是所构置的点、线要素类地质变量图层,新添属性Xi(X9-X14)并统赋值1.
表2 地质变量要素内容及其图层信息量统计表
图2 断层含矿性统计直方图及缓冲区半径与I值关系折线图
初选“破裂”度、断层线密度、断层条数密度等3个曲面数据变量,其中:(A)“破裂”度定义为网格单元内地层、岩体、岩脉等面状要素的“个”数,个数越多说明地质条件越复杂、构造越薄弱,越有利于找矿.(B)断层线密度为单元内所有断层总累积长度,值越大岩石越破碎、薄弱,越有利成矿.(C)断层条数密度为网格内所有断层的段/条数,值越大也反映岩石越破碎、薄弱,越有利成矿.对初选的3个变量,分别以750 m×750 m(图面15 mm×15 mm)网格采样后,经MAPGIS建立DEM后,通过尝试得最合适的阈值二值化(重分类)为“平面”、与g·wt相交叠置得g·and·(Ai-DEM).wt(中心点法识别)、进行属性分析得si,k、ai,k并代入式(2)得Ii,则与最合适阈值对应的g·and·(Ai-DEM).wt即为变量图层,新添属性Xi并统赋值1,见表2的X15-X17.
综合变量包括综合变量、组/复合变量、变量对立变换等,主要根据地质理论及变量的si、ai、Ii、专家意见等进行变量重构,以达到尽可能挖掘各种找矿信息的目的.如表2所示,构置的综合变量主要有X8、X14、X15、X16、X17即Au·and·(Ag·or·Cu)·not·(Mo·or·Sn)、NE与NW向断层交点、“破裂”度、断层线密度、断层条数密度等综合变量;X4、X3、X2即PA·or·δμ·or·γπ·or·λπ·or·λ、J3n1·or·J3n2、J3c·or·P2w·or·P2q2等组/复合变量;not·(Mo·or·Sn)等对立变换变量(未入选模型).and、or、not分别为MAPGIS中的图层相交、合并(添加)、相减叠置,si,k、ai,k统计、Ii计算及赋属性与前述各要素类似,结果见表2.
变量筛选是从构置的地质变量中选择对成矿定量预测评价具重要意义的变量,是多元成矿远景定量预测评价的主要环节,要求既不损失有直接、间接联系的主要信息又不存在冗余信息,达到简化与优化变量组合.信息量法选取正的、信息量较大的变量,按“从大到小提取全部正信息(Ii>0)累计的70%±”原则进行筛选[16],本文认为选取的Ii不限于正值,要求Ii绝对值适中且si不是很小并与ai有一定差值,若Ii绝对值较大而si却较小的一般作为构置组合变量,若Ii为负值进行对立变换.因此,在本区成矿背景、金矿实例、成矿规律等研究的基础上,依据si、ai、Ii等(表2),并兼顾不同变量之间及不同类型变量之间的信息均衡等,从与本区金矿找矿有关的地层、侵入岩、构造、分散流异常等调查成果中,筛选17项变量,则本区图层形式的金矿找矿预测模型及基于信息量法的数字找矿预测模型见表2.
3.2.1地质变量的量表与赋值
地质变量的量测、赋值是定量预测评价的基础,因为地质变量具模糊性、推断性、不精确性、不准确性等特征,因此采用基于变量存在与否的“0-1”二态赋值,即存在赋1、不存在则赋0.在找矿信息挖掘与变量构置过程中,已分别对Xi属性赋值,得相应的“.wt”或“.wb”并保存于指定位置,则“.wb”通过属性管理模块的“连接属性”工具、“.wt”在输入编辑模块应用“Label与区合并”工具,分别将表2各变量的同名属性Xi传递给G.wp,达到对G.wp的“0-1”赋值.全部变量对G.wp赋值后,通过“生成Label点文件”或与g.wt相交叠置,均可实现对g.wt赋值.
3.2.2属性的导出及其定量模型的数据集
G.wp赋值后,导出属性为“.xls”格式,其中变量X1-X17值集记为X,则:
X=[xj,i]s×n=[xj,i]1127×17
(3)
式中:(a)j=1,2……S,为预测研究区划分的网格单元的编号(取与ID、编号相同),S=1127,即共划分1127个单元.(b)i=1,2……n,为初选的变量的编号,n=17,即初选变量17个.(c)xj,i表示j单元i变量Xi的观测值,xj,i=1或0,=1表示j单元存在Xi,=0表示j单元不存在Xi.
3.3.1信息量法综合评价模型
根据有关文献,信息量法的变量综合模型如下:
(4)
式中:(a)i=1,2,…n为Ii>0并按Ii从大到小排列的新序号,即“Ii大于0”是变量筛选依据之一.(b)Ii即Ij,i,见式(2),单元j存在Xi时Ij,i=Ii、否则Ij,i=0.(c)相关文献依据有用信息水平k确定n值(k指按从大到小累积到全部正信息的70%±),舍弃了负信息及正信息小的变量,且有的正信息很大的变量的找矿信息未必真的大(s小并与a的差值很小的变量)并显著影响k值.(d)Ij为单元j的总信息量,表示该单元地质条件与标志对找矿的总有利程度,越大越有利.基于式(3),改进式(4)为:
[Ij]s×1=[xj,i]s×n·[Ii]n×1=[xj,i]1127×17·[Ii]17×1
(5)
图3 金矿成矿远景预测成果示意图
若要剔除某个变量Xk,则令Ik=0即可,方便应用Excel矩阵乘积函数“MMULT()”计算.
3.3.2总信息量计算结果
由表2及式(3)、式(5),得信息量法计算结果(保存为评价单元的属性),并在MapGIS的DTM模块经泛克立格法Kring网格化插值得.Grd格式的DTM(总信息量数字模型),经绘制等值线功能得视觉化预测评价结果,如图3所示.
3.3.3成矿远景划分与靶区圈定
本次采用曲线法、绘图法、“%”法等确定找矿远景下限临界值分别为1.190、<1.3296、1.2169,综合后取1.2169.同时,利用“标准差”法对成矿远景进行进一步的分类/级,设均值μ取样本的算术平均值,标准差s为基于给定样本的标准差,则
(6)
(7)
式中:Ij为j单元的总信息量(式5);若以总单元为样本则m=1127,并记均值为μ总、标准差为s总;以“切头切尾”后的单元为样本是检索符合“μ总+1.645×s总≤xj,l≤μ总-1.645×s总”的单元计数m,并记均值为μ切、标准差为s切.则A类/级临界值=μ总+1.96×s总=2.5799,B类/级临界值=μ切-1.96×s切=1.9625,C类/级临界值取前述综合的远景临界值1.2169,结果见图3.而且,在成矿远景分类/级的基础上,结合地质及构造、已发现的金矿、矿化蚀变等的空间分布特征,圈定找矿靶区(远景区)4处,如图3所示.
综上所述,本文:(A)广泛收集涉及德化-尤溪金矿矿集区的各种地质找矿资料,建立矿集区金矿成矿预测空间数据库,并划分NE-SW向1000 m×1000 m的网格为成矿远景预测评价的基本单元.(B)在领域专家支持下,通过数据驱动建模方法,利用信息量模型通过GIS空间分析进行智能化的找矿信息定量挖掘与变量筛选,构建本区金矿综合信息找矿图层模型、数字模型.(C)对信息量法进行改进,一方面是将“按有用信息水平k(k=70%)从大到小提取正信息量大的变量”的变量筛选方法改进为“综合考虑I、a、s及s与a差值等进行变量筛选与重构”;另一方面是改进了多变量综合的形式(即式5),与之相应,各变量的“I(信息量)、0”二态赋值统一改进为采用“1、0”二值化(是一种规范化值),有利于应用GIS空间分析进行智能化自动化的找矿信息挖掘及变量综合,而且方便变换为采用特征分析法、证据权法等其它定量方法.(D)对德化-尤溪金矿矿集区进行了1∶5万的定量预测评价,圈定金矿找矿靶区(远景区)4处及若干零星的远景区段,为该区进一步部署金矿找矿地质工作提供依据.