基于水化學(xué)數(shù)據(jù)的礦山涌水水源識(shí)別:主成分分析與殘差分析

顧鴻宇，馬鳳山，王東輝，李勝偉*，劉 港，向元英，郭子奇

(1. 中國(guó)地質(zhì)調(diào)查局成都地質(zhì)調(diào)查中心，四川 成都 610081； 2. 中國(guó)科學(xué)院地質(zhì)與地球物理研究所頁(yè)巖氣與地質(zhì)工程重點(diǎn)實(shí)驗(yàn)室，北京 100029； 3. 中國(guó)地質(zhì)調(diào)查局西安地質(zhì)調(diào)查中心，陜西 西安 710054； 4. 德陽(yáng)市環(huán)境監(jiān)測(cè)中心站，四川 德陽(yáng) 618000)

0 引 言

礦井涌水不僅給井建、巷道掘進(jìn)以及機(jī)電設(shè)備的管理帶來(lái)困難，嚴(yán)重時(shí)還可能造成淹井事故和人員傷亡。因此，查明礦山涌水的水源(端元)是一個(gè)重要課題[1-3]。運(yùn)用巖體力學(xué)和水力學(xué)來(lái)研究礦山開(kāi)采過(guò)程中地下水動(dòng)態(tài)是一種間接方法，這些方法都必須給出一些假設(shè)和簡(jiǎn)化，研究尺度往往也局限于有限尺寸，計(jì)算結(jié)果只能給出一些規(guī)律性的結(jié)論，對(duì)礦山開(kāi)采的指導(dǎo)意義不是很大。利用水化學(xué)來(lái)研究地下水動(dòng)態(tài)可覆蓋整個(gè)礦山范圍，是一種直觀且經(jīng)濟(jì)可行的方法[4-6]。地下水水化學(xué)特征往往是地下水混合、水巖反應(yīng)、蒸發(fā)作用及離子交換綜合作用的結(jié)果，利用地下水離子特征可以對(duì)其物理化學(xué)過(guò)程進(jìn)行直觀判斷并得出可靠的結(jié)論[7-9]。

水源識(shí)別的第一步需要確定水源數(shù)量。目前，少有學(xué)者研究如何利用數(shù)據(jù)自約束來(lái)求取水源數(shù)量，更多的研究都是在已知水源數(shù)量和類(lèi)型后，利用訓(xùn)練樣本對(duì)未知樣本進(jìn)行分類(lèi)，在分類(lèi)過(guò)程中人為地引入判別標(biāo)志和閾值。這個(gè)過(guò)程稱(chēng)為監(jiān)督學(xué)習(xí)，如極限學(xué)習(xí)機(jī)[10]、熵權(quán)-模糊可變集理論[11]、支持向量機(jī)[12-13]、聚類(lèi)算法[14]等。實(shí)際情況中，水源數(shù)量并不等于研究區(qū)所有可能水源數(shù)量，基于概念模型來(lái)預(yù)先推測(cè)的水源數(shù)量仍然需要水化學(xué)數(shù)據(jù)來(lái)證明，因此，上述方法在確定水源數(shù)量的問(wèn)題上依然有所不足。采用多元統(tǒng)計(jì)的特征值分析方法來(lái)確定水源數(shù)量是一個(gè)常用方法。Vazquez等提出水源數(shù)量應(yīng)等于特征值大于1的特征值個(gè)數(shù)加1[15]。Long等采用基于主成分分析和聚類(lèi)的方法來(lái)確定水源數(shù)量，將主成分分析結(jié)果投影到坐標(biāo)系中，再利用聚類(lèi)來(lái)判斷聚類(lèi)數(shù)，最后水源數(shù)量就是聚類(lèi)數(shù)[16]，但是這種方法中的聚類(lèi)數(shù)是根據(jù)歐式距離來(lái)確定的，因此，聚類(lèi)數(shù)的多少有很多的人為不確定性。Vazquez等采用反演模型來(lái)計(jì)算混合比[15]，水源數(shù)量的確定是通過(guò)概念模型來(lái)確定，在反演過(guò)程中，模型收斂性問(wèn)題要求水源數(shù)量盡量小，不能出現(xiàn)共線性水源[17]，否則模型效果較差，但是共線性水源在現(xiàn)實(shí)中很可能存在且對(duì)涌水的貢獻(xiàn)不可忽略，這是該反演模型的一個(gè)缺點(diǎn)。Laaksoharju等提出的M3法中關(guān)于水源的選擇是基于各個(gè)水源樣品在主成分上的得分，選擇依據(jù)是結(jié)合實(shí)地調(diào)查選擇盡可能少的水源來(lái)描述所有樣品，即水源組成的多邊形能夠包含所有水源樣品[18]。水源識(shí)別的第二步是要確定水源類(lèi)型，也就是確定水源的水化學(xué)特征。傳統(tǒng)水源類(lèi)型的確定是對(duì)研究區(qū)地質(zhì)、構(gòu)造、水文和水文地球化學(xué)等情況進(jìn)行綜合了解，但地下含水層的各向異性導(dǎo)致這種定性描述并不能讓人信服。

