基于水化学数据的矿山涌水水源识别:主成分分析与残差分析

顾鸿宇，马凤山，王东辉，李胜伟*，刘 港，向元英，郭子奇

(1. 中国地质调查局成都地质调查中心，四川 成都 610081； 2. 中国科学院地质与地球物理研究所页岩气与地质工程重点实验室，北京 100029； 3. 中国地质调查局西安地质调查中心，陕西 西安 710054； 4. 德阳市环境监测中心站，四川 德阳 618000)

0 引 言

矿井涌水不仅给井建、巷道掘进以及机电设备的管理带来困难，严重时还可能造成淹井事故和人员伤亡。因此，查明矿山涌水的水源(端元)是一个重要课题[1-3]。运用岩体力学和水力学来研究矿山开采过程中地下水动态是一种间接方法，这些方法都必须给出一些假设和简化，研究尺度往往也局限于有限尺寸，计算结果只能给出一些规律性的结论，对矿山开采的指导意义不是很大。利用水化学来研究地下水动态可覆盖整个矿山范围，是一种直观且经济可行的方法[4-6]。地下水水化学特征往往是地下水混合、水岩反应、蒸发作用及离子交换综合作用的结果，利用地下水离子特征可以对其物理化学过程进行直观判断并得出可靠的结论[7-9]。

水源识别的第一步需要确定水源数量。目前，少有学者研究如何利用数据自约束来求取水源数量，更多的研究都是在已知水源数量和类型后，利用训练样本对未知样本进行分类，在分类过程中人为地引入判别标志和阈值。这个过程称为监督学习，如极限学习机[10]、熵权-模糊可变集理论[11]、支持向量机[12-13]、聚类算法[14]等。实际情况中，水源数量并不等于研究区所有可能水源数量，基于概念模型来预先推测的水源数量仍然需要水化学数据来证明，因此，上述方法在确定水源数量的问题上依然有所不足。采用多元统计的特征值分析方法来确定水源数量是一个常用方法。Vazquez等提出水源数量应等于特征值大于1的特征值个数加1[15]。Long等采用基于主成分分析和聚类的方法来确定水源数量，将主成分分析结果投影到坐标系中，再利用聚类来判断聚类数，最后水源数量就是聚类数[16]，但是这种方法中的聚类数是根据欧式距离来确定的，因此，聚类数的多少有很多的人为不确定性。Vazquez等采用反演模型来计算混合比[15]，水源数量的确定是通过概念模型来确定，在反演过程中，模型收敛性问题要求水源数量尽量小，不能出现共线性水源[17]，否则模型效果较差，但是共线性水源在现实中很可能存在且对涌水的贡献不可忽略，这是该反演模型的一个缺点。Laaksoharju等提出的M3法中关于水源的选择是基于各个水源样品在主成分上的得分，选择依据是结合实地调查选择尽可能少的水源来描述所有样品，即水源组成的多边形能够包含所有水源样品[18]。水源识别的第二步是要确定水源类型，也就是确定水源的水化学特征。传统水源类型的确定是对研究区地质、构造、水文和水文地球化学等情况进行综合了解，但地下含水层的各向异性导致这种定性描述并不能让人信服。

综上所述，实际研究中混合样品非常容易取得，问题的复杂性和难点在于：①各种水源样品很难取得，因此，水源特征只能根据混合样品来反推；②水源数量往往不清楚；③水源样品的参数指标很可能随着时间和空间发生随机变化；④水化学数据的高维度和复杂性。为此，本文提出了结合主成分分析和残差分析的方法来确定水源数量和水源类型。

1 分析方法

1.1 主成分分析

