臥龍湖礦8101面地下水水質(zhì)分類判別研究

姚輝，李朝陽，張永生

(合肥工業(yè)大學(xué)資源與環(huán)境工程學(xué)院，安徽合肥230009)

地下水系統(tǒng)的水化學(xué)分類以及判別分析是水文地質(zhì)學(xué)的重要內(nèi)容之一，為了快速判別突水水源，實(shí)現(xiàn)煤炭生產(chǎn)的高產(chǎn)高效，許多學(xué)者從不同側(cè)面研究水化學(xué)模型，取得了一定的效果［1－4］，其中最基本方法之一是多元統(tǒng)計(jì)方法［5］，被證實(shí)具有一定的穩(wěn)定性。但每種方法都有局限性，考慮到礦井地下水運(yùn)動(dòng)的空間和時(shí)間變異特性，如何提高水源識(shí)別分析精度仍然是一個(gè)需要長期解決的問題。本文以臥龍湖礦8101工作面為例，先采用多元統(tǒng)計(jì)方法建立水化學(xué)判別模型，主要包括主成分分析、聚類分析、判別分析，建立水化學(xué)判別模型，對(duì)未知水樣做出判定。同時(shí)與Maltab編程方法相結(jié)合［6－8］，對(duì)未知水樣進(jìn)行判斷，二者進(jìn)行對(duì)比。

一礦區(qū)地質(zhì)、水文地質(zhì)條件

臥龍湖煤礦位于安徽省濉溪縣鐵佛、岳集境內(nèi)，東距百善煤礦約15 km，北以省界與河南省永城市毗鄰，南北長約8～9 km，東西寬約3.5～4 km，面積28.85 km2。該煤田位于華北板塊南部，東與郯廬深大斷裂相距較近，屬于嵩箕構(gòu)造區(qū)東緣永夏斷隆帶。本區(qū)屬華北地層區(qū)、魯西分區(qū)、徐州小區(qū)。缺失上奧陶統(tǒng)至下石炭統(tǒng)和三疊系至古新統(tǒng)等地層，區(qū)域內(nèi)除靈壁、泗縣、濉溪、渦陽的石弓、龍山等地有震旦、寒武、奧陶系地層及烈山、白土有石炭、二疊系地層出露外，地表被新近系和第四系松散沉積物覆蓋。地層由老至新有青白口系、震旦系、寒武系、奧陶系、石炭系、二疊系、侏羅系、白堊系、古近系、新近系及第四系。含水層從上到下有第四系、第三系沖積含水層組，煤系地層砂巖含水層組和灰?guī)r巖溶裂隙含水層組。其中孔隙含水層組共有4個(gè)含水層和3個(gè)隔水層。據(jù)抽水試驗(yàn)資料，主采煤層頂?shù)装迳皫r裂隙含水層q=0.004 012 ～0.009 37 L/s·m。井下揭露的突水點(diǎn)變化規(guī)律，一般是開始涌水量較大，隨時(shí)間延長，衰減較快，呈淋水或滴水狀態(tài)。總的來說煤系砂巖裂隙含水層富水性較弱，因此，在穿過堅(jiān)硬砂巖時(shí)，必須注意儲(chǔ)存水量突然潰出。

二多元統(tǒng)計(jì)模型

收集了本礦井各地層地下水樣品常規(guī)水質(zhì)分析數(shù)據(jù)62組，水樣采集于水源井、抽水孔及采掘工作面，時(shí)間主要是2012－2013年，此期間礦區(qū)總體開采規(guī)模不大，地下各含水層未出現(xiàn)較大水位變化，故本次研究從靜態(tài)角度對(duì)水質(zhì)資料進(jìn)行多元統(tǒng)計(jì)分析。

2.1 地下水化學(xué)的主成分分析