綜上所述，實(shí)際研究中混合樣品非常容易取得，問(wèn)題的復(fù)雜性和難點(diǎn)在于：①各種水源樣品很難取得，因此，水源特征只能根據(jù)混合樣品來(lái)反推；②水源數(shù)量往往不清楚；③水源樣品的參數(shù)指標(biāo)很可能隨著時(shí)間和空間發(fā)生隨機(jī)變化；④水化學(xué)數(shù)據(jù)的高維度和復(fù)雜性。為此，本文提出了結(jié)合主成分分析和殘差分析的方法來(lái)確定水源數(shù)量和水源類(lèi)型。

1 分析方法

1.1 主成分分析

主成分分析(PCA)是一種標(biāo)準(zhǔn)的數(shù)據(jù)分析方法，被廣泛應(yīng)用于很多科學(xué)領(lǐng)域，如神經(jīng)系統(tǒng)科學(xué)、計(jì)算機(jī)制圖、人口統(tǒng)計(jì)學(xué)和水化學(xué)數(shù)據(jù)分析等。數(shù)據(jù)的冗余是阻礙水化學(xué)數(shù)據(jù)解譯的重要因素。數(shù)據(jù)冗余度表現(xiàn)在兩個(gè)變量之間就是相關(guān)關(guān)系：對(duì)于相關(guān)性較好的兩個(gè)變量，利用一個(gè)變量就可以預(yù)測(cè)另一個(gè)變量，數(shù)據(jù)冗余度較高；對(duì)于無(wú)相關(guān)性的兩個(gè)變量，數(shù)據(jù)冗余度較低。對(duì)于二維數(shù)據(jù)，冗余度可以通過(guò)擬合數(shù)據(jù)來(lái)判斷；對(duì)于高維數(shù)據(jù)集，則需要通過(guò)協(xié)方差矩陣來(lái)判斷。協(xié)方差矩陣中，正值代表兩個(gè)變量之間成正相關(guān)，負(fù)值代表兩個(gè)變量之間成負(fù)相關(guān)。矩陣中元素的絕對(duì)值代表了數(shù)據(jù)冗余度。為了最小化數(shù)據(jù)冗余度和最大化信噪比，有必要對(duì)協(xié)方差矩陣進(jìn)行優(yōu)化。理想的最優(yōu)協(xié)方差矩陣應(yīng)該滿足以下兩個(gè)條件：①最小化數(shù)據(jù)冗余度，要求協(xié)方差矩陣中非對(duì)角元素的值為0，即協(xié)方差矩陣為對(duì)角陣；②協(xié)方差矩陣中的每一連續(xù)維度應(yīng)按照方差的大小排列。主成分分析提供了一條最簡(jiǎn)便的思路，其假設(shè)所有的基矢量都是標(biāo)準(zhǔn)正交的。在多維(m維)數(shù)據(jù)集中，可以表述為以下步驟：①首先在m維空間中選擇一個(gè)方差最大的方向，將其記為P1；②選擇另一個(gè)和P1正交且方差最大的方向，記為P2；③重復(fù)第②步直到m個(gè)向量選擇完成。P1，P2，…，Pm依次為方差遞減的主成分，其解釋元數(shù)據(jù)集的信息比例也是遞減的。

在進(jìn)行主成分分析的過(guò)程中，應(yīng)當(dāng)注意主成分分析的第一步是要選擇好所用的變量組合，也就是選擇不同的離子組合，當(dāng)然也包括一些物理參數(shù)(電導(dǎo)率及pH值等)，不同的組合對(duì)主成分分析的分析結(jié)果有一定影響?？梢圆捎脙蓚€(gè)指標(biāo)來(lái)評(píng)價(jià)組合的適應(yīng)性：取樣適切性量數(shù)(Kaiser-Meyer-Olkin統(tǒng)計(jì)量)和巴特利(Bartlett)球形度檢驗(yàn)。

1.2 殘差分析

