ILR变换后数据的因子分区标准化在东昆仑东段地球化学异常圈定中的应用

2020-02-27耿国帅杨帆郭建娜

物探与化探 2020年1期

耿国帅，杨帆,郭建娜

(1.中国地质大学(北京) 地球科学与资源学院，北京 100083； 2.中国人民武装警察部队黄金地质研究所，河北廊坊 065000； 3.北京矿产地质研究院，北京 100012； 4.中国地质调查局土地质量地球化学调查评价研究中心，河北廊坊 065000； 5.河北省廊坊市自然资源和规划局，河北廊坊 065000)

0 引言

地球化学背景与异常划分是地球化学找矿中信息提取的关键问题。元素地球化学背景受地质环境、景观条件等诸多因素的制约，在一定区域内、同一景观条件下，地质环境是影响元素地球化学背景的主要因素[1-3]。以往按 1∶20万图幅确定统一异常下限的做法，导致低背景区矿化信息被掩盖，而高背景区出现了较多的非矿异常，背景与异常的划分存在一定的不合理性。为此，勘查地球化学工作者使用移动平均法、分区背景校正法[4]、子区中位数衬值滤波法[5-6]、归一化法[7]、小波分析法[8]、趋势面法[9-10]、 C型转换法[11]、泛克里格法和局部奇异性分析法[12]等来处理区域化探数据，在很大程度上提高了异常找矿信息的可靠性，有其优越性，但仍然存在一定的局限性。这些方法要么注重控制元素含量分布的地质因素，要么把地球化学场看作一个连续变化的曲面而采用各种方法去拟合。

时艳香提出了地球化学单元的概念，利用R型因子分析方法对样品进行地球化学分区[13]，并在实际应用中取得了一定的效果，使地球化学分区标准化得到了很好的应用和推广[14-19]。因子分区标准化的方法充分考虑了地质背景与元素组合之间的关系，有一定的客观性。但地球化学元素的含量数据是成分数据，具有闭合效应。Filzmoser指出，地球化学数据在进行计算时，如果不对数据进行变换，即使是做平均值和标准方差这样简单的统计或者是做直方图，都可能得出错误的结果；用原始数据统计标准差和方差几乎是没有意义的[20]，方差矩阵是做因子分析和主成分分析的基础，方差矩阵的不准确会使因子分析的结果与实际结果不相符，从而影响其效果。因此，笔者试图采用ILR变换后的因子分区标准化法来圈定地球化学异常，旨在进一步消除高背景地质体对异常圈定的影响。结果表明，与传统方法相比，利用ILR变换后的因子分区标准化方法所圈定的地球化学异常可以去除地质体中高背景值的影响，与矿床空间位置对应较好，且该方法受人为因素影响较小，实际应用效果显著。

1 研究区概况

青海省东昆仑东段跨越了古亚洲和特提斯两大构造域，位于南北两大构造域结合部位，区内可划分为秦祁昆、特提斯2个一级单元，东昆仑造山带、巴颜喀拉造山带2个二级构造单元及 4 个三级构造单元(表1，图1)。

表1 东昆仑东段大地构造分区

图1 东昆仑东段矿床分布及构造分区

不同的沉积环境具有不同的岩石组合，研究区昆北主要为海相沉积的碳酸盐岩和火山岩组合；昆中主要为古老变质岩和花岗岩组合，昆南主要为洋盆环境下形成的火山岩、碎屑岩和碳酸盐岩组合，而北巴则主要为三叠纪的复理石沉积建造。不同的岩石组合造成了不同的地球化学特征和矿床组合。

2 成分数据的概念和变换