地下水重分布最廣，含量較高的7種離子:Na+、K+、Ca2+、Mg2+、Cl－、SO42－、CO32－，它們的不同濃度，組成反映了不同的地下水水質(zhì)特征，其濃度差異是地下水化學(xué)分類的主要依據(jù)。這些離子并非完全獨(dú)立地存在于地下水體統(tǒng)中，它們之間總有些相關(guān)性，離子之間的相關(guān)意味著其反映的地址信息重合。若不作處理，就進(jìn)行聚類分析，相當(dāng)于給這些離子變量進(jìn)行了加權(quán)，夸大了某些地址信息因素的作用，使計(jì)算結(jié)果發(fā)生畸變。為克服這種現(xiàn)象，本文先對(duì)離子進(jìn)行主成分分析，找出基本的數(shù)據(jù)結(jié)構(gòu)，用新生成的主成分得分作為新變量進(jìn)行聚類分析。利用spss軟件的Na+、K+、Ca2+、Mg2+、Cl－、SO42－和CO32－的相關(guān)系數(shù)矩陣(表1)。由表1得出:Ca2+和 Mg2+、Na+、K+與Cl－直接存在較強(qiáng)的相關(guān)關(guān)系，其他離子變量直接也存在一定相關(guān)關(guān)系，需要對(duì)原始變量數(shù)據(jù)進(jìn)行因子分析。

利用因子分析中的主成分分析法，從離子相關(guān)系數(shù)矩陣中提取初始因子，得到了相關(guān)系數(shù)矩陣的6個(gè)特征值(主成分所解釋的方差)。按照提取主成分的個(gè)數(shù)一般要求其累計(jì)方差貢獻(xiàn)超80%的原則，提取前3個(gè)主成分，其累計(jì)方差貢獻(xiàn)率達(dá)89.661%(見表2)。

選取的因子初始載荷矩陣一般需要進(jìn)行方差極大正交旋轉(zhuǎn)，使各離子變量在同一因子上的載荷區(qū)別明顯。表3為采用方差極大正交學(xué)專方法迭代3次后的因子載荷矩陣，每個(gè)離子變量的信息被保留的程度體現(xiàn)在各自的公因子方差上。

表1相關(guān)系數(shù)矩陣Tab.1 Correlation matrix

表2解釋的總方差Tab.2 Total variance analysis

表3旋轉(zhuǎn)因子載荷矩陣Tab.3 Rotated component matrix

表4主成分得分Tab.4 The score of main component score

2.2 水質(zhì)樣品的系統(tǒng)聚類分析

聚類分析是研究事物分類的一種方法，將一批樣本或變量按照他們?cè)谛再|(zhì)上的親疏程度加以分類。為獲取礦區(qū)地下水各子系統(tǒng)典型水樣，利用系統(tǒng)聚類法中的最短距離法、最長距離法、類平均法、重心法和離差平方和法等5種方法對(duì)水樣進(jìn)行聚類，距離測(cè)度使用平方歐式距離，聚類變量即為正交因子(TDS因子和堿度因子)。在上述主成分分析的基礎(chǔ)上，利用旋轉(zhuǎn)后的因子回歸計(jì)量值，即主成分1、主成分2和主成分3的得分記為f1、f2、f3，代替原始水樣數(shù)據(jù)用于聚類分析進(jìn)行系統(tǒng)聚類分析。用歐式距離表示主成分得分值之間的相似性，即

其中，i、k均為樣品序號(hào)，i=1，2，3，…62;j為樣品指標(biāo)數(shù)，j=1，2，3。主成分得分值見表4。

利用spss軟件對(duì)62個(gè)水樣的主成分得分進(jìn)行了聚類分析，聚類結(jié)果如圖2所示。若將并類距離d設(shè)為5，則可以將62組水樣分為4類，即1，2，3為一類，4，5，6分別自成一類。如圖兩個(gè)橢圓區(qū)域所示，這兩個(gè)橢圓區(qū)域中，8101工作面的水樣為本礦10煤和太灰。根據(jù)這聚類結(jié)果可判定，臥龍湖8101面部分出水與本礦10煤和太灰水質(zhì)相似性高，地下水水質(zhì)類型較為相似。