主成分分析的分析結(jié)果往往是由幾個(gè)主成分構(gòu)成。那么，在實(shí)際應(yīng)用過(guò)程中就會(huì)面臨著保留幾個(gè)主成分就足以代表原始數(shù)據(jù)的問(wèn)題。大多數(shù)準(zhǔn)則是保留特征值大于1的主成分，而舍去特征值小于1的主成分。這種準(zhǔn)則可能遺漏一些有價(jià)值的符合實(shí)際的水源。因此，本文用殘差分析的方法來(lái)檢驗(yàn)主成分保留個(gè)數(shù)的合理性[19-20]。通常有用的信息往往存在著內(nèi)在關(guān)系，主成分分析將這些具有相關(guān)性的信息按線性組合的方式進(jìn)行壓縮。因此，先提取任意多個(gè)主成分對(duì)離子濃度進(jìn)行重構(gòu)，再得到原始數(shù)據(jù)和重構(gòu)數(shù)據(jù)之間的差值。當(dāng)殘差表現(xiàn)出結(jié)構(gòu)性信息時(shí)，說(shuō)明提取的主成分個(gè)數(shù)不足以代表原始數(shù)據(jù)信息；而當(dāng)殘差表現(xiàn)出隨機(jī)性時(shí)，這時(shí)候殘差就代表了數(shù)據(jù)中的隨機(jī)噪聲，包括離子實(shí)驗(yàn)誤差、采樣誤差等。殘差分析步驟如下：首先對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，以消除地下水離子組分因濃度數(shù)量級(jí)差異造成的不同離子權(quán)重不相同的現(xiàn)象。其表達(dá)式為

x＇ij=(xij-j)/Sj

(1)

式中:x＇ij為標(biāo)準(zhǔn)化后的離子濃度;xij為原始離子濃度;j為第j種離子濃度均值;Sj為第j種離子濃度的標(biāo)準(zhǔn)差。

將樣品數(shù)據(jù)以濃度的形式表示在主成分分析的前m維子空間上，其表達(dá)式為

(2)

式中:v表示前m個(gè)特征向量;ij表示重構(gòu)的離子濃度。

那么殘差(R)就能表示為

R=ij-xij

(3)

同時(shí)，由于各種離子濃度的絕對(duì)值相差很大，利用偏差(RB)和相對(duì)平方根偏差(RRMSE)來(lái)度量重構(gòu)的精確性。其表達(dá)式為

(4)

(5)

式中：bj為第j種離子濃度的偏差；rj為第j種離子濃度的相對(duì)平方根偏差。

2 實(shí)例分析

2.1 礦山地質(zhì)概況

本實(shí)例中礦山為濱海金礦，西北及北方向?yàn)l臨渤海，東南及南方向?yàn)殛懙?圖1)。該區(qū)域地勢(shì)低洼，大部分地方地形平坦，平均海拔為1.2～4.5 m，最高的3座山峰位于斷層F1下盤(pán)，海拔約67 m。礦區(qū)部分位于海底。

第一含水層主要由中砂、粗砂組成，局部地段有細(xì)砂和礫石出現(xiàn)，地表0～2 m范圍內(nèi)，粒度均勻，向下出現(xiàn)貝殼等有機(jī)物；該含水層厚度為3.50～7.29 m，其富水性因含泥質(zhì)的多少而有較大差異，富水性較強(qiáng)區(qū)域的單位涌水量達(dá)15.27 L·s-1，滲透系數(shù)為117.46 m·d-1；該含水層主要接受大氣降水和海水潮汐補(bǔ)給；水位埋深為0.5～6.0 m，礦區(qū)南部埋深淺；礦化度從東至西、從北至南逐漸增高。第一隔水層位于第一含水層下，埋深為5.5～9.0 m，主要為砂質(zhì)黏土、黏質(zhì)砂土和含鈣質(zhì)結(jié)核砂質(zhì)黏土，部分地方夾砂及礫石的含水透鏡體；該層厚度變化不大，一般為7～8 m。第二含水層位于第一隔水層之下，該層不連續(xù)，主要分布在斷層F3上部以及東北方向；該層巖性主要為中砂、粗砂和礫石；厚度由北往南逐漸變厚，一般為3～4 m，最大厚度達(dá)11.9 m；該層地下水具有一定承壓性，可接受海水和上部第一含水層的補(bǔ)給。第二隔水層位于風(fēng)化殼上，埋深為7.8～25.0 m，厚度穩(wěn)定，一般為3～5 m，最厚達(dá)19.6 m，主要為黃棕色含礫砂質(zhì)黏土以及紅棕色黏土，比第一隔水層的隔水性能更好[21]。

圖1 礦山地質(zhì)概況及采樣點(diǎn)分布Fig.1 Mine Geological Background and Distribution of Sampling Sites