成分数据是一组部分占整体的比例数据，只携带相对信息。目前成分数据不再被认为是具有定和的比值向量，而扩展为比值向量的等价类。一个成分数据是一个代表任意元素组合的等价类，其分量的和可以不是一个常数。元素含量作为相对含量，其和是一个常数，是典型的成分数据，各元素间不是独立的，而是相互影响的，这种影响对常量元素尤其严重，例如样品中SiO2含量的增加，会造成其他元素含量的总体降低，这就会造成SiO2与其他元素间通常呈负相关关系。因此在数据处理之前，如果不进行一定的变换，就可能得出错误的结论，尤其是相关性计算和多元统计分析[21-24]。

目前成分数据的处理与分析有3种变换：ALR变换(additive log-ratio transformation)、CLR变换(centered log-ratio transformation)和ILR变换(isometric log-ratio transformation)。3种变换的共同点在于变换形式都是变量的对数比值，因而统称为对数比变换。本质上，ALR变换和CLR变换不能消除成分数据的闭合效应，而ILR变换可以；原始成分数据经ALR变换和ILR变换之后，减少一个变量；CLR变换之后的数据与原始数据的变量数相同。

从3种变换的结果来看，ILR变换无疑是最好的，但由于变换结果减少了一个变量，增加了结果解释难度，而CLR变换能够保持变换前后变量数不变，因此Filzmoser利用标准正交基，把ILR变换后的因子分析结果在CLR空间下表示出来，既克服了成分数据的闭合效应，又突显了元素组合的地球化学意义。CLR变换和ILR变换以及与之联系的标准正交基简要介绍如下。

以一个原始数据矩阵Xij(m×n，其中m为样品数，n为元素个数)为例，分别介绍2种变换。

1) CLR变换

j=1,2,…,n；g(xi)为第i个样品的几何平均值。

2) ILR变换

ILR变换是在CLR变换和标准正交基vj的基础上得来的，标准正交基[25]

j=1,2,…,n-1。

因此CLR变换和ILR变换可以通过标准正交基联系起来。

clr(X)=V×ilr(X), ilr(X)=clr(X)×V,

V=[v1,v2,…,vn-1]。

笔者利用研究区所收集到的4 001件水系沉积物样品中11种常量元素(或氧化物)的含量数据，将其未经变换和经ILR变换后的数据作主成分分析(图2)。结果表明，未经ILR变换的元素(或氧化物)主要分布在主成分双标图的右半部分，除SiO2，其余10个元素或氧化物在第一主成分的载荷都是正数；在第二主成分中，除CaO、MgO、P、Mn等元素和氧化物的载荷为负外，其余元素和氧化物的载荷为正； SiO2和其他元素间隔较远，显示SiO2与其他元素的相关性较差，成独立一簇。

对数据作ILR变换，计算其主成分并变换在CLR 空间下显示其结果。经ILR变换的主成分双标图显示，11个元素(或氧化物)较均匀地分布在一～四号象限中，第一象限中是CaO和MgO，代表钙碱性元素的氧化物，第二象限Fe2O3、Mn、P、Ti等，代表的是基性元素组合，第三象限是Al2O3、K2O氧化物，代表碱性氧化物，而第四象限是SiO2、Na2O，代表的是酸性元素组合。

a—经ILR变换后的第一、第二主成分双标图；b—未经变换的第一、第二主成分双标图；图中红色圆点为第一和第二主成分的得分

同时，笔者又利用Ag、As、Au、Bi、Cd、Co、Cr、Cu、Hg、Mo、Ni、Pb、Sb、Sn、W、Zn等17种成矿元素经ILR变换和未经变换的数据分别做主成分分析双标图，结果显示，未经变换的数据双标图中，17种元素都分布在第一主成分的正载荷部分，而经变换后的数据，则较均匀地分布在各个象限(图3)。

综上对比结果可知，对原始化探数据作主成分分析有2个重大缺陷：①由于成分数据存在闭合效应，原始化探数据的第一主成分的载荷几乎都为正数，且这一现象普遍存在；②承载最大信息量的第一主成分通常是难以解释的，而经过ILR变换后的第一主成分则具有较为明确的地球化学意义。

a—经ILR变换后的第一、第二主成分双标图；b—未经变换的第一、第二主成分双标图；图中红色圆点为第一和第二主成分的得分

