基于浓度-面积分形模型和模糊证据权的铜矿资源潜力评价：以江西九瑞地区为例

2024-01-22林俞亨王立立欧阳永棚李增华曾闰灵陈祺邓友国

西北地质 2024年1期

林俞亨，王立立，欧阳永棚,，*，李增华，曾闰灵，陈祺，邓友国

（1.东华理工大学地球科学学院，江西南昌 330013；2.江西省地质局第十地质大队，江西鹰潭 335001；3.鹰潭市紧缺及优势矿产勘查与研究重点实验室，江西鹰潭 335001）

铜具有优良的导热、导电性能和较强的延展性，因其优良的特性，被广泛应用于各个领域并在日常基础设施和技术中发挥着不可替代的作用（Ayres et al.，2003）。世界铜产量在过去70 年稳步增长，其中，中国是对世界铜产量做出重大贡献的国家之一，具有丰富的铜矿资源（Mudd et al.，2018）。此外，长江中下游成矿带是中国最具悠久历史的多金属成矿带之一（王兆强等，2023），位于该成矿带西部的江西省九瑞地区是中国重要的铜产区，该区域有较好的铜找矿潜力（徐耀明，2014）。另一方面，如地球化学测绘、证据权重法等矿床资源定量评价的方法，能够结合地质、地球物理、地球化学、遥感等数据进行信息综合，进行预测远景区的圈定，对指导找矿勘查具有重要意义（欧阳渊等，2023）。因此，可以利用已有的地球化学数据和地质数据等其他信息来进行处理，对九瑞地区指导铜矿产勘查（薛琮一等，2020；李文明等，2021；Li et al.，2022）。

地球化学微量元素作为成矿的直接物质来源，一直是地质工作者研究的重点，如何快速、高效的识别地球化学异常则是处理地球化学数据的一大难题。其中，分形/多重分形模型是识别地球化学异常的强有力的工具（Zuo et al.，2011，2012，2013）。自从Mandelbrot（1982）提出分形理论以来，已经有许多成果利用分形/多重分析理论从背景中提取出了异常。例如Cheng 等（1994）提出了浓度-面积（C-A）分形模型，并成功应用在米切尔-硫酸盐区金-铜异常的提取。各种分形方法，包括能谱-面积分形模型（S-A）、浓度-数量分形模型（C-N）、分形奇异值分解法（MSVD）等都能较好的用于不同矿化系统，并成功区分出地球化学异常和背景（Cheng，2004；Ali et al.，2007；Hassanpour et al.，2013）。

源自于没有空间意义的医学诊断上发展而来的证据权模型（Weight of Evidence，简称WofE）被国内外学者广泛应用于多元信息综合和空间决策支持系统，它可以整合多源地质信息作为诸多证据因子，从而形成潜在矿产分布的后验概率图（Agterberg，1989；Agterberg et al.，1994；Keller，1995；张津瑞等，2023）。Agterberg 等（1990）首次将证据权方法用于成矿预测中。在此基础上，Cheng 等（1999）发展出了模糊证据权重法。模糊证据权重（FWofE）是一种用于支持决策预测目的的数据集成方法,是一种空间决策支持方法。相比于普通证据权方法，模糊证据权重法通过引入模糊集合理论和隶属度函数，可以有效地减少由于遗漏或错误标记的数据引起的不确定性和提高地质数据的可信度，在确定模糊隶属度函数以代替后验概率计算的缺失数据方面更加灵活，为地质分析和决策提供更可靠的基础，其已被广泛应用于矿产潜力预测中（Cheng et al.，2002）。成秋明等（2007）采用模糊证据权等技术对老王寨及其邻区的金矿资源进行潜力评价，成功圈定有利成矿地段。Xiao 等（2020）综合普通证据权、模糊证据权和逻辑回归3 种方法，利用构造、岩性、地球物理和地球化学等数据，绘制新疆东天山造山带的斑岩型铜钼多金属矿床前景并进行定量评价。Zhang 等（2016）运用模糊证据权法得到预测区域仅占总研究区的5% 且成功识别了65% 的研究区内的铁矿床。Porwal 等（2006）开发的一种基于模糊证据权重方法和模糊逻辑方法相结合的子模型，该子模型在确定高矿产潜力区域方面是有效的，可以作为矿产勘探的工具。在基于多元信息和信息综合技术的矿产资源定量评价与成矿预测领域，模糊证据权重法是行之有效且应用广泛的方法（Cheng，2008a；黄秀等，2010；Zhang et al.，2014；陈风河等，2015；武进，2017；邓军等，2021；Huo et al.，2022；黄鑫怀等，2023；霍雨佳，2023；王佳营等，2023）。