主成分分析(PCA)是一种标准的数据分析方法，被广泛应用于很多科学领域，如神经系统科学、计算机制图、人口统计学和水化学数据分析等。数据的冗余是阻碍水化学数据解译的重要因素。数据冗余度表现在两个变量之间就是相关关系：对于相关性较好的两个变量，利用一个变量就可以预测另一个变量，数据冗余度较高；对于无相关性的两个变量，数据冗余度较低。对于二维数据，冗余度可以通过拟合数据来判断；对于高维数据集，则需要通过协方差矩阵来判断。协方差矩阵中，正值代表两个变量之间成正相关，负值代表两个变量之间成负相关。矩阵中元素的绝对值代表了数据冗余度。为了最小化数据冗余度和最大化信噪比，有必要对协方差矩阵进行优化。理想的最优协方差矩阵应该满足以下两个条件：①最小化数据冗余度，要求协方差矩阵中非对角元素的值为0，即协方差矩阵为对角阵；②协方差矩阵中的每一连续维度应按照方差的大小排列。主成分分析提供了一条最简便的思路，其假设所有的基矢量都是标准正交的。在多维(m维)数据集中，可以表述为以下步骤：①首先在m维空间中选择一个方差最大的方向，将其记为P1；②选择另一个和P1正交且方差最大的方向，记为P2；③重复第②步直到m个向量选择完成。P1，P2，…，Pm依次为方差递减的主成分，其解释元数据集的信息比例也是递减的。

在进行主成分分析的过程中，应当注意主成分分析的第一步是要选择好所用的变量组合，也就是选择不同的离子组合，当然也包括一些物理参数(电导率及pH值等)，不同的组合对主成分分析的分析结果有一定影响。可以采用两个指标来评价组合的适应性：取样适切性量数(Kaiser-Meyer-Olkin统计量)和巴特利(Bartlett)球形度检验。

1.2 残差分析

主成分分析的分析结果往往是由几个主成分构成。那么，在实际应用过程中就会面临着保留几个主成分就足以代表原始数据的问题。大多数准则是保留特征值大于1的主成分，而舍去特征值小于1的主成分。这种准则可能遗漏一些有价值的符合实际的水源。因此，本文用残差分析的方法来检验主成分保留个数的合理性[19-20]。通常有用的信息往往存在着内在关系，主成分分析将这些具有相关性的信息按线性组合的方式进行压缩。因此，先提取任意多个主成分对离子浓度进行重构，再得到原始数据和重构数据之间的差值。当残差表现出结构性信息时，说明提取的主成分个数不足以代表原始数据信息；而当残差表现出随机性时，这时候残差就代表了数据中的随机噪声，包括离子实验误差、采样误差等。残差分析步骤如下：首先对数据进行标准化，以消除地下水离子组分因浓度数量级差异造成的不同离子权重不相同的现象。其表达式为

x＇ij=(xij-j)/Sj

(1)

式中:x＇ij为标准化后的离子浓度;xij为原始离子浓度;j为第j种离子浓度均值;Sj为第j种离子浓度的标准差。

将样品数据以浓度的形式表示在主成分分析的前m维子空间上，其表达式为

(2)

式中:v表示前m个特征向量;ij表示重构的离子浓度。

那么残差(R)就能表示为

R=ij-xij

(3)

同时，由于各种离子浓度的绝对值相差很大，利用偏差(RB)和相对平方根偏差(RRMSE)来度量重构的精确性。其表达式为

(4)

(5)

式中：bj为第j种离子浓度的偏差；rj为第j种离子浓度的相对平方根偏差。

2 实例分析

2.1 矿山地质概况

本实例中矿山为滨海金矿，西北及北方向濒临渤海，东南及南方向为陆地(图1)。该区域地势低洼，大部分地方地形平坦，平均海拔为1.2～4.5 m，最高的3座山峰位于断层F1下盘，海拔约67 m。矿区部分位于海底。

