基于GIS的矿床空间信息成矿预测模型应用

矿产资源定量预测是在对矿床与地质条件之间的关系进行分析的基础上实施的，借用地质单元将矿床与地质信息联系起来建立预测模型，得以对已知的和潜在的矿床及其相应的资源量做出预测及评价（朱裕生等，1997）。赵鹏大等（1999）认为“在传统的矿床统计分析方法中，须用样本的观测结果来描述总体特征和确定远景区”。因此定量预测首先应保证抽样的随机性，还要保证样品的代表性。为此，通常选择一定大小的网格将整个研究区划分为面积相等、形状相同的“单元”。“单元”犹如地质取样中的样品，用作统一预测和取值范围的基本单位，同时也是进行成矿远景计算、比较、评价的基本单位。所以单元的大小和形状对预测效果有很大的影响（王於天，1990；李新中等，1998；池顺都等，1998；陈石羡，1998；陈永清等，1999）。

本次研究根据空间信息单元化定位预测方法，采用特征分析方法开展云南老君山矿田及其外围地区成矿条件及找矿预测研究，对地质、地球物理、地球化学、遥感等多元信息进行有机综合，建立矿床空间信息模型，为找矿靶区优选提供依据。

2.5.2.1 信息统计单元的划分

目前，在国内外的成矿预测中应用最广的是规则网格单元划分法，基本思想是运用统计学分析原理，在一定比例尺条件下选择一定大小的网格将整个研究区划分为面积相等、形状相同的单元，用作统一观测和取值的基本单位，通过样本的观测来描述总体，并遵循抽样的随机性及样品的代表性原则。这里如何确定最佳的网格单元大小是关键问题，矿点空间分布统计模型与单元面积大小也有直接关系（张振飞等，1999；吴红星等，2002）。

目前之所以主要采用规则网格单元划分方法，是因为它能在统一观察和定量的前提下，把众多的地质变量所包含的矿产资源信息量最大限度地反映出来，这有利于矿与非矿地质特征的判断，并且给矿产预测的计算机网格化带来了方便，尤其是在GIS支持下，网格单元的划分及单元中信息的提取非常便利。通常对预测单元的划分应考虑的因素有：预测比例尺和精度要求；预测区地质条件复杂程度、矿点数及空间分布特征；研究区范围大小及保证统计分析所需的单元数；地质特征的空间变异性等。

从矿点分布的方差S与其均值X的比例变化可知，单元面积越大，单元矿点分布模型越接近负二项分布；单元面积越小，则越接近泊松分布。因此，单元大小反映了不同的抽样观测条件。条件不同，则会影响统计分析的结果，单元面积的大小目前尚无明确的划分准则，但常用的经验算法如下：

1）经验性最佳面积S＝2×预测区总面积/矿点总数。

2）单元大小能保证当矿点的分布为随机型时，落入单元内的期望矿点数等于或小于实际落入单元矿点数标准差的3倍，即：

动态成矿作用与找矿

式中：E为落入单元内矿点数的数学期望；δ为实际落入单元内的矿点数x的标准差；S为单元面积（%）；n为矿点总数。若δ/E＝1/3，则S＝9/（9＋n）/100。

3）相应比例尺单元大小的参考数据区间，根据相应比例尺的地质图用1～4km²的面积为基本单元的大小，如对于1：5万地质图单元大小则0.25～1km²的面积比较适宜。

4）智能单元面积。根据算法1）、2）、3）中的因素和经验公式，系统将以上专家知识经验形式化、具体化，采用对话框提问方式，通过对预测区基本地质特征的询问，如预测区长度范围、地质图面上的矿点数目、比例尺信息，推理计算出最佳单元面积，然后转换成相应的网格图形叠加于地质图上。

5）用户单元面积。它主要根据用户直接提供的单元大小参数，进行屏幕图形的网格单元确定。单元大小及网格数的多少根据用户的经验和知识随意缩放、旋转，直到用户满意为止。

本次研究统计单元的划分主要考虑对矿化的显示，同时又考虑了统计计算、地质信息变量的选取和空间分析等要素。根据研究区的实际情况和统计计算的处理能力，采用规则网格法在1：5万的都龙老君山矿田地质图上按1km×1km的网格将研究区划分为1680个信息统计单元（图2.23）。

2.5.2.2 预测区地质信息变量的确定及编码

（1）地质信息变量的确定

根据2.1节的分析总结，可知老君山矿田的形成是老君山花岗岩体、田蓬组（

A.构造信息变量

地质地球物理资料显示，区域性的北东向文山-麻栗坡大断裂、马关-都龙大断裂等为继承基底构造，并具有间歇性活动演化的深大断裂，它们的长期发展演化及地壳的拉张作用造成了古断拉谷的形成，使该区经历了早期的强烈拉张下陷、晚期挤压隆起的发展过程，成为燕山期花岗岩浆侵位通道和就位空间。