江西九瑞地区位于长江中下游铜多金属成矿带中，该地区长期的构造作用、岩浆活动和成矿作用形成了丰富的铜等多金属矿床，具有勘探铜矿的潜力。但区域成矿地质条件复杂，致使传统的地质数据提取找矿信息难度较大。笔者基于九瑞地区1∶20 万水系沉积物地球化学采样数据、地质构造和控矿岩性等数据，在ArcGIS、GeoDAS 等软件平台上，综合利用因子分析、多重分形反距离权重插值、C-A 分形模型和模糊证据权等方法，构建九瑞地区模糊证据权模型，研究九瑞地区铜矿成矿潜力靶区，圈定可靠的成矿预测远景区，给予未来矿产勘探提供一定的参考意义。

1 地质背景

长江中下游成矿带位于扬子克拉通北缘，秦岭-大别造山带和华北克拉通以南。成矿带周边发育有3 条大型断裂，分别为东西向横贯成矿带南侧的阳新-常州断裂，西北缘呈NW 向的襄樊-广济断裂以及东北缘的郯庐断裂。研究区九瑞地区位于江西省北部，属长江中下游成矿带之九瑞铜金矿集区（图1）。

图1 长江中下游成矿带简易地质图（据Pan et al.，1999 修）Fig.1 Simplified geological map of the Mid-Lower Yangtze metallogenic belt

九瑞地区地层发育较齐全，以古远古界为结晶基底，中远古界为褶皱基底，晚震旦世以来的地层为盖层，由老到新为奥陶系至第四系，其中主要出露奥陶系至三叠系，研究区内矿床的形成、分布与地层关系密切，区内铜矿床的主要赋矿地层为泥盆系五通组、石炭系黄龙组和二叠系—三叠系大冶组（孔凡斌，2014；徐耀明，2014）。

区域内曾经历过强烈的岩浆活动，主要活动时期为燕山期，岩石类型以中酸性侵入岩为主（蒋少涌等，2013），产状多为小型岩株、岩脉和岩枝。研究区内岩体种类繁多，以花岗闪长斑岩和石英闪长玢岩为主，另分布有石英斑岩、花岗细晶岩、闪长岩、辉绿岩等。其中，燕山期花岗闪长斑岩是研究区内最主要的铜成矿母岩。

研究区内褶皱、断裂构造发育，构造系统分级明显，严格控制了岩浆岩空间展布与矿床的分布。区域构造主要为NW 向与NEE 向，其次为NNW 与NNE向，断裂及其构成的网格状构造节点，控制和影响着岩浆活动与成矿作用（图2）。

九瑞地区成矿岩体及其围岩普遍经受热液蚀变作用，矿床岩体多以岩体中心向外形成环状蚀变分带，与铜矿化相关的蚀变主要有矽卡岩化、硅化、绿泥石化。其中，据徐耀明(2014)制作的遥感蚀变光谱图可以得知，在空间上绿泥石化与已知矿床表现的较为密切。

矿体以层状、似层状、透镜状为主，铜矿主要类型有矽卡岩型、层控硫化物型、斑岩型等，区内截至2010 年累计探明的铜储量为0.3 万t（所颖萍，2013）。

在模糊证据权模型的构建的过程中，需要与成矿相关的地质要素转化为证据图层参与模型的建立。根据九瑞地区的地质背景，选择铜矿化的主要赋存地层、铜成矿相关的花岗闪长斑岩、影响铜矿空间分布区域断裂以及与铜矿空间分布的有关的绿泥石化蚀变作为证据图层参与模型的构建。这些地质要素与铜矿的分布和潜力密切相关，有助于模糊证据权模型更精确地评估矿产资源的分布和潜力。此外，通过对水系沉积物中提取的地球化学数据进行分析，可以识别出与铜矿化相关的地球化学异常，通常这些异常指示潜在的铜矿化区域。因此，地球化学异常也可以作为一类证据图层用于模型的建立，以进一步提高对矿产资源潜力的评估准确性。