第一含水层主要由中砂、粗砂组成，局部地段有细砂和砾石出现，地表0～2 m范围内，粒度均匀，向下出现贝壳等有机物；该含水层厚度为3.50～7.29 m，其富水性因含泥质的多少而有较大差异，富水性较强区域的单位涌水量达15.27 L·s-1，渗透系数为117.46 m·d-1；该含水层主要接受大气降水和海水潮汐补给；水位埋深为0.5～6.0 m，矿区南部埋深浅；矿化度从东至西、从北至南逐渐增高。第一隔水层位于第一含水层下，埋深为5.5～9.0 m，主要为砂质黏土、黏质砂土和含钙质结核砂质黏土，部分地方夹砂及砾石的含水透镜体；该层厚度变化不大，一般为7～8 m。第二含水层位于第一隔水层之下，该层不连续，主要分布在断层F3上部以及东北方向；该层岩性主要为中砂、粗砂和砾石；厚度由北往南逐渐变厚，一般为3～4 m，最大厚度达11.9 m；该层地下水具有一定承压性，可接受海水和上部第一含水层的补给。第二隔水层位于风化壳上，埋深为7.8～25.0 m，厚度稳定，一般为3～5 m，最厚达19.6 m，主要为黄棕色含砾砂质黏土以及红棕色黏土，比第一隔水层的隔水性能更好[21]。

图1 矿山地质概况及采样点分布Fig.1 Mine Geological Background and Distribution of Sampling Sites

矿区主要出露的岩浆岩为中生代郭家岭花岗闪长岩和玲珑二长花岗岩，局部有以煌斑岩和辉绿岩为主的基性岩脉发育。矿区围岩有强烈蚀变和分带特点，总体上从矿体向围岩依次为黄铁绢英质碎斑岩、绢英岩、钾硅化蚀变花岗岩、新鲜花岗岩等。

矿区同时存在着3条断层,分别为F1、F2和F3。断层F1为矿区的控矿断裂，矿体分布在断层下盘，宽度为50～200 m,平面上呈S型，总体走向为NE(40°)，倾向SE，倾角约45°，主裂面呈波状舒缓，具扭压性，断层面发育5～10 cm厚灰黑色断层泥。

在断层F1西侧有一条规模较小的伴生断层F2。该断层走向为280°，倾角为85°，上盘北移，下盘南移，具有扭性断层特点。断层两侧裂隙发育特征差异明显，西侧NE向裂隙较发育，NW向裂隙少见，东侧主要发育NW向裂隙，靠近断层F2裂隙发育密集。根据物探结果分析，断层F2表现为明显的低阻特性，表明其具有良好的导水性[22]。

断层F3与F1近直交并切断断层F1(图1)，将断层F1错断10～20 m。断层F3为一条横穿整个矿区的区域性断层，即三元—陈家大断裂，总体走向为300°～310°，倾向NE，倾角近90°。断层F3主裂面位于南侧，有泥质及断层角砾填充，南侧岩体完整性较好，主要发育NW—SE向节理，节理面平直。断层北侧发育有NE—SW向、SE—NW向及水平节理，岩体破碎，节理开度大。根据目前采矿勘探，其深度应大于850 m，破碎带宽度为15～35 m，属于张性断裂。断层面未发现充填物，具有良好的导水性[23]。

2.2 样品采集与结果分析

2.3 水源数量和水源类型

表1 主成分分析结果

图2 所有样品在PC1-PC3平面和PC2-PC3平面上的投影Fig.2 Projections of All Samples on the PC1-PC3 and PC2-PC3 Planes

c(·)为某离子浓度图3 Ca-Cl和Mg-Cl离子相关性Fig.3 Relationships of Ca-Cl and Mg-Cl

表2 不同主成分个数重构离子浓度和氢氧同位素的偏差和相对平方根偏差

图4 保留第一主成分后残差与原始离子浓度、氢氧同位素的相关性Fig.4 Relationships Between Residuals and Measured Ion Concentrations, Hydrogen-oxygen Isotopes After Retaining the First Principal Component

图5 保留前2个主成分后残差与原始离子浓度、氢氧同位素的相关性Fig.5 Relationships Between Residuals and Measured Ion Concentrations, Hydrogen-oxygen Isotopes After Retaining the First Two Principal Components

图6 保留前3个主成分后残差与原始离子浓度、氢氧同位素的相关性Fig.6 Relationships Between Residuals and Measured Ion Concentrations, Hydrogen-oxygen Isotopes After Retaining the First Three Principal Components