在老君山矿田内，外接触带主要储矿体构造为纵向断裂及裂隙带组合，以及缓倾褶皱带内层间剥离、破碎带、裂隙带组合，矿体总体呈南北向展布；内接触带控矿构造为花岗岩边缘及内部的东西向、南北向、北东向裂隙带。因此，地表构造行迹可作为判断深部隐伏岩体形态变化及成矿有利程度的重要依据，尤其断裂密集区、断裂交会部、构造转折部及褶皱、断裂构造的复合部是成矿的有利部位。

B.地层信息变量

本区主要矽卡岩型矿床均赋存于田蓬组（

C.燕山期花岗岩信息变量

围绕老君山岩体突起周边分布有多处成矿带，总体上构成一环状带，并且这些成矿带延伸稳定，而老君山边部派生出来的小岩体或次级突起往往与次级褶皱和断裂构造有关，对成矿十分有利，因此，老君山岩体主突起周边及次级突起部位是寻找隐伏矿体的重要标志。

D.矿化蚀变信息变量

花岗岩，尤其是矿点附近的花岗岩中普遍发育云英岩化、电气石化、萤石化、黄铁矿化蚀变，而外接触带围岩中则发育矽卡岩化、绿泥石化、硅化、褐铁矿化等蚀变，由于断裂系统的贯通，矿化和蚀变的范围会远远超出岩体所在的位置，其影响范围甚至可达地表附近。

图2.23 老君山矿田信息统计预测单元划分

1.硅化蚀变；2.遥感铁化蚀变；3.遥感泥化蚀变；4.褶皱构造；5.环状构造；6.地层线及符号；7.断层；8.矿点；9.矽卡岩

因此，地表矿化蚀变现象及其类型、强度和规模，对隐伏矿床的预测同样有效。

E.遥感异常信息变量

由于老君山矿田岩石出露程度较高，对遥感蚀变信息的提取十分有利，加之蚀变岩石与广泛出露的碳酸盐岩地层及未蚀变花岗岩在光谱特征上反差明显，因此所提取的硅化、铁化和泥化蚀变区域与现有矿床分布区具有高度的吻合性，可作为隐伏矿体成矿预测的有效标志。

综上所述，控制和影响老君山矿田成矿的地质信息变量非常复杂，但在提取与找矿有关的信息时，有些信息是定性而不是定量的，因而无法进行统计计算。在统计变量选择的过程中，既要考虑选择的变量便于计算机进行空间分析与计算，又要注重选择有利于成矿预测的代表性变量。因此，本次研究为了建立空间定位预测模型，为了使参与叠加分析的所有图层都包含在研究的预测范围内，特选择以下变量进行空间分析，其中包括岩体、矽卡岩、绿片岩、北东向断裂、北西向断裂、东西向断裂、南北向断裂、构造交会部位、岩体与围岩接触带、遥感硅化蚀变、遥感铁化蚀变、遥感泥化蚀变、遥感构造解译、

（2）地质信息变量编码设计

为了有效地组织和管理上述地质信息变量，需要依据变量实体之间不同的特征、相似的特征以及不同变量实体的组合特征来对地质变量进行编码。对地质信息变量的编码设计是在分类的基础上进行的。我们选择的17个变量基本上可以分为五类：地层、构造、岩体、蚀变、矿点。在编码过程中要注意对整个系统的数据进行系统设计、统筹安排，使系统数据编码具有较强的系统性。综合考虑以上原则，结合地学空间数据的特点，参考有关国家标准，本书编码体系见表2.13～表2.18。

表2.13地层信息编码表

表2.14 断裂信息编码表

表2.15 岩体信息编码表

表2.16 蚀变信息编码表

表2.17 褶皱信息编码表

表2.18 矿点信息编码表

2.5.2.3地质信息变量赋值

地质变量在GIS中是以层的形式存储在数据库中的，因此针对不同的地质变量其属性也是不同的，但基本上可以分为三类：点文件、线文件和区文件。地质变量取值的实质是统计各网格单元内是否有点分布在网格内；是否有线通过网格；是否有某个层位的区文件覆盖网格。处理这样的变量在以往的研究中通用的取值方法是二态赋值法，即如果该地质变量在某一划分网格中存在，则其值为1，否则值为0。在这里，我们采用MapGIS软件中的空间分析模块对各个地质变量进行叠加分析取值。对不同的地质变量，其空间分析方法亦不同（池顺都，1998；刘春学等，2003），具体如下：

1）点变量取值：判断某个网格单元内是否有点分布，如矿点等。在MapGIS空间分析模块中空间分析菜单下，用区空间分析的区对点相交分析，就可得到含矿单元的区文件。这个区文件中就包含了所有取值为1的预测矿点。

2）线变量取值：判断某个网格单元内是否有线通过，如断裂等。在MapGIS空间分析模块中空间分析菜单下，用区空间分析的区对线相交分析，就可得到有断裂通过单元的区文件。