2 研究方法

2.1 因子分析

因子分析是常用的多元统计分析方法之一，主要用于降维，并可以通过确定用于解释所观察到的相关矩阵所需的公共因子的适当数量和性质来计算出被测变量的更简约结构（Fabrigar et al.，1999）。同时也可以用于解决复杂的地质成因及矿化叠加问题，对叠加地球化学场也有较好的分解作用（董庆吉等，2008）。因子分析步骤包括确定因子载荷、计算变量共同度、因子旋转、方差贡献率和计算因子得分。其中，因子得分是观测变量的线性组合，它可以代替原来的变量用于后续的分析。

2.2 多重分形反距离权重法

在地球化学异常提取分析的过程中常有将点样本的属性转换为连续面的数值插值的处理过程。传统的插值方法（如普通克里金和反距离权重法）是空间插值分析中常用的方法。Cheng（2000）所提出的多重分形反距离权重插值（MIDW）在对地球化学数据处理在局部奇异性中表现出不错的优势。相比于传统的插值方法，多重分形方法能增强地球化学特征的预测能力，能更好的确定矿产勘探目标，其插值结果更能对弱异常进行精细识别（Parsa et al.，2017）。有关局部奇点指数和MIDW 的更多详细信息见Cheng（2007，2008b）和Zheng 等（2014）。MIDW 方法广泛应用于地球化学异常找矿中，可以增强地球化学异常区域，有助于地球化学勘探（Cheng，2004，2008b，2015）。

MIDW 的方程式可以表示为公式（1）：

式中：ω是多重分形权重；ε是单位单元大小的尺度（ε ＜1）；α是奇异性指数。

2.3 C-A 分形模型

C-A 模型是地球化学异常识别的强有力的工具，是人们常用的“地球化学异常建模的基本技术”（Carranza，2009）。浓度-面积分形法由Cheng 等（1994）首先提出，是用来分选地球化学异常的。浓度-面积分形法（C-A）是基于浓度值（量值的频率）的从背景中分离异常值的方法，这些量值也可以是地球化学指标的空间和几何特征。双对数坐标下的C-A 图可以用来建立浓度值大于s 的面积A（≥ s）和浓度值s 之间的幂律定律关系。从而可绘制出高、中、低异常区域地图。

C-A 模型表达方式见公式（2）：

式中：A(ρ) 表示浓度大于或等于轮廓值 ρ的区域；ν 是阈值；α1和 α2是大于零的分形维数。这2 个分形参数可以从 A(ρ) 对 ρ的对数-对数图中最佳拟合直线的斜率来估计。

2.4 模糊证据权重

FWofE 将证据层定义为具有多值隶属函数（0 ≤µ（A）≤ 1）而不是普通证据权重法中通常涉及的二元或三元集合。二元和三元模式成为模糊集的特例，例如µ（A）=1 或0 对应于二元模式，µ（A）=0、0.5 和1对应于三元模式。具有隶属函数1 和0 的模糊证据的2 个子集在FWofE 的实现中特别重要：A1={A :µ（A）=1}或A2={A : µ（A）=0}。第一组A1是证据的一部分，通常是勘探良好的区域，大部分矿床已被发现。第二组是发现的矿床很少，发现更多矿床的机会较小。假设与A1和A2相关的点事件（E）发生的条件概率为P[E|A1]、P[E|Ā1]、P[Ē|A1]、P[Ē|Ā1]、P[E|A2]、P[E|Ā2]、P[Ē|A2]和P[Ē|Ā2]。

E 与任何A 相关的模糊条件概率计算过程见公式（3）与公式（4）：

其中带“bar”的字母表示没有发生事件或没有证据。

模糊证据的权重计算过程见公式（5）：

式中：W+和W-是普通证据权重法中涉及的2 个普通权重；W0是先验权重。公式（5）具有以下特性：W[µ（A）]=W+，假设µ（A）=1 和W[µ（A）]=W-，假设µ（A）=0 和W[µ（A）]=W0，假设µ（A）=1/ 2 和A1+A2=Ω（整个研究区域）。

在条件独立于不同证据A 和B 的假设下关于事件E 的发生，给定A 和B 的E 的logit 计算过程见公式（6）：

