0 引言

成矿元素是成矿系统中最基本的单元，利用其时空分布信息不仅有利于了解地质作用过程产生的物质分异和就位机制，而且可以直接圈定成矿物质富集区，厘定分布范围和边界，为找矿提供直接信息^[1-3]。因此，剖析成矿元素分布是研究矿床成因、矿化富集和赋存规律的有效途径之一。由于地质作用过程的长期性，矿化过程往往呈多期次性。成矿事件形成于特定的地质构造环境，不同的成矿期对应于不同成矿事件，而每一次成矿事件都会形成与其相对应的地球化学响应，导致复杂的成矿元素品位变化。从非线性和复杂性理论角度看，成矿过程是一种奇异性地质过程，可在较短的时间或空间范围内发生能量的巨量释放或物质的超常富集和堆积，表现为元素品位空间分布具有自相似嵌套结构，传统的时间序列和统计方法已经无法有效提取隐藏的非线性动力学的相关信息^[4-12]。

多重分形是非线性理论的一个重要分支，指在分形结构上由多个标度指数所组成的集合，其谱函数可以描述分形结构上不同的局域条件或分形结构在演化过程中不同层次所导致的特殊结构行为与特征^[13]。传统计算多重分形谱函数的方法是配分函数法，但该方法对受趋势影响的非平稳数据的计算结果会出现不正确的现象。多重分形去趋势移动均值(multifractal detrended moving average，MFDMA) 法是有效刻画非平稳序列分形结构的新方法，它是通过去除移动窗口的平均值趋势后进行拟合计算。该方法比配分函数法更为稳健^[14]，已在金融、环境和物理等领域有较好的应用^[15-19]。

斑岩型铜矿床中地质过程及其形成的地质现象是复杂的，在空间上的分布具有明显的不均一性，即便是各种条件相近的斑岩型铜矿床，其矿化的类型、强度、矿化分带特征也均有不同，甚至有些斑岩体中成矿元素品位很高却不能形成矿体。随着斑岩型多金属矿的发现和开采，以及现代分析测试技术的进步，斑岩型矿床成矿复杂性的相关研究越来越深入，对岩浆-成矿作用过程及含矿斑岩的动力学背景的探索也越来越重视^[20-22]。

本研究将运用MFDMA法，研究云南普朗斑岩型铜矿床4号勘探线多钻孔Cu元素品位序列的多重分形特征，并探索其对矿化强度的指示作用，为进一步建立区域成矿预测和资源潜力评价模型提供理论支持。

1 理论与方法 1.1 MFDMA算法步骤

令长度为N的序列为{x_i} (i=1, 2, …, N)，多重分形去趋势移动均值法的主要步骤如下^[23]。

第1步，计算累积和y(i)，构造新的序列：

第2步，在尺度为s的移动窗口中，计算第i个移动窗口的平均函数，即

式中：表示≤x的最大非负整数；表示≥x的最小非负整数；θ是在区间[0, 1]之间取值的位置参数，代表移动平均值在移动窗口中的位置。

第3步，计算残差序列ε(i)：

第4步，把得到的残差序列根据相等的长度s划分为N_s=int(N/s)(int(·)取整数)个不重叠的区间，再将序列反向分割一次，得到2N_s个区间段。定义第v个子区间的局部波动变量F_v(s)：

这里ε_v(i)=ε(l+i)，l=(v-1)s。

第5步，计算q (q为整数)阶全局波动函数F_q(s)：

当q=0时，

第6步, 根据不同的划分长度s，波动函数F_q(s)与s若具有幂律关系，则

式中，h(q)是广义Hurst指数，可通过最小二乘法拟合得到其估计值。h(q)与多重分形质量指数τ(q)相关，可表示为

根据Legendre变化，得到奇异指数α(q)以及谱函数f(α)满足：

多重分形谱α~f(α)形状呈钩或钟罩状，序列不同分布结构的多重分形谱曲线呈不同的形状，对应着不同的物理意义。多重分形谱也称奇异谱，其谱高差为