3）面变量取值：判断某个网格单元内是否有面通过，如地层。首先用空间分析模块的条件检索功能，根据地层代号（若无此属性字段，可在编辑模块中根据地层颜色参数统改层号、改当前层、存当前层等功能）将地层分布图分解成几个区文件，每个文件只包含一个地层单位，有几种地层单位（或岩体）就分为几个区文件。判断某个单元内是否有某一地层出露，可用空间分析模块中的检索菜单下的区域内检索功能，在对话框中选择区域条件文件为地层区文件，被检索文件为网格单元区文件，就可生成有某一地层通过单元的区文件。有几个地层区文件就做几次区域内检索，并生成相应数量的区文件。

利用我们已划分好的网格进行地质变量取值，这样就形成m（1680）×n（17）数据矩阵，m表示网格数，n 表示变量数。作为参考，这里仅列出10个已知矿点的变量取值表（表2.19）。

表2.19 变量取值表

2.5.2.4 成矿有利度法的数学描述及其确定

（1）数学描述

成矿有利度法是希腊和德国地质学家和数学地质学家合作推出的，该方法在1986年意大利国际数学地质讨论会上受到了各国数学地质工作者的好评。其数学表达式为：

动态成矿作用与找矿

式中：f为成矿有利度；w_i为第i个找矿标志的权系数；c_i为第i个找矿标志；p（c_i）为第i个找矿标志出现的概率；n为参加估计的找矿标志个数。

从式（2.6）可以看出，在成矿有利度法的数学表达式中，各找矿标志的权系数的确定是建模的关键。

变量权系数w可根据下列矩阵方程求得：

（CC^T）w＝λw （2.7）

这里λ是（CC^T）的最大特征值，C是m×n矩阵，代表n个地质变量在m个网格单元上的取值，C^T是C的转置矩阵。

地质变量以二态赋值方式赋值，即预测单元内出现为1，否则为0；数值型变量则以实际数值归一化后赋值。地质变量型找矿预测标志出现的概率以统计方法估计，数值型找矿预测标志的概率以归一化数值替代。根据矩阵表2.20，应用10个已知矿点组成的数据矩阵，采用MATHCAD数学软件就可以计算出权系数w，代入成矿有利度公式。据此就可以确定找矿预测标志的权系数。然后将各找矿预测标志的权系数经正规化变换，使其和为1。由此可建立都龙老君山找矿预测数学模型。其主要过程包括：

1）建立地质变量矩阵C（表2.20），并求得转置矩阵C^T；

表2.20地质变量矩阵

2）根据地质变量矩阵C和转置矩阵C^T，加入中间变量B；B＝（CC^T），求得（CC^T）矩阵（表2.21）；

3）调用eigenval（）函数求得特征矩阵，再调用max（eigenval（），0）函数求得最大特征值λ；

4）最后求得对应最大特征值λ的特征向量eigenvec（B，λ），即权系数向量W。

表2.21地质变量（CC^T）矩阵

5）根据权系数向量W，结合成矿有利度公式（2.6）求得找矿预测数学模型表达式：

动态成矿作用与找矿

这样，我们就可以计算出每个网格单元中的成矿有利度。

（2）信息统计单元成矿有利度的确定

本次在都龙老君山研究区内共划分了1680个网格信息单元，将信息单元的成矿有利度值按0.1的值域划分为9个信息数据组，并分别统计每组数据的频率。根据表2.22绘制成矿有利度频率分布图（图2.24），从图上的频率分布点，可确定预测单元的找矿信息临界值为0.6。在全区的1680个单元中，有111个单元的成矿有利度≥0.6，其中有47个单元为已知有矿单元（赵鹏大等，1999；谢贵明等，2000；曹瑜等，2003）。

表2.22 成矿有利度分级表

对信息统计单元成矿有利度数据进行多重分形研究，采用本章前述的C-N分形模型，得出信息统计单元成矿有利度分形曲线（图2.25）。曲线连续性很好，从成矿有利度来说，表明矿田成矿也有很好的连续性，区内成矿条件优越。

2.5.2.5 矿床空间信息成矿预测模型的实现

预测成果输出有两种，数据输出和图形输出。数据输出是利用预测模型在数学软件MATHCAD中计算后，将已预测结果写入到属性数据表中，通过查找数据表即可得到。再利用生成的信息单元数据成图（图2.26）。

图2.24 成矿有利度频率分布图

图2.25 成矿有利度频率对数图

在综合分析成矿条件、成矿机制基础上，依据成矿规律和找矿标志圈定靶区（王雄军等，2008）。按成矿有利度0.5、0.6和0.7为异常分界点，将预测单元分为三级，即A级、B级和C级，其中A级预测单元（大于0.7）为成矿条件最有利，找矿标志明显，并具有寻找大型多金属矿床的潜力；B级预测单元（介于0.6和0.7之间）为成矿条件比较有利，找矿标志较明显，具有寻找中型多金属矿床的潜力；C级预测单元（介于0.5和0.6之间）为成矿条件较一般，但仍有成矿可能，具有寻找小型多金属矿床的潜力。

图2.26 老君山矿田找矿信息统计预测图

1.A级预测单元；2.B级预测单元；3.C级预测单元