以应用反演变换将logit 转换回后验概率P[E|µ（A）,µ（B）]，也可以计算后验概率的标准偏差，以测试结果的显著性。模糊证据权重方法的实施步骤概括为6 个步骤：ϕ 确定研究目标，如预测给定类型但未发现矿产的空间位置。κ 确定与目标相关的空间层（属性）,以便他们可以用于评估目标。λ 提取与目标相关的证据层，用隶属度函数表示证据层的可信度。µ计算模糊证据层的权重（衡量各个证据层对评估目标的重要性）。ν 综合多个模糊证据层，计算后验概率形成矿产资源潜力分布图。o 模型检验和预测区解释。

2.5 接受者操作特征曲线（ROC）

接收者操作特征曲线（Receiver Operating Characteristic Curve，简称ROC）自被创建以来已被用于各个学科，是一种可以用来评价分类器性能的可视化方法（Egan，1975；Zou et al.，2007；Swets，2014）。此外，ROC 也可以用于测试空间预测模型的有效性（Nykänen，2008；Chauhan et al.，2010；Rossi et al.，2010）。其主要原理是将异常空间分布模式与成矿有利要素或与已知矿床（点）进行比较，并通过相关统计指标来度量两者间的空间对应关系（周光锋，2021）。ROC 曲线通过反映真阳性（True Positive，TP，模型将阳性分为阳性）和假阳性（False Positive，FP，模型将阴性分类为阳性）以及真阴性（True Negative，TP，模型将阴性分类为阴性）和假阴性（False Negative，FN，模型将阳性分类为阴性）的权衡关系作为一种评估分类器性能的方法，ROC 曲线是纵轴上的灵敏度（真阳性率：TP/（TP+FN））与横轴上的1-特异性（假阳性率：FP/（FP +TN））的比较图，其曲线下面积（Area Under the Curve，简称AUC）可用作诊断测试准确性的度量，也可用于测量空间预测模型的性能。AUC 值可能在0 到1 之间变化，如果AUC 值等于0.5，则分类性能相当于完全随机猜测；而如果AUC 值等于1，则分类性能完美，即分类器可以正确分类所有样本（图3）。分类器的AUC 值通常介于0.5 和1 之间，其值越高代表的预测效果越好（Nykänen et al.，2015；Chen et al.，2016）。

图3 ROC 曲线示例图Fig.3 Example plot of ROC curve

3 数据与模糊证据权模型构建

3.1 数据来源

本研究用于建立预测模型的数据集有地质数据、化探异常数据和遥感解译数据，数据包括铜矿床空间分布位置数据、九瑞地区断裂数据、主要赋矿地层数据、花岗闪长斑岩分布数据、绿泥石化蚀变数据、水系沉积物化探异常空间分布数据，这些数据经过处理后作为证据图层输入模糊证据权模型中。其中，断裂、地层、花岗闪长斑岩的数据提取自1 : 25 万的九瑞地区构造图。绿泥石化蚀变数据来自徐耀明（2014）对九瑞地区遥感蚀变光谱图的解译结果。化探异常数据来自1 : 20 万水系沉积物地球化学数据的分析结果。

此外，1 : 20 万水系沉积物地球化学采样数据由江西省地质局第十地质大队按照采样密度每4 km2一个样品点位采集所得，且225 件样品点基本均匀覆盖整个研究区域（图4）。各微量元素检出限的测定包括：电感耦合等离子体光谱法测定Be、Co、Cr、Cu、Ba、K2O、La、Li、Mn、Nb、Ni、P、Sr、Th、V、Y、Ti、Zn、Zr 的浓度；掩蔽光谱法测定As、Bi、Cd、Sb 和Pb 的浓度；用比色法测定U 的浓度；催化极谱法测定Mo 和W 的浓度；Ag 和Au 的浓度通过原子吸收光谱法测定B、F 和Hg 的浓度分别通过发射光谱法、离子选择电极法和冷原子法测定。所有32 个元素的检测限制如表1 所示。

表1 32 种元素的检出限表Tab.1 Detection limits of 32 elements

图4 水系沉积物样品采样点分布图Fig.4 Steam sediment sample location map

3.2 数据预处理