多重分形强度可以h(q)的极差表示:

多重分形谱宽度及其左右分支可定义为：

式中：α(0)是对应q=0的特殊值；α_max、α_min分别为α(q)的最大值和最小值；q_max、q_min分别为q的最大值和最小值；Δh、Δα及Δα_L和Δα_R分别反映总体及最大、最小测度子集多重分形的强度，值越大分布越不均匀；多重分形谱高差Δf(α)反映最大、最小测度子集数目间的比例，其物理意义是对分形结构上非均质性的量度^[13-14]。

1.2 算法检验与参数选择

在MFDMA法中θ是确定如何划分子区间位置的参数。利用de Wijs模型，对不同位置参数所估计的广义Hurst指数与其理论值的误差和谱图进行对比分析，探索适合地质数据的位置参数。de Wijs模型是H.J.de Wijs在20世纪50年代研究岩石中矿物分布时提出的，亦称为矿化富集模型。该模型蕴含着多重分形特征，已成为模拟成矿元素或矿物富集分布的经典模型^[24]。其一维迭代公式为

式中：n(k)表示数值k在二进制中取1的个数；N=n_max；d(0＜d＜1)为富集参数。此迭代过程不断重复，使分割的两个子集一个不断相对贫化，而另一个相对富集，d值越大指示相对贫化或富集的程度越高。

令h_th(q)表示de Wijs模型Hurst指数的理论值，则MFDMA法估计的广义Hurst指数h(q)与h_th(q)的误差记为

图 1是de Wijs模型不同位置参数的ΔH和多重谱分形曲线，其中d=0.3，θ为0、0.3、0.5、0.7、1.0，数据容量N=2¹²，尺度范围[5，N/4]，q为-5~5，步长0.5。从图 1a中可以看出，θ= 0时ΔH绝对值最小且最接近零，其次为θ=1.0，而θ为0.3，0.5，0.7时误差更大；图 1b中也是θ=0时的谱图最接近理论值，即θ=0时计算结果更可靠，称θ=0时为后向移动平均法(backward detrended moving average，BDMA)。下面的实例分析中我们均选择位置参数θ=0进行计算。

图 2是θ= 0，d为0.10、0.25、0.50、0.75、0.90时h(q)随q变化的曲线，以此检验模型MFDMA法对不同富集参数多重特征识别的有效性。可见不同富集参数de Wijs模型的h(q)对q存在依赖性，且随d的增加h(q)对q的依赖性增强，亦即富集参数增大使奇异性增大。这与地质中成矿元素或矿物富集意义一致，表明MFDMA法能有效识别不同富集参数de Wijs模型的多重分形特征。因此，该方法适用于成矿元素或矿物富集的地质数据。

2 实例分析 2.1 地质背景和数据来源

云南普朗矿床是中甸弧内近年发现的规模最大的斑岩型铜矿床，目前已探明普朗铜矿区内铜资源储量达650万t，属义敦-中甸印支期斑岩成矿带，构造上位于义敦岛弧带。其南北断续延伸500余km，介于西侧中咱-中甸陆块与东侧甘孜-理塘蛇绿构造混杂带之间，主要为晚三叠世火山-沉积岩系和印支晚期—燕山期花岗岩基。成矿作用发生于普朗复式斑岩体内，矿化与蚀变相伴生长，同时，由蚀变带中心向外铜矿化逐渐减弱，且斑岩体从中心向外，依次出现钾硅化带—靠外侧的绢英岩化带—外缘的青磐岩化带，矿化呈现出“面型”—“线型”—“稀疏细脉或星点状” ^[25-30]。