礦區(qū)主要出露的巖漿巖為中生代郭家?guī)X花崗閃長(zhǎng)巖和玲瓏二長(zhǎng)花崗巖，局部有以煌斑巖和輝綠巖為主的基性巖脈發(fā)育。礦區(qū)圍巖有強(qiáng)烈蝕變和分帶特點(diǎn)，總體上從礦體向圍巖依次為黃鐵絹英質(zhì)碎斑巖、絹英巖、鉀硅化蝕變花崗巖、新鮮花崗巖等。

礦區(qū)同時(shí)存在著3條斷層,分別為F1、F2和F3。斷層F1為礦區(qū)的控礦斷裂，礦體分布在斷層下盤(pán)，寬度為50～200 m,平面上呈S型，總體走向?yàn)镹E(40°)，傾向SE，傾角約45°，主裂面呈波狀舒緩，具扭壓性，斷層面發(fā)育5～10 cm厚灰黑色斷層泥。

在斷層F1西側(cè)有一條規(guī)模較小的伴生斷層F2。該斷層走向?yàn)?80°，傾角為85°，上盤(pán)北移，下盤(pán)南移，具有扭性斷層特點(diǎn)。斷層兩側(cè)裂隙發(fā)育特征差異明顯，西側(cè)NE向裂隙較發(fā)育，NW向裂隙少見(jiàn)，東側(cè)主要發(fā)育NW向裂隙，靠近斷層F2裂隙發(fā)育密集。根據(jù)物探結(jié)果分析，斷層F2表現(xiàn)為明顯的低阻特性，表明其具有良好的導(dǎo)水性[22]。

斷層F3與F1近直交并切斷斷層F1(圖1)，將斷層F1錯(cuò)斷10～20 m。斷層F3為一條橫穿整個(gè)礦區(qū)的區(qū)域性斷層，即三元—陳家大斷裂，總體走向?yàn)?00°～310°，傾向NE，傾角近90°。斷層F3主裂面位于南側(cè)，有泥質(zhì)及斷層角礫填充，南側(cè)巖體完整性較好，主要發(fā)育NW—SE向節(jié)理，節(jié)理面平直。斷層北側(cè)發(fā)育有NE—SW向、SE—NW向及水平節(jié)理，巖體破碎，節(jié)理開(kāi)度大。根據(jù)目前采礦勘探，其深度應(yīng)大于850 m，破碎帶寬度為15～35 m，屬于張性斷裂。斷層面未發(fā)現(xiàn)充填物，具有良好的導(dǎo)水性[23]。

2.2 樣品采集與結(jié)果分析

2.3 水源數(shù)量和水源類(lèi)型

表1 主成分分析結(jié)果

圖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離子相關(guān)性Fig.3 Relationships of Ca-Cl and Mg-Cl

表2 不同主成分個(gè)數(shù)重構(gòu)離子濃度和氫氧同位素的偏差和相對(duì)平方根偏差

圖4 保留第一主成分后殘差與原始離子濃度、氫氧同位素的相關(guān)性Fig.4 Relationships Between Residuals and Measured Ion Concentrations, Hydrogen-oxygen Isotopes After Retaining the First Principal Component

圖5 保留前2個(gè)主成分后殘差與原始離子濃度、氫氧同位素的相關(guān)性Fig.5 Relationships Between Residuals and Measured Ion Concentrations, Hydrogen-oxygen Isotopes After Retaining the First Two Principal Components

圖6 保留前3個(gè)主成分后殘差與原始離子濃度、氫氧同位素的相關(guān)性Fig.6 Relationships Between Residuals and Measured Ion Concentrations, Hydrogen-oxygen Isotopes After Retaining the First Three Principal Components

圖8 保留前5個(gè)主成分后殘差與原始離子濃度、氫氧同位素的相關(guān)性Fig.8 Relationships Between Residuals and Measured Ion Concentrations, Hydrogen-oxygen Isotopes After Retaining the First Five Principal Components

在主成分分析結(jié)果中，各個(gè)因子的荷載有正有負(fù)，因子對(duì)某一主成分的物理意義就顯得模糊不清。為了解釋每個(gè)主成分在地下水水化學(xué)信息演化過(guò)程中的作用，利用逐漸增加主成分個(gè)數(shù)來(lái)對(duì)比重構(gòu)離子濃度和原始離子濃度之間的差別是一種有效的方法。然后再結(jié)合離子相關(guān)性及離子比率分析，可以確定每一主成分代表的地下水演化過(guò)程(圖9)。

圖(b)中，a～h分別對(duì)應(yīng)采樣周期2009年、2011～2017年；kf代表2015年冬季采樣圖9 2009～2017年樣品地下水演化過(guò)程Fig.9 Evolution Processes of Underground Water for Samples in 2009-2017

