本发明涉及矿产资源勘查领域，具体而言，涉及一种基于层次分析法的成矿预测方法。
背景技术：
：随着我国经济的快速发展，矿产资源的消耗日益剧增，经济发展与目前已开采矿产资源不足的矛盾日益明显。目前低海拔近地表的矿产资源日益减少，找矿逐渐由浅层找矿向深部找矿发展，找矿难度不断增加，都给找矿工作带来极大的挑战。目前，用于成矿预测的参考因素主要有遥感地质解译、遥感蚀变异常提取和地球化学异常提取。现有的成矿预测模型中，一种是仅参考了上述三大因素中的两种甚至一种；另一种是考虑了所有参考因素，但各参考因素的权重是通过模糊的主观判断敲定。所以现有的这两种成矿预测模型均使预测结果存在较大的误差。技术实现要素：本发明的目的在于提供一种基于层次分析法的成矿预测方法，所述方法不仅考虑了遥感地质解译、遥感蚀变异常提取和地球化学异常提取这三大因素，而且利用层次分析法更准确地确定各因素的权重，有效地减小了预测结果的误差。为了实现上述目的，本发明提供以下第一技术方案：一种基于层次分析法的成矿预测方法，所述方法包括以下步骤：基于遥感影像数据进行地质特征解译，得到地质特征解译结果图；基于遥感影像数据进行蚀变异常区圈定，得到蚀变异常区圈定结果图；基于已有的地质资料，通过地球化学水系沉积物测量，进行蚀变带内的地球化学异常区圈定，得到地球化学异常区圈定结果图；将地质特征解译结果、蚀变异常区圈定结果和地球化学异常区圈定结果作为成矿影响因子，利用层次分析法计算各成矿影响因子的权重；将地质特征解译结果图、蚀变异常区圈定结果图和地球化学异常区圈定结果图叠加，根据各成矿影响因子的权重，计算叠加图中各像素点的综合权重值，并根据所述综合权重值，圈定成矿预测区。基于本发明的第一技术方案，第一种实施方式为：在进行地质特征解译和蚀变异常区圈定前，所述第一技术方案还包括预处理所述遥感影像数据，具体包括：地质特征解译前，融合不同时段、不同空间分辨率的遥感影像数据；对融合后的遥感影像进行拉伸处理，增加图像的可视性；遥感蚀变异常区圈定前，统一遥感影像数据中各波段的空间分辨率，并将各波段合并；对统一空间分辨率后的各波段进行大气校正；对大气校正后的各波段进行干扰地物剔除。基于本发明的第一技术方案，和所述第一技术方案的第一种实施方式，第二种实施方式为：所述基于遥感影像数据进行地质特征解译，具体包括：在遥感影像上分别建立断裂、褶皱、环形构造和岩体接触带的解译标识，结合已有的地质图和地形图资料，通过目视解译对解译标识进行归纳整理，得到地质特征解译结果图。基于本发明的第一技术方案，和所述第一技术方案的第二种实施方式，第三种实施方式为：所述基于遥感影像数据进行蚀变异常区圈定，具体为：选用aster遥感影像数据作为所述遥感影像数据，通过主成分分析法进行蚀变异常区圈定，具体包括：(1)基于蚀变带内的主要蚀变矿物，从aster1至aster9波段中选取的四个波段及四个波段分别对应的各像素点，建立协方差矩阵；(2)求出所述协方差矩阵的特征值和特征向量；利用所述特征向量作为行，求出变换矩阵；(3)利用所述变换矩阵，再对每个像素点进行变换，求得图像的四个主分量pc1、pc2、pc3和pc4；所述每个主分量均为四维向量；(4)根据蚀变矿物的波谱特征，从四个主分量中选定与所述蚀变矿物相对应的主分量；(5)统计所述选定的主分量的像素灰度平均值和标准离差；(6)根据所述像素灰度的平均值和标准离差，确定蚀变矿物的一级、二级和三级蚀变异常区的阈值，根据所述阈值分别圈定蚀变矿物的一级、二级和三级蚀变异常区，得到蚀变异常区圈定结果图。所述第三种实施方式中，从aster1至aster9中选取波段的方法如下：铁染蚀变矿物选择aster1、aster2、aster3、aster4波段；铝羟基蚀变矿物选择aster1、aster3、aster4、aster6波段；镁羟基蚀变矿物选择aster1、aster3、aster4、aster8波段。所述第三种实施方式中，建立所述协方差矩阵的公式如下：协方差矩阵协方差矩阵元素其中xi(k,l)为像素点在i波段的值，xj(k,l)为像素点在j波段的值；为各像素点在i波段的平均值，为各像素点在j波段的平均值。所述第三种实施方式中，从四个主分量中选定与所述蚀变矿物相对应的主分量，其选择依据如下：铁染蚀变矿物所对应的主分量的四个向量元素中，向量元素一与向量元素三的符号为负；向量元素二与向量元素四的符号为正；铝羟基蚀变矿物所对应的主分量的四个向量元素中，向量元素三的符号为正，向量元素四的符号为负；镁羟基蚀变矿物所对应的主分量的四个向量元素中，向量元素三的符号为正，向量元素四的符号为负。基于本发明的第一技术方案，和所述第一技术方案的第三种实施方式，第四种实施方式为：在基于已有的地质资料，通过地球化学水系沉积物测量，进行蚀变带内的地球化学异常区圈定完成之后，还包括步骤：基于所述地球化学异常区圈定的结果，通过比较各地球化学异常区的异常信息显著程度，将多个地球化学异常区分为一级、二级和三级地球化学异常区。基于本发明的第一技术方案，和所述第一技术方案的第四种实施方式，第五种实施方式为：所述将地质特征解译结果、蚀变信息提取结果和地球化学异常圈定结果作为成矿影响因子，利用层次分析法计算各成矿影响因子的权重，具体包括：(1)建立成矿影响因子的层次结构模型；(2)根据所述层次结构模型建立判断矩阵p；(3)根据所述判断矩阵p计算各成矿影响因子的权重。基于本发明的第一技术方案，和所述第一技术方案的第五种实施方式，第六种实施方式为：所述建立成矿影响因子的层次结构模型，具体包括：(1)将地质特征解译结果、蚀变信息提取结果和地球化学异常圈定结果作为一级成矿影响因子；(2)将断裂、褶皱、环形构造、岩体接触带、一级蚀变异常区、二级蚀变异常区、三级蚀变异常区、一级地球化学异常区、二级地球化学异常区和三级地球化学异常区作为二级成矿影响因子；(3)将断裂、褶皱、环形构造和岩体接触带分别对应的0～0.3km缓冲区、0.3～0.6km缓冲区和0.6～1.0km缓冲区作为三级成矿影响因子。基于本发明的第一技术方案，和所述第一技术方案的第六种实施方式，第七种实施方式为：所述根据所述层次结构模型建立判断矩阵p，具体包括：其中p为一级、二级或三级成矿影响因子所对应一级、二级或三级判断矩阵；所述pij为判断矩阵p中的第i个成矿影响因子与第j个成矿影响因子的影响程度之比，所述pij根据已有的权重数据资料得出；其中i＝1、2...n，j＝1、2...n。基于本发明的第一技术方案，和所述第一技术方案的第七种实施方式，第八种实施方式为：所述根据所述判断矩阵p计算各成矿影响因子的权重，具体包括：(1)计算各级判断矩阵p的最大特征值和最大特征值对应的特征向量；(2)利用所述最大特征值进行判断矩阵p的一致性检验，若所述判断矩阵p未通过一致性检验，则重新建立所述判断矩阵；(3)若判断矩阵p通过一致性检验，则根据所述特征向量，为各成矿影响因子分配权重。基于本发明的第一技术方案，第九种可实施方式为：所述将地质特征解译结果图、蚀变异常区圈定结果图和地球化学异常区圈定结果图叠加，根据各成矿影响因子的权重，计算图中各点的综合权重值，并根据所述综合权重值，圈定成矿预测区，具体包括：(1)将叠加图划分为若干栅格，每个栅格中包括若干像素点；(2)计算每个栅格中各像素点的综合权重值，并计算每个栅格中若干像素点的综合权重值之和，所述综合权重值之和为每个栅格的总权重；(3)根据每个栅格的总权重的大小，圈定成矿预测区。本发明所提供的技术方案不仅全部考虑了遥感地质解译、遥感蚀变异常提取和地球化学异常提取这三大因素，而且还利用层次分析法更准确地为每种成矿影响因子确定权重，避免了模糊地判断。由于本发明所提供的方法考虑的成矿影响因子更全面，各成矿影响因子的权重计算更准确，所以成矿预测结果的误差更小。附图说明为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。图1所示为实施例提供的基于层次分析法的成矿预测方法的流程图。图2所示为实施例提供的地质特征解译结果图。图3所示为实施例提供的基于遥感影像数据进行蚀变异常区圈定的流程图。图4所示为实施例提供的蚀变异常区圈定结果图。图5所示为实施例提供的地球化学异常区圈定结果图。图6所示为实施例提供的利用层次分析法计算各成矿影响因子的权重的流程图。图7所示为实施例提供的成矿预测区圈定结果图。具体实施方式为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。考虑到目前现有的成矿预测模型中，考虑的成矿影响因素不全面，或成矿影响因素的计算偏于主观，导致成矿预测结果的误差较大。基于此，本发明提供的一种基于层次分析法的成矿预测方法，能够使成矿预测结果更加准确，减少误差。本发明运用所述方法对西藏洛扎地区进行成矿预测研究，基于所述研究，提供以下实施例。请参阅图1所示，实施例所述方法包括：s101：基于遥感影像数据进行地质特征解译，得到地质特征解译结果图。本实施例选用两幅spot-6影像(时相：2015-2-5、2015-11-26，分辨率为6m)，和一幅etm+影像(时相：2000-11-17，分辨率etm+1～7为30m，etm+8为15m)作为用于地质特征解译的遥感影像数据。在进行遥感地质特征解译前，需预处理所述遥感影像数据，具体为：融合不同时段、不同空间分辨率的遥感影像数据；对融合后的遥感影像进行拉伸处理，增加图像的可视性。该预处理方法为本领域技术人员所熟知的方法，此处不再详细描述。在预处理后的遥感影像上建立断裂、褶皱、环形构造和岩体接触带的解译标识，结合已有的1:25万地质图和1:5万地形图，通过目视解译对解译标识进行归纳整理，得到地质特征解译结果图。请参阅图2，共解译出26条断裂，褶皱共15处，环形构造2处，以及侵入岩体6处。图2中：1-中新世二长花岗岩；2-中新世花岗闪长岩；3-岩体界线；4-推测断层；5-断层；6-环形构造；7-背斜；8-向斜；9-地名。s102：基于遥感影像数据进行蚀变异常区圈定，得到蚀变异常区圈定结果图。本实施例选用三幅aster影像(时相均为2015-12-14)，在进行遥感蚀变异常区圈定前，需预处理所述遥感影像数据，具体为：(1)将可见/近红外aster1-aster3波段的分辨率重采样至30m后，与短波红外aster4-aster9波段合并，统一为30m的空间分辨率，本实施例借助envi平台的layerstacking指令实现该步骤。(2)借助envi平台的flaash大气校正模块对统一空间分辨率后的各波段进行大气校正。(3)对大气校正后的aster1-aster9波段的遥感影像进行干扰地物剔除，具体请参阅表1。表1各类干扰地物的去除方法干扰地物类型去除方法阴影aster9/aster1低端切割水体(aster1-aster3)/(aster1+aster3)低端切割植被(aster3-aster2)/(aster3+aster2)第四系aster1高端切割请参阅图3，对预处理后的aster遥感影像数据进行蚀变异常区圈定，具体包括：s201：基于蚀变带内的主要蚀变矿物，从aster1至aster9中选取的四个波段及四个波段分别对应的各像素点，建立协方差矩阵。结合美国地质勘探局提供的蚀变矿物的光谱曲线图，对于铁染蚀变类的黄铁绢云岩化，选取aster1、aster2、aster3和aster4这四个波段。黄铁绢云岩化对应的aster1波段和aster3波段处于吸收谷，黄铁绢云岩化对应的aster2波段和aster4波段处于反射峰。对于镁羟基蚀变类的绿泥石化，选取aster1、aster3、aster4和aster8这四个波段。绿泥石化对应的aster8波段处于吸收谷，绿泥石化对应的aster4波段处于反射峰。本实施例中，提供以下建立协方差矩阵的方法：协方差矩阵协方差矩阵元素其中xi(k,l)为像素点在i波段的值，xj(k,l)为像素点在j波段的值；为各像素点在i波段的平均值，为各像素点在j波段的平均值。s202：求出所述协方差矩阵的特征值和特征向量。协方差矩阵c的特征值为λ0，λ0使特征方程|c-λe|＝0成立。求出协方差矩阵c的特征值为λ1、λ2、λ3、λ4，要求λ1>λ2>λ3>λ4，相应的特征向量为μ1、μ2、μ3、μ4。利用所述特征向量作为行，求出变换矩阵a。s203：利用所述变换矩阵a，再对每个像素点进行变换，求得图像的四个主分量pc1、pc2、pc3和pc4；所述每个主分量均为四维向量。本实施例中，所述四个主分量pc1、pc2、pc3和pc4，请参阅表2、表3：表2黄铁绢云岩化的四个主分量aster1aster2aster3aster4pc10.5968190.6076520.5222630.042515pc20.6226000.024990-0.714687-0.317754pc3-0.0903910.396990-0.4244390.808752pc40.498002-0.6874110.1905650.493097表3绿泥石化的四个主分量aster1aster3aster4aster8pc1-0.719658-0.653981-0.060764-0.225368pc2-0.5162530.2792200.2275920.776995pc3-0.4632800.7021370.094988-0.532309pc4-0.0307690.037829-0.9672060.249269s204：根据蚀变矿物的波谱特征，从四个主分量中选定与所述蚀变矿物相对应的主分量。其选择依据如下：考虑到黄铁绢云岩化对应的aster1波段和aster3波段处于吸收谷，则aster1波段和aster3波段为负贡献，则向量元素一与向量元素三的符号为负；黄铁绢云岩化对应的aster2波段和aster4波段处于反射峰，则aster2波段和aster4波段为正贡献，则向量元素二与向量元素四的符号为正。考虑到绿泥石化对应的aster8波段处于吸收谷，则aster8波段为负贡献，则向量元素四的符号为负；绿泥石化对应的aster4波段处于反射峰，则aster4波段为正贡献，则向量元素三的符号为正。s205：统计所述选定的主分量的像素灰度平均值和标准离差。其中，像素灰度平均值的计算公式如下：其中xn为主分量图像中每个像素点的灰度值，为主分量图像的像素灰度平均值。本实施例中，黄铁绢云岩化对应的xn的计算结果为0；绿泥石化对应的xn的计算结果为0。其中，像素灰度的标准离差的计算公式如下：本实施例中，黄铁绢云岩化对应的σ的计算结果为58；绿泥石化对应的σ的计算结果为20。s206：根据所述像素灰度的平均值和标准离差，确定蚀变矿物的一级、二级和三级蚀变异常区的阈值，根据所述阈值分别圈定蚀变矿物的一级、二级和三级蚀变异常区。具体为依据相应地质资料和已知蚀变带内的实际情况，选择1.5σ、2.5σ、3.5σ将黄铁绢云岩化和绿泥石化蚀变异常割裂为一、二和三级，黄铁绢云岩化对应的阈值分别为87、145、203；绿泥石化对应的阈值分别为30、50、70。在mapgis6.7平台中将黄铁绢云岩化和绿泥石化的蚀变异常区叠加，将所述两种蚀变矿物的同级别蚀变异常区重叠的区域圈定出来，得到一级遥感蚀变异常区6个、二级遥感蚀变异常区7个、三级遥感蚀变异常区6个，请参阅图4所示。图4中：1-一级蚀变区；2-二级蚀变区；3-三级蚀变区；4-地名。s103：基于已有的地质资料，通过地球化学水系沉积物测量，进行蚀变带内的地球化学异常区圈定，得到地球化学异常区圈定结果图。本实施例收集西藏洛扎地区已有的化探资料，通过本领域技术人员所熟知的地球化学水系沉积物测量，并利用mapgis操作平台开展了以cu、au、pb、zn、ag、w、mo、sb、cr、ti、ni、co、u、cs为主元素的地球化学综合异常圈定，共圈定了26个综合异常区，请参阅图5所示。图5中：1-au为主的综合异常；2-cu为主的综合异常；3-u为主的综合异常；4-cr为主的综合异常；5-pb为主的综合异常；6-ni为主的综合异常；7-sb为主的综合异常；8-ti为主的综合异常；9-zn为主的综合异常；10-ag为主的综合异常；11-co为主的综合异常；12-mo为主的综合异常；13-w为主的综合异常；14-cs为主的综合异常；15-地名。其中，以cu、au、pb、zn、ag为主元素的综合异常信息显著，异常浓度分带明显，单元素异常套合好，异常可能由矿化引起，作为最可能成矿区，即一级地球化学异常区。w、mo、sb为主元素的综合异常信息较为显著，单元素异常套合较好，做为中等可能成矿区，即二级地球化学异常区。cr、ti、ni、co、u、cs为主元素的异常信息较为分散，做为可能成矿区，即三级地球化学异常区。s104：将地质特征解译结果、蚀变异常区圈定结果和地球化学异常区圈定结果作为成矿影响因子，利用层次分析法计算各成矿影响因子的权重。请参阅图6所示，利用层次分析法计算各成矿影响因子的权重，具体包括：s301：建立成矿影响因子的层次结构模型。本实施例中将地质特征解译结果(yd)、蚀变信息提取结果(ys)和地球化学异常圈定结果(hy)作为一级成矿影响因子；将断裂(dl)、褶皱(zz)、环形构造(hg)、岩体接触带(yt)、一级蚀变异常区(ys1)、二级蚀变异常区(ys2)、三级蚀变异常区(ys3)、一级地球化学异常区(hy1)、二级地球化学异常区(hy2)和三级地球化学异常区(hy3)作为二级成矿影响因子；将断裂(dl)、褶皱(zz)、环形构造(hg)和岩体接触带(yt)分别对应的0～0.3km缓冲区、0.3～0.6km缓冲区和0.6～1.0km缓冲区作为三级成矿影响因子。表4成矿影响因子的层次结构模型s302：根据所述层次结构模型建立判断矩阵p。具体包括：其中p为一级、二级或三级成矿影响因子所对应一级、二级或三级判断矩阵；所述pij为判断矩阵p中的第i个成矿影响因子与第j个成矿影响因子的影响程度之比，所述pij根据已有的权重数据资料和专家的评价意见得出。本实施例中，所述判断矩阵的结果，请参阅表5。表5各级判断矩阵ps303：根据所述判断矩阵p计算各成矿影响因子的权重。(1)计算各级判断矩阵p的最大特征值和最大特征值对应的特征向量。(2)利用所述最大特征值进行判断矩阵p的一致性检验，若所述判断矩阵p未通过一致性检验，则重新建立所述判断矩阵。本实施例中，所有判断矩阵的一致性比率cr均等于0，即cr<0.1，通过一致性检验。(3)判断矩阵p通过一致性检验后，根据所述特征向量，为各成矿影响因子分配权重。本实施例中各成矿影响因子的权重结果，请参阅表6所示。表6各成矿影响因子的权重为了更清楚地表达成矿影响因子权重的计算方法，下面以一级成矿影响因子为例，其权重计算如下：(1)一级判断矩阵其最大特征值λ为3.0000，最大特征值λ对应的特征向量为p＝[0.54700.57630.6072]t。(2)一致性指标经过查saaty表得知ri＝0.58，则ci/ri＝0<0.1，满足一致性检验。(3)根据所述特征向量中元素的比重分配权重，即为0.3161(yd)，0.3330(ys)，0.3509(hy)。s105：所述将地质特征解译结果图、蚀变异常区圈定结果图和地球化学异常区圈定结果图叠加，根据各成矿影响因子的权重，计算叠加图中各像素点的综合权重值，并根据所述综合权重值，圈定成矿预测区。具体包括：将叠加图划分为若干栅格，每个栅格中包括若干像素点；计算每个栅格中各像素点的综合权重值，并计算每个栅格中若干像素点的综合权重值之和，所述综合权重值之和为每个栅格的总权重；根据每个栅格的总权重的大小，圈定成矿预测区。本实施例中，计算出每个栅格的总权重的大小后，选取一个值为0.076447477的权重点，将每个栅格的总权重除以所述权重点，即在叠加图上得到每个栅格中所述权重点的个数和密度。根据所述权重点的分布密度，圈定成矿预测区，请参阅图7所示。图7中：1-地名；2-权重点；3-新发现矿化点；4-向斜；5-背斜；6-岩体接触带；7-断层；8-环形构造；9-一级地球化学异常区；10-二级地球化学异常区；11-三级地球化学异常区；12-一级遥感蚀变区；13-二级遥感蚀变区；14-三级遥感蚀变区；15-成矿远景区。ⅰ号成矿预测区权重点分布最多且最密集，遥感影像上断裂构造极其发育，且解译有环形构造，中新世二长花岗岩体侵位其中，遥感蚀变异常强烈，发育以pb为主元素伴有zn、cu、au等元素的地球化学异常。显示出该区岩体、构造、地球化学、遥感等找矿指示因子相互耦合，成矿条件优良，具有寻找铅、锌、铜等多金属矿床的潜力。ⅱ号成矿预测区权重点大量分布且密集，断裂构造分布较多，发育中新世二长花岗岩体和中新世花岗闪长岩体，遥感蚀变异常强烈，分布1个以pb为主元素、2个以u为主元素的地球化学异常，是寻找铅、铀、铋矿床的有利地区。ⅲ号成矿预测区权重点分布较多，遥感影像上褶皱、断裂等构造较发育，断裂以北东向为主，遥感蚀变较强烈，分布有以cu为主元素、以au为主元素的地球化学异常，具有寻找铜金矿床的潜力。ⅳ号成矿预测区权重点大量分布且较密集，遥感影像上褶皱、断裂等构造十分发育，褶皱构造轴向为北西向，断裂构造为北西向和近东西向，遥感蚀变异常强烈，分布有以au为主元素、以sb为主元素、以ni为主元素的地球化学异常，是寻找金矿的有利地区。ⅴ号成矿预测区权重点分布较多且密集，遥感影像上推测有北东向断裂存在，遥感蚀变异常十分强烈且集中，区内分布有以au为主元素、以ti为主元素的地球化学异常，周边还分布着以cu为主元素、以co为主元素的地球化学异常，具有较好的铜金成矿背景。ⅵ号成矿预测区权重点分布众多且较密集，遥感影像上断裂、褶皱等构造都十分发育且均为北西向，遥感蚀变强度一般，区内地球化学异常集中，分布以au、ag、cr、sb为主元素的地球化学异常，且周边还分布有以au、ag、mo为主元素的地球化学异常，是寻找金银矿床的有利地区。以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本
技术领域：
的技术人员，在本发明揭露的技术范围内，可轻易想到变化或替换，都应该涵盖在本发明的保护范围内。当前第1页12