3 ILR变换的数据因子分区标准化

主量组分是构成岩石的主体，主量组分的组成和相对含量决定了岩石性质，支配了微量元素的地球化学行为及其在岩石中的分配，因此主量组分可作为地球化学分区的数据基础。

3.1 数据处理步骤

首先，选择研究区水系沉积物样品中常量元素的氧化物(Al2O3、CaO、Fe2O3、K2O、MgO、Na2O、SiO2)和含量相对较大的元素(Mn、P、Ti、Zr)共11种成分，并对数据进行ILR变换。其次，确定分区类型，即对样品中所选择的成分作因子分析，根据因子载荷矩阵中所反映的不同组合来确定分区的地球化学类型。然后，确定子区边界，因子得分是标准化值, 可以进行大小比较。当某样品在因子k中的得分值的绝对值最大, 说明该样品在元素组合k中所占份额最大, 故将其归为第k类(其中正值为k1类，负值为k2类)

根据分区结果对各变量进行分区标准化，圈定地球化学异常。标准化公式如下：

其中Zik表示标准化后的值；Xik为原始分析值；Q3k为第k类样品的上四分位值；IQRk为第k类样品的内散度，其值为上分四位点的值减去下四分位点的值。

分区标准化是基于EDA(exploratory data analysis)非常规化探数据处理方法，EDA技术以稳健统计学为基础，不需要设定假设条件，也无需对原始数据作任何处理，而是根据数据本身的特点来提取异常。EDA方法的统计参数包括中位数、2个极值点和上下四分位点，该方法认为数据集中有25%的数据是无效的，这些数据对整个数据集的中位数和上下四分位点不产生影响。绝大多数情况下，EDA方法提供的中位数比平均值能更好地估计数据的中心位置，内散度比标准离差能更好地估计数据的离散程度。它们都能抵抗特高值和特低值的干扰，比传统的方法能更有效地处理化探数据，因此EDA法受到了地球化学家的关注[26-32]。

根据EDA方法以Q3+1.5倍IQR作为地球化学异常下限，并根据上述标准化计算公式，就能很方便地圈定异常，即标准化值大于1、小于2是弱异常点，大于2是强异常点。

3.2 结果分析

基于主成分分析方法进行(初始)因子提取，首先对ILR变换后的11种成分数据进行主成分分析，计算出相关系数矩阵的特征根和解释率(见表2)；根据因子特征值(>1)和累计解释率(>75%)的要求，选择前4个主成分作为主因子，进行因子旋转(见表3)；然后计算各样品的主因子得分，制作各主因子得分计量图(图4)；根据各样品因子得分绝对值的大小，确定样品的地球化学类型，制作地球化学分区图(图5)。

由表3、图4和图5可知，F11为P、Ti、Zr、Fe2O3组合，反映的是与基性、超基性岩分布有关的元素组合，主要分布于昆南带内，此外还分布于昆南断裂带上，与该区中基性岩的分布较一致；F12为Na2O、CaO、SiO2、K2O组合，反映的是与中酸性岩有关的元素组合，主要分布于昆北及昆中带上，与中酸性花岗岩的分布范围较一致。F21为CaO、MgO组合，反映的是与生物沉积为主的灰岩、白云岩等元素组合，分布于昆中和昆南带上，与碳酸盐岩的分布范围较一致；F22为Al2O3、K2O组合，反应的是与碎屑沉积为主的砂泥岩的元素组合，分布于北巴颜喀拉，与北巴颜喀拉地层的板岩分布范围较一致。F31为K2O、MgO、Na2O组合，反映的是与河流的蒸发作用有关的元素组合，主体分布于柴达木盆地南缘；F32为SiO2、Zr组合，反映的是与河流的搬运作用有关的元素组合，主体分布于昆南断裂带以南。F31和F32是表生地球化学作用因子，在柴达木盆地的南缘，地表由于蒸发作用，可溶性组分结晶形成各种盐类，形成K2O、MgO、Na2O等易溶元素组合，东昆仑断裂带以南，蒸发作用弱，而与河流的搬运为主，形成难溶的Si、Zr等元素组合。F41为Na2O、K2O、Zr、CaO组合，反映的是与河流碎屑沉积为主的元素组合，主要分布于昆中带上，与河流碎屑沉积环境有关；F42为F2O3、Mn组合，反映的是与海相页岩、泥质岩为主的元素组合，主要分布于北巴，与海相的泥、页岩沉积环境有关。