2.3 地下水的Bayes判別分析

利用系統(tǒng)聚類得到的地下水各含水層的典型水樣，先剔除水樣類別錯(cuò)誤的，選擇剩下的44組水樣，建立判別函數(shù)。Bayes線性判別一般要求母體變量觀測(cè)值必須服從多元正態(tài)分布，而且要求各母體的協(xié)方差陣無顯著差別。在應(yīng)用中，當(dāng)各母體的樣本規(guī)模比較接近時(shí)，如最大一組的樣本容量不超過最小一組樣本容量的1.5倍時(shí)，違反協(xié)方差陣相等的假設(shè)條件，影響也不太大。表5為函數(shù)系數(shù)。

其中y6、y7分別表示礦區(qū)第四系和煤系水的判別函數(shù)。判斷結(jié)果與聚類結(jié)果一致，接近100%的準(zhǔn)確率，說明建立的判別函數(shù)能很好的預(yù)測(cè)突水水樣的類別。判別函數(shù)代入 Na+、K+、Ca2+、Ma2+、Cl－、SO42－、CO32－的濃度即可。功能已足夠趕超任何其他專用的統(tǒng)計(jì)軟件，Matlab統(tǒng)計(jì)工具箱幾乎包括了數(shù)理統(tǒng)計(jì)方面的書有概念、理論、方法、算法以及實(shí)現(xiàn)。在應(yīng)用上，Matlab具有其他軟件不可比擬的操作簡(jiǎn)單，接口方便，擴(kuò)充能力強(qiáng)等優(yōu)勢(shì)。因此可以預(yù)見該軟件在統(tǒng)計(jì)應(yīng)用上的重要地位，利用臥龍湖礦及周邊礦的水樣以實(shí)例給出Matlab在主成分分析、聚類分析、判別分析的應(yīng)用。

3 Matlab軟件在多元統(tǒng)計(jì)分析的應(yīng)用

在Matlab較早的版本中，統(tǒng)計(jì)功能還不是那么強(qiáng)大，而在Matlab6.x版本中，僅在統(tǒng)計(jì)工具箱(Statitic Toolbox)中的功能函數(shù)就達(dá)到200多個(gè)，

表5函數(shù)系數(shù)Tab.5 The function coefficient

3.1 主成分分析

利用Matlab6.5中的princomp命令實(shí)現(xiàn)，具體程序如下:

x=［9.432 20.265 13.177 17.836 14.179 11.015 45.792 43.938 24.186 8.055 8.251 14.362 40.611 21.273 23.434 23.179 35.058 23.147 18.077 28.821 23.625 36.627 12.464 19.793 37.334 28.993 ….…14.44 27.55 10.19 9.62 9.16 11.79 10.27 8.85 6.92 6.84 6.62 6.91 8.98 9.10 9.06 6.76 7.22 5.99 7.29 7.07 6.55 5.70 6.31 2.66 20.57 12.06 17.32 4.72 10.15 3.05 15.48 1.66 2.75 2.51］

stdr=std(x);%求各變量標(biāo)準(zhǔn)差

［n，m］=size(x);

sddata=x./stdr(ones(n，1)，:);% 標(biāo)準(zhǔn)化變換

［p，princ，egenvalue］=princomp(sddata)%調(diào)用主成分分析程序

p3=p(:，1:3)%輸出前三個(gè)主成分系數(shù)

sc=princ(:，1:3)%輸出前三個(gè)主成分得分

egenvalue%輸出特征根

per=100*egnvalue/sum(egenvalue)%輸出各個(gè)主成分貢獻(xiàn)率