对水系沉积物样品的32 个元素进行分析处理，用以提取九瑞地区的铜矿化异常分布图作为构建模糊证据权模型的证据图层。首先对标准化后的数据进行因子分析。因子分析方差最大正交旋转因子载荷表2 显示了各种因素的组成，其中F1 因子中的Cr、Co、Ni、Ti、V 为亲铁元素，主要存在于基性岩、超基性岩中，可能表明了铁的矿化。因子F2 含有Cu、Au、Bi、Mo、Pb、Sb、W 是研究区主要的成矿元素，相应的Cu、Au、Mo 元素属于矽卡岩-斑岩型铜钼金硫成矿亚系列，而Pb、Sb、Bi 为亲铜元素，W 则常产出于矽卡岩型铜矿中，因此F2 因子可以作为研究区主要铜成矿元素空间聚集的综合指标。因子F3-F7 与区域成矿无明显关系，对于成矿的指示性不强。

表2 R 型因子分析的正交旋转因子载荷矩阵表Tab.2 Orthometric rotating factor loading matrix for R-factor analysis

通过GeoDAS 软件，使用多重分形IDW 方法对F2 因子得分（Fac2）的空间分布图进行插值（图5），发现其现有的矿点与异常区域吻合程度较差，说明只进行插值不足以从背景中提取异常，需要进一步分析来提高异常区域提取的准确性。

图5 因子2 得分的多重分形反距离权重插值结果图Fig.5 Multifractal IDW result for the factor 2 scores

将浓度-面积（C-A）模型应用于因子F2 多重分形IDW 的结果进行地球化学异常与背景分离。为了研究元素含量与对应面积之间的分形特征值，通过GeoDAS 软件处理，做出对应的浓度-面积双对数散点图。在双对数图上背景值和异常值可以对应于不同的斜率值，不同线段所对应的分界值可作为区分背景和异常的临界值（图6）。经过分析，利用最小二乘拟合将图上的离散点拟合成4 条直线，从左往右分别表示为背景、弱异常、异常和高异常，最终得到研究区Cu 矿化异常分级图（图7）。

图6 浓度与面积的双对数图Fig.6 Log-log plot of concentration versus area

图7 C-A 分形模型识别的异常图Fig.7 Anomaly map identified by the C-A fractal model

此外，地质特征一定范围缓冲区与成矿是有关联的，并可以提供重要的勘探指南，所以可以对这些地质要素“延伸”一定的距离作为研究的证据图层（Zhang et al.，2016）。对于地层、断裂、花岗闪长斑岩和绿泥石化蚀变等地质要素，可以构建相应的缓冲区作为证据图层。根据已知矿床与各地质要素的空间关系，当主要赋矿地层（泥盆系五通组、石炭系黄龙组和二叠系—三叠系大冶组）的多环缓冲区为1 km、花岗闪长斑岩的多环缓冲区为0.5 km、断裂的多环缓冲区为1.3 km、绿泥石化蚀变的多环缓冲区为0.8 km 时，平均有85%以上的矿床落入各个地质要素的多环缓冲区内。在这些缓冲区范围外，随着距离的增加，所包含的矿床的数量增加极少。因此，这些距离可以视为最有利距离，即在这些缓冲距离内与已知矿床的空间相关性最大。经过数据预处理后，共得到了5 个建立预测模型所需的证据层：①花岗闪长斑岩多环缓冲区证据层。②断裂多环缓冲区。③主要赋矿地层多环缓冲区证据层。④地球化学异常证据层。⑤绿泥石化蚀变多环缓冲区（图8）。

3.3 模糊证据权模型

在构建模糊证据权模型前，需要在GeoDAS 软件对研究范围和相关参数进行定义，研究范围可以分为训练区域和预测区域进行定义，一般而言训练区域和预测区域为同一块区域。此外，为了确保每个栅格单元大小包含并仅包含一个矿床，设置研究区的面积单元大小为400 000 地图单元，选择11 个已知矿床作为训练图层，其先验概率为0.003 6。