表2 研究区主成分分析特征值和累积方差贡献率

表3 研究区正交旋转因子载荷

图4 研究区常量元素经ILR变换后4个因子得分计量

图5 研究区地球化学分区

为了解各地球化学分区中常量元素的含量分布，统计了8类地球化学分区中各常量组分的中位数(表4)。结果显示，P、Ti、Zr中位数的最大值都处于F11中，而最小值都处于F12中；CaO中位数的最大值处于F21中，而最小值处于F22中；Al2O3、K2O中位数的最大值都处于F22中；Na2O、MgO中位数的最大值及Al2O3、K2O中位数的次大值都处于F31中，而最小值都处于F32中；SiO2的最大值也处于F32中；Fe2O3、Mn的最大值都处于F42中，次小值都处于F41中；该结果与因子载荷所反映的一致(表3)。

根据本文分区标准化公式和EDA确定异常下限的方法，以1、2、4为异常的外、中、内带下限值圈定分区标准化异常图(图6a)；将剔除特高值后的原始数据平均值+2倍标准离差作为异常下限，以2、4倍异常下限为中、内带圈定传统异常分布图(图6b)。两种方法圈定的金异常点位分布对比如图7所示，异常点位特征对比列于表5。

表4 研究区各分区中常量组分中位数统计

a—ILR变换后因子分区标准化圈定的Au异常; b—传统化探数据方法圈定的Au异常

a—因子分区标准化圈定的金异常点及金矿床点；b—传统方法圈定的金异常点及金矿床点；c—两种方法都圈定出的金异常点及金矿床点；d—只在其中一种方法中圈定的金异常点

表5 两种方法圈定的金异常点位统计

从图6可知，用IIR变换后数据分区标准化圈定的异常和直接用原始数据圈定的异常在西部矿床(点)上都有良好的指示，而东北部指示效果不太理想，其原因推测与该区的表生地球化学环境有关，该区受风成砂影响严重。

直接用原始数据所圈定的异常较零星，浓集趋势相对较差，而用分区标准化所圈定的异常浓集趋势明显，与矿床点的对应关系也较好，因此用分区标准化圈定的异常更具有找矿指导意义。

结合图7和表5可知，用传统方法所圈定的金异常点数占总样品数的10.4%，而用分区标准化所圈定的异常点数只占总样品数的8.1%，两种方法都圈定出的异常点数占总样品数的7.7%，且金矿床(点)几乎都处于两种方法都圈定出的异常点位处，表明分区标准化圈定的异常具有更高的可信度。此外，传统方法还在图幅的东南端圈定了大量的非矿异常点。东南端处于北巴子区中，出露地层主要为三叠系巴颜喀拉群，其砂泥质类复理石沉积建造是Au、Sb、As、Hg的重要矿源层[33]。麻多幅和扎陵湖幅地球化学调查结果显示，巴颜喀拉山群变砂岩中Au的背景含量达到1.85×10-9[34]。将本文采用的4 001件水系沉积物测量数据按昆北、昆中、昆南和北巴4个子区分别统计背景值，结果显示，Au的背景值在昆南和北巴(昆南1.59×10-9，北巴1.54×10-9)远远高于昆北和昆中(昆北0.76×10-9，昆中0.96×10-9)，因此推断异常与该区的高背景有关。而分区标准化可以剔除高背景的影响，从而使异常圈定更合理。