執(zhí)行后得到所要的結(jié)果，這里事前三個(gè)主成分、主成分得分、特征根，即+0.217 0 ×4+0.640 3 ×5 －0.026 9 ×6;Z3= －0.176 6 ×1 －0.213 0 ×2 －0.190 1 ×3 －0.566 0×4+0.304 5 ×5 －0.688 7 ×6;

第一主成分貢獻(xiàn)率為50.70%，第二主成分貢獻(xiàn)率23.23%。第三主成分貢獻(xiàn)率為15.73%，前三個(gè)主成分累計(jì)貢獻(xiàn)率達(dá)89.66%。與spss軟件運(yùn)算得到方差貢獻(xiàn)率(89.661%)基本一致。

3.2 聚類分析

利用Matlab6.5中的cluster命令實(shí)現(xiàn)，具體程序如下

y=pdist(xx);%計(jì)算個(gè)樣本間距離(這里為歐氏距離)

z=linkage(y);%進(jìn)行聚類(這里為最短距離法)

t=cluster(z，4)%將全部樣本分為4類

find(t==2);%找出屬于第2類的樣品編號(hào)h=dendrogram(z);%畫聚類譜系圖

執(zhí)行后得到所要結(jié)果，聚類譜系圖見圖3。

4 結(jié)論

1)在保證原始數(shù)據(jù)信息損失最小的前提下，通過主成分分析法，前3個(gè)主成分(50.70%，23.23%，15.73%)，其 累 計(jì) 方 差 貢 獻(xiàn) 率 達(dá)89.661%。

2)聚類分析方法，得出8101工作涌水水源與本礦10煤砂巖水質(zhì)類型一致。

3)因子分析方法，得出臥龍湖礦各地層地下水水質(zhì)的首要因子為 Ca2+、Mg2+、SO42－，體現(xiàn)太灰水的外來補(bǔ)給;第二因子是陽離子交替吸附作用，其涌水體現(xiàn)了地下水靜儲(chǔ)量消耗。

4)建立的Bayes線性判別模型，具有計(jì)算簡(jiǎn)便、誤判率低、穩(wěn)定性高等特點(diǎn)。同時(shí)運(yùn)用Maltab程序，解釋的結(jié)果與多元統(tǒng)計(jì)模型研究結(jié)果總體上一致，且與該礦區(qū)水文地質(zhì)條件相吻合。

模型的建立依賴于水樣的數(shù)量、變量的選擇、總?cè)芙夤腆w物、硬度、堿度等，建議今后采集更多的地下水樣品，從中選擇時(shí)間跨度較小且空間分布均勻的水樣作為樣本，模型需不斷調(diào)整和完善。

［1］尹國勛，楊 娜，賀玉曉，等.焦作市市區(qū)地下水水質(zhì)現(xiàn)狀評(píng)價(jià)［J］.環(huán)境工程，2004，22(4):66 －69.

［2］賁旭東，郭黃海，解奕偉，等.模糊綜合評(píng)判在礦井突水水源判別中的應(yīng)用及探討［J］.礦業(yè)安全與環(huán)保，2006，33(3):57 － 59.

［3］王廣才，王秀輝，李競(jìng)生，等.平頂山礦區(qū)礦井突(涌)水水源判別模式［J］.煤田地質(zhì)與勘探，1998，26(3):47 － 50.

［4］閆志剛，杜培軍，郭達(dá)志.礦井涌水水源分析的支持向量機(jī)模型［J］.煤炭學(xué)報(bào)，2007，32(8):842 －847.

［5］陳桂明.MATLAB數(shù)理統(tǒng)計(jì)(6.x)［M］.北京:科學(xué)出版社，2002.

［6］劉則毅.科學(xué)計(jì)算技術(shù)與Matlab［M］.北京:科學(xué)出版社，2001.

［7］李濤.Matlab工具箱應(yīng)用指南一應(yīng)用數(shù)學(xué)篇［M］.北京:電子工業(yè)出版，2000.