在设置的训练参数基础上，对上述5 种证据图层分别进行模糊权重计算来构建模糊证据图层。普通证据权方法一般是以对比度C 与C 的标准偏差S（C）的比值t（分布t 值，t=C/ S（C））的极大值作为各证据层二值化的阈值，该阈值将各证据图层分成截然不同的二值图，即阈值之前的分类部分对成矿有利，赋值为1，阈值之后的分类部分对成矿不利，赋值为0。而对于模糊证据权，证据层由隶属度函数（MSF）进行重分类，隶属度函数是描述某个元素或时间对于某个特定集合（或模糊集合）的隶属程度或归属度，在地质中可以描述地质数据中的不确定性和模糊性并对地质数据进行度量，因此可以使用MSF 在闭区间[0,1]上任意取值对证据图层的确定性进行度量，根据已知矿点与不同分类的关系赋予不同MSF，一般的根据矿床在不同分类上的分布情况确定各类的MSF 取值，各证据图层的MSF 取值如图9 所示，其计算的模糊权重值如表3 所示。

表3 各证据层隶属度表（MSF）及模糊证据权重计算表Tab.3 Table of membership of each evidence layer (MSF) and calculation of fuzzy weights of evidence

图9 MSF 分类的证据图层Fig.9 MSF for evidences

4 结果与讨论

4.1 预测结果

每一个模糊证据图层刻画了与已知矿点的空间关系，构建模糊证据图层的过程是训练和学习的过程，综合5 个模糊证据图层创建后验概率图是达到预测和评价目标的必须过程。根据各证据图层的模糊证据权重，综合5 个证据图层，进行后验概率的计算，得到后验概率图（图10）。从后验概率图的结果来看，已知的11 个矿床的位置都具有较高的后验概率，城门山、武山、邓家山等矿床的矿体均位于后验概率高值区域。这表明后验概率高值区与已知铜矿床的分布存在显著的一致性，进一步验证了模糊证据权模型预测结果的准确性。高后验概率值的区域并不局限于已知矿床的分布范围内，表明在已知矿床范围内或周边可能具有成矿潜力。根据所得的后验概率图的后验概率值的大小及矿床的分布情况，将研究区划分3个一级预测远景区和4 个二级预测远景区，一级远景区与已知矿床在空间上的分布高度吻合。

图10 九瑞地区找矿后验概率及远景区分级图Fig.10 Map of posterior probability and hierarchical prospective in the Jiurui region

4.2 不确定性分析

经过邻域分析，得出距离成矿最有利的距离为7.5 km，在该距离内的铜矿床密度达到65%，并建立了相应的有利成矿区域（图11 白色区域），缓冲区以外的绿色区域代表负样本选择区域，并通过随机选点的方式，选取了与已知铜矿床相等数量的非矿点作为负样本。以研究区内已知的铜矿床作为正样本（True Positive，TP），在非矿点区域随机生成与已知铜矿床相等数量的非矿点作为负样本（True Negative，TN）。基于这两类正样本点和负样本点对所得的后验概率预测模型进行ROC 曲线分析，得到其ROC 曲线图以及AUC 值为0.983（图12），说明已知矿床与后验概率高值部分在空间上存在较强的相关性，表明了模糊证据权模型对九瑞地区铜矿床预测具有良好的预测性能。

图11 正、负样本点空间分布图Fig.11 Distribution map of positive and negative samples

图12 预测模型ROC 曲线Fig.12 ROC curve of the predictive model

此外，前人对九瑞地区矿床控制因素的研究表明，九瑞地区内铜矿除了受到区域性深大断裂的组成的菱形网格结点控制外，还受“四带一面”控制，即岩体与围岩接触带、层间破碎带、裂隙带与断裂带及不同岩性界面或沉积间断面控制，说明九瑞地区的岩体与断裂控制铜矿床的空间位置。将花岗闪长斑岩和断裂叠加在后验概率图上表明（图10），预测的远景区中花岗闪长斑岩发育，说明花岗闪长斑岩是研究区内最主要的铜成矿母岩，可以将预测远景区中花岗闪长斑岩发育区域纳入下一步的Cu 勘探中进一步调查。同时，圈定的二级远景区域表明即使有花岗闪长斑岩分布的区域，但是断裂发育较差，其后验概率会较低，说明九瑞地区铜成矿是受到花岗闪长斑岩和断裂共同影响的。模糊证据权重计算得到的高值后验概率区域和圈定的3 个一级远景区中，普遍有断裂和花岗闪长斑岩分布，表明基于浓度-面积分形和模糊证据权相结合的模型，具有一定的找矿指导意义。