1. 引言
矿产资源是国家经济可持续发展的重要物资基础和支持力量。地球化学异常识别是找矿工作的重要一环。在成矿作用的影响下,元素浓度的地表分布与成矿特征量相关,通过元素在地表的分布与地下空间矿体的相关性进行找矿预测 [1]。地质过程的复杂性使得地球化学元素含量的分布具有复杂的空间结构特征,准确有效地刻画地球化学元素分布模式对于提取地球化学异常也有积极的意义。因此,地球化学异常信息提取方法和模型研究是勘查地球化学领域的一个研究热点和前沿。
为了识别地球化学异常,以前的研究通过探索地球化学数据的频率分布或空间分布规律,在地球化学异常提取工作上取得了可喜的进展。单变量统计方法如标准差、概率图等是最简单直接的传统统计方法,在对单一元素进行概率分布的研究基础之上,研究元素的统计分布规律,确定地球化学背景和下限。多变量数据处理方法如主成分分析 [2]、聚类分析 [3]、判别分析 [4]、因子分析 [5] 等方法,以多元素或多元素组合为基本特征探索地球化学数据的规律,捕捉地球化学元素的频率分布特征并分离地球化学异常。
地球化学数据是一个典型的具有高维元素属性的地理空间数据。因此,许多方法通过分析地球化学数据中的基本空间结构和空间自相关的特点,在地球化学异常识别上取得较好的成果 [6] [7]。如克里格法、地理加权回归法 [8]、趋势面分析法等方法 [9]、局部差距统计学 [10]、局部邻域统计 [11]、空间因素分析 [12]、移动平均技术 [13] [14] 等方法被用来探索地球化学数据的空间结构特征。基于分形与多重分形理论的信息提取方法根据地球化学数据中的空间尺度不变量特性、几何特性和与空间面积有关的幂律关系来分离地球化学异常现象。如浓度–面积分形模型,频谱–面积模型(spectrum-area),局部奇异性分析(LSA) [15] [16] 等方法在地球化学识别上都获得了较好的效果。
最近,将机器学习方法引入到地球化学异常识别中也是一个研究热点。包括单类支持向量机 [17],连续限制波尔兹曼机 [18],孤立森林 [19],高斯混合模型 [20],深度自动编码器网络 [21],深度变异自动编码器 [22],GANomaly网络 [23] 等方法已被利用来学习地球化学元素的复杂模式和非线性关联。多卷积自动编码器 [24] 空间约束多自动编码器、空间约束多自动编码器 [25],深度卷积神经网络 [26] 堆叠卷积去噪自动编码器 [27] 等方法能够保留空间特征来分离地球化学异常现象。
上述地球化学异常识别方法都能有效提取异常,取得显著效果。然而,专注于从地球化学元素的频率分布特征或者元素相关特征中识别异常的这类方法,往往忽略了地球化学元素的空间结构,掩盖地球化学数据中一些与空间分布相关的潜在规律 [7]。而由于处理高维变量的局限性,关注地球化学元素空间分布的方法大多依赖于对高维地球化学数据的降维方法。这样的地球化学数据处理可能会导致地球化学特征的缺失。
为了综合考虑地球化学数据的元素相关性和空间分布规律,本文提出了一种基于多重主成分分析的方法识别地球化学异常,并选取山东省胶东半岛西北部为案例进行分析。使用多重主成分分析,将地球化学数据视为张量,并从张量的空间维度和元素维度上分别使用主成分分析方法提取信息,既考虑了地球化学数据的空间分布又顾及了元素之间的相关性。
常规的主成分分析方法在地球异常识别工作中一般作为多元统计方法出现,使用主成分分析提取地球化学元素组合信息,将多个元素变量降维,再针对单独的一个或者几个主成分,进一步使用数学方法或理论模型识别地球化学异常。作为多元统计方法的主成分分析仅在元素维度进行降维或者信息融合。而本文利用了主成分分析方法的信息提取能力,不仅对元素维度进行主成分提取,同时还对空间维度进行主成分信息提取,捕捉到地球化学数据的元素关联信息和空间分布信息。在主成分分析方法的转换下,大部分地球化学样本的信息都被保留下来,而异常代表的小部分信息未包含在主成分中。因此,通过多重主成分分析提取研究区水系沉积物地球化学数据的空间–元素联合信息,重建地球化学数据。并计算重建数据与原始数据之间的距离,提取地球化学异常。
2. 研究区地质概况
胶东地区做为我国最大的金矿集区,拥有200多个金矿点,已查明的黄金储量超过5000吨 [28] [29]。胶东地区的大地构造位置位于华北克拉通东南缘和大别–苏鲁超高压变质带东北端。研究区位于山东省胶东半岛西北部地区(见图1)。主要岩性是前寒武纪元古代火山岩和褐铁矿–闪长岩片麻岩及中生代花岗岩,局部覆盖有古新世火山岩、碎屑沉积物和第四纪沉积物。
区域内重要的控矿构造主要由三条控矿断裂自西向东构成:三山岛断裂、焦家断裂、招远–平度断裂。这几条北北东向和北东向断裂构成了胶西北地区最突出的线性构造,是与金矿密切相关的断裂构造。近90%的金矿资源都与北北东–北东向断裂有关,分布于主要断裂成矿带之间的区域。
区内金矿床的成矿条件十分复杂,具有明显的多期性、多因素控矿性。区域内大型–超大型金矿的形成主要受前寒武纪变质基底岩系、中生代燕山期构造–岩浆活动、北东、北北东向韧–脆性构造三大因素复合控制。主要金矿类型有焦家式的破碎带蚀变岩型金矿、玲珑式的含金石英脉型金矿和河西式的网脉型金矿。不同的金矿类型在伸展构造中出现的位置不同。其中,焦家式金矿分布于靠近主断面的黄铁绢英岩化碎裂岩带和黄铁绢英岩化花岗质碎裂岩带中,代表性矿床为焦家金矿、三山岛金矿和大尹格庄金矿。河西式金矿出现于断裂下盘焦家式金矿外围的黄铁绢英岩化花岗质碎裂带和黄铁绢英岩化花岗岩带中,玲珑式金矿则赋存于伸展构造下盘远离主构造带的黄铁绢英岩化花岗岩带和正常花岗岩内,比如玲珑金矿 [30]。
![](//html.hanspub.org/file/11-1771421x7_hanspub.png?20220602091954271)
Figure 1. Geological map of the Jiaodong peninsula (modified after Liu et al., 2021 [31])
图1. 胶西北地质特征图(改自Liu et al., 2021 [31]))
3. 多重主成分分析方法框架
我们将39种元素的地球化学数据视为一个整体张量。为了从背景中分离出地球化学异常,使用多重主成分分析,对地球化学元素数据张量从空间维度x-、空间维度y-和元素维度这三个维度提取其信息。然后通过三个维度的主成分信息来重新构建地球化学张量。根据重建值和原始值之间的距离来计算异常值,识别地球化学异常。
3.1. 多重主成分分析
首先,对于高维地球化学数据数据
,将其分成空间上重叠的长方体块,我们可以将这种从完整的数据中分出来的长方体块称之为“斑块”。斑块大小为
。每个斑块就是一个小型张量,表示为
,其中上标(i)表示的是第i个斑块。
对于一个三阶张量
,其n-模式展开矩阵由
的所有n-模式展开向量进行排列而成(
)。我们将地球化学数据三阶张量的1-模式、2-模式、3-模式展开矩阵称为x-模式、y-模式、元素模式展开矩阵,它们分别表示沿x、y和元素维度的展开矩阵。因此,斑块
可以被展开为三种模式的矩阵,即x-模式展开矩阵
、y-模式展开矩阵
、元素模式展开矩阵
。
对于每个模式下的展开矩阵,我们分别利用主成分分析方法从中提取信息和特征,形成主成分信息矩阵。主成分分析(principal component analysis, PCA)是将多个变量指标综合为少数几个综合性变量指标的一种统计分析方法。主成分分析可以通过正交变换,将一组可能存在相关性的变量转换为一组线性不相关的变量,转换后的这组变量叫主成分。因此,主成分能够保留原始变量的绝大部分信息,而且主成分之间是互不相关的。
对于斑块
的x-模式展开矩阵
,我们可以对其进行主成分分析,则
可以表示成下列分解形式:
(3-1)
其中,主成分数据矩阵
,
表示第i个主成分组成的向量。
是矩阵
的协方差阵的特征值,也叫做载荷(loading)矩阵。载荷矩阵是正交的,即:
。由此,空间x-维度的主成分矩阵
可以由x-模式展开矩阵与载荷矩阵的矩阵乘法得到:
(3-2)
我们取PCA结果的前p个主成分(
),也就是
的前p列,形成主成分信息矩阵
。即:
。
可以代表空间x-维度上的重要信息。
接下来,我们对斑块
的y-模式展开矩阵进行主成分分析,并且得到空间y-维度的主成分矩阵
:
(3-3)
其中,主成分数据矩阵
的每一个列是一个主成分向量,
表示载荷矩阵。与
类似,空间y维度信息矩阵
可以通过取前p个主成分(
)组成,
,
表示第i个主成分向量。由此得到代表空间y维度重要信息的矩阵
。
与空间维度上提取主成分信息矩阵过程相同,元素维度信息矩阵可以通过PCA结果中前q个主成分获得。元素维度主成分矩阵
表示为:
(3-4)
其中,主成分数据矩阵
的一列表示元素维度上的一个主成分向量,
表示载荷矩阵。元素维度信息矩阵
由PCA结果中的前q个主成分排列组成,包含着元素维度上的主要信息。
3.2. 多维信息重建
得到了空间维度上的信息矩阵
,
和元素维度信息矩阵
后,我们可以得到一个联合的信息矩阵:
(3-5)
其中符号
表示的是Kronecker乘积运算 [32]。联合的空间–元素维度信息矩阵能表示张量
的x、y和元素维度的联合特征。我们能够根据空间–元素维度联合特征,还原重建出张量的主要信息。我们将张量
的向量化形式写为
,则有
,其中
表示主成分信息矩阵
对应的重建核心系数,包含着矩阵
中各个主成分之间的关系。即,存在一组线性变换,能够将主成分信息矩阵
的各列(主成分)通过变换得到
。核心系数
可以通过最小化下列公式估计得到:
(3-6)
公式中,符号
表示的是
-norm。这里我们使用使用交替乘法(ADMM) [33] 求解出
。
主成分信息矩阵
包含着元素相关性和空间结构的信息。因此,核心系数
可以对斑块
的元素特征和空间特征进行编码,来重建地球化学背景。重建的过程中,我们需要使重建的数据
接近原始的数据
。我们先分别给定原始张量
和重建张量
的向量化形式
和
。则求上述问题可以表述为:
(3-7)
式中,
是一个由0和1组成的对角矩阵,用于从
中提取斑块
。在公式中,第一项确保关于
的信息来自于空间元素信息矩阵
和核心系数
的信息,公式第二项可以约束重建数据的向量
逼近原始数据向量
。
是加权系数,控制
与
之间的接近程度,在我们的工作中,
。
公式(3-7)中的优化是一个简单的线性最小二乘法问题,对于
的封闭式公式为:
(3-8)
是单位阵。最后,重建的地球化学数据
被转换称张量形式,即:
。
得到重建数据之后,我们为每一个样本点计算重建值与原始值之间的距离,作为其异常得分,计算公式如下: [21] [27] [34]
(3-9)
其中,
和
分别表示
和
在第k个元素位置
上的值。
4. 地球化学异常识别结果
我们以研究区39种元素地球化学数据为原始数据。地球化学元素数据是一种成分数据,我们使用等距对数比变换(isometric-logratio transformation, ilr)打开成分数据,剔除掉元素间的伪相关关系,消除成分数据结构中的“闭合效应”。然后再对ilr变换后的数据进行网格化,形成地球化学数据张量。
ilr变换后的地球化学张量被分割成16 × 16 × 38大小的97个斑块。使用PCA分析中的主成分构成空间信息矩阵和元素信息矩阵,重建地球化学张量,得到地球化学异常得分。本文用接收者操作特征(receiver operating characteristic, ROC)曲线和ROC下曲线面积(area under the receiver operating characteristic curve, AUC)来量化异常得分与已知金矿之间的空间相关性 [35] [36],评价异常提取结果。
4.1. 参数的选择
主成分分析后得到的主成分拥有者最大的信息量。在使用多重PCA分析来提取地球化学张量的空间信息矩阵和元素信息矩阵时,我们需要确定所取得主成分的数量,即确定空间信息矩阵和元素信息矩阵的列数。因此,我们用GridSearch策略寻求最佳的列数。空间特征的主成分个数搜索范围是从2到16,元素特征的主成分个数搜索范围是从2到20,使用AUC对结果对模型结果进行评价。从图2中可以看到,元素信息矩阵列数的大小为3时,AUC的值普遍比较高。而固定元素信息矩阵大小时,空间信息矩阵的列数对AUC的影响不大。
在分别学习x-与y-信息矩阵和元素信息矩阵的基础上,为了找到合适的列数,我们进一步探索空间信息矩阵和元素信息矩阵列数(即主成分个数)的影响,选取不同列数的空间信息矩阵提取元素和空间的信息,然后评估模型的结果。首先我们固定元素信息矩阵的大小,选取不同大小的空间信息矩阵进行主成分分析实验。根据Grid Search的结果(图2),我们选取列数为3的元素信息矩阵,空间信息矩阵的列数范围从2到8。如图显示,当空间信息矩阵列数为3时,基于多重PCA信息矩阵的方法能够产生最大的AUC值0.794。空间信息矩阵列数大于3时,增加主成分个数对结果并没有正向的影响。
![](//html.hanspub.org/file/11-1771421x102_hanspub.png?20220602091954271)
Figure 2. Grid-search results of the optimized atomic sizes for elemental informationmatrix and spatial informationmatrix
图2. 空间信息矩阵和元素信息矩阵的主成分个数网格搜索结果
我们再次固定空间信息矩阵的大小,选取不同大小的元素信息矩阵进行实验,使用AUC评估基多重PCA方法的结果。根据不同空间信息矩阵列数的实验结果(图3(a)),我们固定空间信息矩阵的列数为3,采用范围从2到8的主成分元素信息矩阵列数进行实验。如图3(b)显示,当使用3个主成分构成元素信息矩阵时,基于多重PCA方法的AUC达到最高值。说明在元素域3个主成分包含的信息量最适合提取地球化学异常。
下面我们采用元素信息矩阵列数(主成分个数)为3,空间信息矩阵列数(主成分个数)也为3的设置进行多重主成分分析,计算地球化学异常得分。即:p = 3,q = 3。按照这种列数大小的设置,空间–元素联合信息矩阵大小为9,728 × 27,其中9,728 = 16 × 16 × 38,对应于“斑块”的尺寸,27 = 3 × 3 × 3对应每个维度的主成分个数的乘积。即,由多重主成分分析提取的空间–元素联合信息矩阵的列数为27。
![](//html.hanspub.org/file/11-1771421x103_hanspub.png?20220602091954271)
![](//html.hanspub.org/file/11-1771421x104_hanspub.png?20220602091954271)
Figure 3. The ROC curves from spatial information matrix (a) and the elemental information matrix (b) with different numbers of principal components
图3. 不同主成分个数的空间信息矩阵(a)和元素信息矩阵(b)产生的ROC曲线
PCA方法能够变换出含有最大信息的变量。为了提升模型的性能,使主成分能充分代表性每个斑块的特征,我们用Kmeans算法,提前将97个斑块进行聚类分析,然后再对每一类的斑块进行信息矩阵学习。那么,同一个类的斑块将共享空间信息矩阵。Kmeans算法中,聚类的类别个数也是至关重要的,这里我们画出了聚类的类别个数K从1到10范围下模型结果的AUC值。从图4可以看到,当聚类数K = 5的时候,AUC值最大,ROC曲线更偏向于左上角。并且聚类的操作对AUC的提升比较大,对结果有明显优化作用。
![](//html.hanspub.org/file/11-1771421x105_hanspub.png?20220602091954271)
Figure 4. ROC curves of multiple PCA methods with different number of clusters
图4. 多重PCA方法在不同聚类数下的ROC曲线图
4.2. 地球化学异常得分及评价
为了得到最优效果,最终我们使用p = 3,q = 3的列数设置进行了地球化学异常识别,并且提前将97个“斑块”按照类别数K = 5进行聚类,分别对每一类进行信息矩阵提取。
我们根据公式(3-9)估计了异常得分。基于多重主成分分析的方法得出的异常得分显示在图5中,大部分的金矿点都处在地球化学异常得分高的区域。为了评估已确定的地球化学异常和矿化之间的空间联系,我们绘制ROC曲线和一个显示累积异常区域与累积金矿的关系的图(如图6)。
从异常分数图(图5)可以看到,识别出的异常得分与金矿点大部分是吻合的。图6(a)中的ROC曲线靠近左上角,并且ROC曲线下面积AUC的值为0.856。从图6(b)中看到,近80%的已知金矿床位于异常得分较高的前35%的区域。因此我们可以认为,基于多重主成分分析计算出的异常得分,与研究区的已知金矿点密切相关。
![](//html.hanspub.org/file/11-1771421x106_hanspub.png?20220602091954271)
Figure 5. Geochemical anomaly scores based on multiple principal component analysis
图5. 基于多重主成分分析提取的地球化学异常分数
![](//html.hanspub.org/file/11-1771421x107_hanspub.png?20220602091954271)
Figure 6. ROC curve (a), and predictive probability curve (b) based on multiple principal component analysis
图6. 基于多重主成分分析的ROC曲线(a),以及预测度曲线(b)
5. 结论
1) 本研究以中国山东省胶东半岛西北部地区为案例,对该区域的地球化学元素数据进行分析,基于多重PCA方法提取地球化学数据的空间–元素结构联合信息,由此重建地球化学高维数据。根据地球化学原始数据与高维数据之间的欧氏距离计算异常得分,实现研究区地球化学异常提取。
2) 基于多重PCA的方法从张量的角度出发,兼顾地球化学元素的空间维度和元素维度。将39种元素的地球化学数据视为张量,利用PCA提取空间维度和元素维度的主要信息,同时考虑地球化学数据的空间结构以及39种元素关联性。
3) 本研究使用AUC指标量化地球化学异常得分与已知金矿点之间的空间关联度,实现对多重PCA方法提取结果的评价。模型提取的地球化学异常得分分数高的区域,与金矿点分布相吻合;地球化学异常得分的AUC值为0.856。评价结果表明:基于多重主成分分析计算的异常点得分与已知金矿床之间存在密切的空间相关性。因此我们可以认为,基于多重主成分分析的方法能够有效提取出地球化学异常。