普朗铜矿床的成矿元素Cu品位数据源于云南省地质调查局，采用示波极谱法和碘氟法在云南省地矿局第三地质大队实验室完成测定。本文选取普朗铜矿床4号勘探线上8个钻孔(ZK402、ZK403、ZK405、ZK406、ZK407、ZK411、ZK415、ZK419)的Cu元素为分析对象，采样位置见表 1。ZK402、ZK403、ZK405、ZK407及ZK411位于钾化带中，ZK406与ZK415位于钾化带与绢英岩化带的过渡部位，ZK419主体位于绢英岩化带与青磐岩化带的过渡部位，并且在钻孔的深部揭露了角岩(图 3)。不同蚀变带中，Cu元素品位的分布存在一定的差异性，Cu的边界品位为0.4%，根据元素品位大于边界品位的比例，将钻孔矿化强度分布由弱到强分为3个等级：(Ⅰ)弱矿化区元素品位大于边界品位的比例小于20%，且矿体难以圈定；(Ⅱ)中等矿化区元素品位大于边界品位的比例介于20%~50%之间，矿化区域间断出现；(Ⅲ)强矿化区元素品位大于边界品位的比例大于50%，矿体部分连续。据此4号勘探线上的8个钻孔可分为：ZK406、ZK411、ZK415、ZK419属于弱矿化孔(记为等级Ⅰ)，ZK407属于中等矿化孔(记为等级Ⅱ)，ZK402、ZK403及ZK405属于强矿化孔(记为等级Ⅲ)。

① 云南地质调查局. 云南省迪庆普朗铜矿区勘探地质报告. 昆明：云南地质调查局，2009.

2.2 数据统计分析与检验

采用Jarque Beratest (JB test) 和Augmented Dickey-Fuller (ADF) 检验研究区Cu元素钻孔数据的正态性和平稳性，基本统计与检验结果见表 1。在5%的显著水平下，JB统计量均大于对应的临界值5.82，表明Cu元素品位分布拒绝了零假设。图 4是样本分位数与标准正态分位数对比双分位数图，可以看出散点两端都偏离虚实线，结合JB检验和双分位数图可得出云南普朗铜矿床4号勘探线钻孔品位均呈非正态分布。此外，通过平稳性检验，钻孔ZK402、ZK403、ZK405的ADF统计量值均大于1%显著性的临界值-3.49，表明ZK402、ZK403、ZK405品位为非平稳序列，其余钻孔Cu品位的ADF统计量小于临界值，即为平稳序列。因此，基于正态分布假设和平稳序列的传统统计方法就难以准确描述其非线性的复杂特征。

2.3 多重分形参数与矿化强度指示

根据公式(1)-(6)计算波动函数F_q(s)，其中q=-5~5，步长为0.5。图 5分别是q=-5，0，5时，ZK403、ZK406、ZK407对应的F_q(s)与s的双对数函数图，可以看出ln F_q(s)与ln s具有明显线性相关性，且相关程度较高，其他钻孔也有类似结果。由公式(7)，利用最小二乘法估计其斜率，即可得到广义Hurst指数h(q)。

由公式(8)-(13)计算相应的多重分形参数，其结果见表 2。其中：Δh和Δα是反映多重分形强度的指标，其分形参数值越大可表示多重分形强度越大；Δα_L对应q < 0，Δα_R对应q>0的情况；Δf的大小和符号反映了谱函数的高差及左、右侧开口的分布模式，可指示元素高品位值部分比低品位值部分在不同级次子区域上的平均密度变化差异；W反映多重谱对称性，当W=1时，多重谱图呈钟罩状。可以看到，Δh和Δα均不为零，表明所有钻孔Cu元素品位序列存在多重分形性，但其强度的大小有差异，且W=Δα_R/Δα_L≠1，表示多重谱图为非对称状；强矿化钻孔ZK402、ZK403、ZK405的对应值Δα_L < Δα_R，W>1，Δf>0，而中弱矿钻孔ZK406、ZK407、ZK411、ZK415、ZK419的Δα_L>Δα_R，且W < 1，Δf < 0，表明强矿化钻孔的高品位局部聚集，随着矿化程度的减小多重分形谱图由左钩转变为右钩，即地球化学模式以富集为主转向贫化；a(0)对应奇异谱的峰值，值越大代表分割子集越少，强矿化和弱矿化的a(0)值大于中等矿化，表明中等矿化钻孔元素序列的分布非均匀性大于强矿化和弱矿化钻孔。