第一主成分解釋了54.836%的原始變量信息，表現(xiàn)出對(duì)Na+和Cl-的有效解釋(圖4)。這兩種離子在地下水中的離子濃度高度相關(guān)，且隨著深度的增加而增加。在采礦過(guò)程中低Na和低Cl的淺部地下水向深部高濃度區(qū)補(bǔ)給，地下水發(fā)生Na+與Ca2+、Mg2+的離子交換可以忽略，Na+和Cl-可視為地下水中的穩(wěn)定離子。因此，第一主成分的物理意義可以理解為地下水的大部分原始來(lái)源是海水，地下水形成過(guò)程中主要發(fā)生蒸發(fā)作用。

第二主成分解釋了同位素和K+的大部分信息，但仍有部分結(jié)構(gòu)化信息(圖5)。注意到這3種元素在海水中的濃度很高，且高于地下水中的濃度。因此，第二主成分很好地解釋了海水的特征。同時(shí)，第二主成分的方差貢獻(xiàn)率較高(22.845%)，說(shuō)明海水對(duì)地下水的補(bǔ)給貢獻(xiàn)較大。

第三主成分解釋了10.312%的原始變量信息，包含了Ca2+和Mg2+的演化信息(圖6)。從圖3可以看出，Ca2+、Mg2+在Cl-濃度為20 000 mg·L-1時(shí)出現(xiàn)了明顯的Ca富集和Mg貧化。從多期離子濃度演化過(guò)程(圖9)可以得出,該主成分代表的化學(xué)演化為Ca2+和Mg2+的離子交換作用，富Mg的基巖裂隙水逐漸演化為富Ca的裂隙水[24-25]。

第四主成分主要消除了δD值和K+濃度殘差的結(jié)構(gòu)性，而對(duì)其他離子濃度殘差影響不大(圖7)。這個(gè)主成分可以認(rèn)為是第二主成分的補(bǔ)充，其對(duì)地下水演化過(guò)程具有一些相同的作用。海水和第四系孔隙水的化學(xué)特征有著一些相似之處，如K+濃度較高，同時(shí)第四系孔隙水和海水具有成因聯(lián)系，第四系孔隙水在形成過(guò)程中海水的補(bǔ)給比例較大。因此，第四主成分應(yīng)代表地下水向第四系孔隙水方向演化的作用。

3 結(jié) 語(yǔ)

(1)主成分分析在水化學(xué)數(shù)據(jù)的處理上具有良好的效率和正確性，其提取出的結(jié)構(gòu)化信息(各個(gè)主成分)能夠很好地解釋地下水演化中的化學(xué)過(guò)程和物理過(guò)程。值得注意的是，這并不代表在分析水化學(xué)數(shù)據(jù)的過(guò)程中可以直接忽略傳統(tǒng)的水化學(xué)分析而直接進(jìn)行主成分分析。傳統(tǒng)的水化學(xué)分析可以提供水化學(xué)演化過(guò)程的類(lèi)型和程度，可以把其分析結(jié)果看作是主成分分析的先驗(yàn)知識(shí)，利用這些先驗(yàn)知識(shí)，確定主成分分析的主成分個(gè)數(shù)，驗(yàn)證主成分分析結(jié)果的正確性，將數(shù)據(jù)不失真地進(jìn)行壓縮，從而為在較小維度上的模型搭建提供方便。

(2)殘差分析巧妙地將主成分分析結(jié)果以重構(gòu)的離子濃度展現(xiàn)出來(lái)，并將重構(gòu)離子濃度和原始離子濃度的差值與原始離子濃度進(jìn)行相關(guān)分析。當(dāng)所有離子濃度殘差表現(xiàn)出足夠的隨機(jī)分布時(shí)，表明有效信息的完全提取，這時(shí)保留的主成分個(gè)數(shù)即為所有的涌水水源數(shù)量。同時(shí)，通過(guò)逐漸增加主成分個(gè)數(shù)，有利于對(duì)每一主成分進(jìn)行有效和合理的解釋?zhuān)创_定地下水的演化過(guò)程及水源類(lèi)型。

(3)將基于主成分分析和殘差分析的方法應(yīng)用到實(shí)際礦山確定涌水水源數(shù)量和水源類(lèi)型問(wèn)題上，得到了比傳統(tǒng)方法更加合理的結(jié)果。本實(shí)例確定了研究礦山的涌水水源數(shù)量為5個(gè)，水源類(lèi)型為海水、第四系孔隙水、富鈣基巖水、富鎂基巖水和淡水。