图8 保留前5个主成分后残差与原始离子浓度、氢氧同位素的相关性Fig.8 Relationships Between Residuals and Measured Ion Concentrations, Hydrogen-oxygen Isotopes After Retaining the First Five Principal Components

在主成分分析结果中，各个因子的荷载有正有负，因子对某一主成分的物理意义就显得模糊不清。为了解释每个主成分在地下水水化学信息演化过程中的作用，利用逐渐增加主成分个数来对比重构离子浓度和原始离子浓度之间的差别是一种有效的方法。然后再结合离子相关性及离子比率分析，可以确定每一主成分代表的地下水演化过程(图9)。

图(b)中，a～h分别对应采样周期2009年、2011～2017年；kf代表2015年冬季采样图9 2009～2017年样品地下水演化过程Fig.9 Evolution Processes of Underground Water for Samples in 2009-2017

第一主成分解释了54.836%的原始变量信息，表现出对Na+和Cl-的有效解释(图4)。这两种离子在地下水中的离子浓度高度相关，且随着深度的增加而增加。在采矿过程中低Na和低Cl的浅部地下水向深部高浓度区补给，地下水发生Na+与Ca2+、Mg2+的离子交换可以忽略，Na+和Cl-可视为地下水中的稳定离子。因此，第一主成分的物理意义可以理解为地下水的大部分原始来源是海水，地下水形成过程中主要发生蒸发作用。

第二主成分解释了同位素和K+的大部分信息，但仍有部分结构化信息(图5)。注意到这3种元素在海水中的浓度很高，且高于地下水中的浓度。因此，第二主成分很好地解释了海水的特征。同时，第二主成分的方差贡献率较高(22.845%)，说明海水对地下水的补给贡献较大。

第三主成分解释了10.312%的原始变量信息，包含了Ca2+和Mg2+的演化信息(图6)。从图3可以看出，Ca2+、Mg2+在Cl-浓度为20 000 mg·L-1时出现了明显的Ca富集和Mg贫化。从多期离子浓度演化过程(图9)可以得出,该主成分代表的化学演化为Ca2+和Mg2+的离子交换作用，富Mg的基岩裂隙水逐渐演化为富Ca的裂隙水[24-25]。

第四主成分主要消除了δD值和K+浓度残差的结构性，而对其他离子浓度残差影响不大(图7)。这个主成分可以认为是第二主成分的补充，其对地下水演化过程具有一些相同的作用。海水和第四系孔隙水的化学特征有着一些相似之处，如K+浓度较高，同时第四系孔隙水和海水具有成因联系，第四系孔隙水在形成过程中海水的补给比例较大。因此，第四主成分应代表地下水向第四系孔隙水方向演化的作用。

3 结 语

(1)主成分分析在水化学数据的处理上具有良好的效率和正确性，其提取出的结构化信息(各个主成分)能够很好地解释地下水演化中的化学过程和物理过程。值得注意的是，这并不代表在分析水化学数据的过程中可以直接忽略传统的水化学分析而直接进行主成分分析。传统的水化学分析可以提供水化学演化过程的类型和程度，可以把其分析结果看作是主成分分析的先验知识，利用这些先验知识，确定主成分分析的主成分个数，验证主成分分析结果的正确性，将数据不失真地进行压缩，从而为在较小维度上的模型搭建提供方便。

(2)残差分析巧妙地将主成分分析结果以重构的离子浓度展现出来，并将重构离子浓度和原始离子浓度的差值与原始离子浓度进行相关分析。当所有离子浓度残差表现出足够的随机分布时，表明有效信息的完全提取，这时保留的主成分个数即为所有的涌水水源数量。同时，通过逐渐增加主成分个数，有利于对每一主成分进行有效和合理的解释，即确定地下水的演化过程及水源类型。

(3)将基于主成分分析和残差分析的方法应用到实际矿山确定涌水水源数量和水源类型问题上，得到了比传统方法更加合理的结果。本实例确定了研究矿山的涌水水源数量为5个，水源类型为海水、第四系孔隙水、富钙基岩水、富镁基岩水和淡水。