图 6是4号勘探线剖面上按钻孔位置从左到右(ZK419—ZK406)对应Δh和Δα的变化图，从钾化带—绢英岩化带，Δh和Δα整体呈一致的波浪变化趋势，在其内侧(ZK403/ZK405)达到波谷，在外侧与角岩化/青磐岩化带又缓慢增大，ZK419最大，ZK406次之，反映在钾化带与绢英岩化带中，高品位分布相对集中，利于成矿。

2.4 多重分形来源识别

引起成矿元素品位序列的多重分形原因一般有两种：一种是由概率密度函数或厚尾产生的，称为分布多重分形；另一种是由序列的长相关性产生的多重分形，称为相关多重分形^[31]。重排前、后的广义Hurst指数分别记为h_orig(q)和h_shuf(q)。对这两种多重分形的区分，可采用数据随机重排前、后广义Hurst指数的变化来识别。如果属纯分布多重分形，则重排后序列的密度分布没有改变，即有h_orig(q) =h_shuf(q)；如果属纯相关多重分形，重排后序列的长相关性被毁坏，则有h_shuf(q)≈0.5，呈接近完全随机状态。

对8个钻孔Cu元素品位序列随机重排20次，图 7是重排前后的广义Hurst指数曲线对比分析，〈h_shuf(q)〉是20次重排后h_shuf(q)的均值，可以看出所有钻孔的h_orig(q)和〈h_shuf(q)〉均随q的变化单调减少。其中：ZK415的h_orig(q)与〈h_shuf(q)〉最为接近，ZK402、ZK403、ZK405的〈h_shuf(q)〉 < h_orig(q)，ZK406、ZK407、ZK411、ZK415、ZK419的〈h_shuf(q)〉与h_orig(q)呈现前段减小后段增大的剪刀叉状，但总体h_shuf(q)曲线变得更平缓，且多重谱宽度Δh_shuf(q)明显变小(表 2)；两条曲线没有完全重叠，表明Cu元素品位的多重分形是由分布概率密度和长相关性共同引起，其多重分形结构主要归于长相关性，其次是概率密度或厚尾。

3 结论

1) 通过矿化富集模型验证了MFDMA法是识别成矿元素多重分形结构的有效方法之一，且当位置参数θ= 0时，计算结果具有较高的可靠性，为进一步地质数据的应用提供了理论依据。

2) 研究发现云南普朗斑岩型铜矿4号勘探线多钻孔的Cu元素品位的分布具有明显多重分形特征，其多重强度和谱函数存在差异，多重分形强度随矿化强度等级减弱而增加。根据剖面上钻孔位置从左到右分析可知，在钾化带—绢英岩化带内侧钻孔多重分形强度低于外侧与角岩化/青磐岩化带，表明在其内侧强矿化钻孔Cu元素品位分布相对均匀聚集，而在外侧与角岩化/青磐岩化带中、弱矿化钻孔Cu元素高品位分布相对分散，显示多重分形特征对矿化强度有指示意义。

3) 利用数据重排技术分析普朗Cu元素品位的多重分形结构来源是由分布概率密度和相关性共同引起，其多重结构主要贡献归于长相关性，其次是分布概率密度或厚尾，这可能是由斑岩型铜矿的硫化、蚀变叠加等控制因素的作用所致。

对研究区成矿元素的多重分形特征及来源与矿化强度指示意义的讨论，有助于查明构造环境对斑岩矿床成矿元素分布特征的影响，揭示成矿元素在成矿系统形成和演化过程中离散与聚集的机制以及时空分布格局，对认识和理解超大型斑岩铜矿成矿体系的复杂现象和形成机理具有较重要科学意义。

致谢: 对提供基础数据的云南省地质调查局表示感谢！