王甜甜，靳德武，劉 基，楊 建，王心義，趙 偉

(1.煤炭科學(xué)研究總院，北京 100013; 2.中煤科工集團(tuán)西安研究院有限公司，陜西 西安 710077; 3.陜西省煤礦水害防治技術(shù)重點實驗室，陜西 西安 710077; 4.河南理工大學(xué) 資源環(huán)境學(xué)院，河南 焦作 454000)

煤炭資源安全高效開采是關(guān)乎我國經(jīng)濟和發(fā)展的大事，然而，我國煤礦水文地質(zhì)條件復(fù)雜，煤層開采時礦井水害頻發(fā)[1]。據(jù)統(tǒng)計從2000年至2014年，4 500人死于礦井突水事故，造成了巨大的人員傷亡與經(jīng)濟損失[2-3]。礦井一旦發(fā)生突水，快速準(zhǔn)確地判斷突水水源，是治理水害的關(guān)鍵[4]。

礦井突水水源識別通常使用地下水化學(xué)、同位素、水溫、地下水位動態(tài)觀測等方法,地下水化學(xué)因能反應(yīng)地下水的本質(zhì)特征，且能準(zhǔn)確、快速、經(jīng)濟地進(jìn)行水源識別，在礦井突水水源識別中得以廣泛應(yīng)用[5-6]。從20世紀(jì)90年代起眾多學(xué)者就開始利用水化學(xué)特征進(jìn)行突水水源判識，孫福元[7]、陳朝陽[8]、王廣才等[9]利用礦井水質(zhì)特征及變化規(guī)律與目標(biāo)水源進(jìn)行對比，成功地進(jìn)行了突水水源判識，為水化學(xué)識別突水水源奠定了深厚的基礎(chǔ)。21世紀(jì)以來，隨著基礎(chǔ)理論以及計算機科學(xué)的發(fā)展，逐步發(fā)展出基于數(shù)學(xué)方法(多元統(tǒng)計學(xué)方法、灰色系統(tǒng)法、模糊數(shù)學(xué)法)的判別模型和結(jié)合計算機技術(shù)的其他判別法(BP神經(jīng)網(wǎng)絡(luò)法、GIS法、SVM 法、可拓識別法)[10-14]。王心義[15]、周健[16]等使用距離判別法識別了礦井突水水源;武強等[17]借助層次分析法進(jìn)行煤礦突水水源判識;孫亞軍等[18]借助計算機，基于GIS識別了突水水源。這些方法提高了礦井突水水源的準(zhǔn)確性，但是仍存在一定的缺陷，例如距離判別法將所有識別因子同等看待，擴大或縮小了不同識別因子對整體的影響，層次分析法權(quán)重計算一定程度受人為因素的影響，且需要一致性校驗，計算過程復(fù)雜繁瑣;GIS法依據(jù)模糊綜合理論，其評價因子權(quán)重及最終隸屬度難以確定，直接影響水源識別的準(zhǔn)確性?；疑到y(tǒng)、模糊數(shù)學(xué)法等均存在以上類似問題或缺陷。

1 動態(tài)權(quán)-集對分析模型

1.1 集對分析原理與改進(jìn)

集對分析是處理不確定性問題的系統(tǒng)分析方法，其核心思想是將事物中客觀存在的不確定性，以辯證分析(同、異、反)表示，即以某種聯(lián)系度來描述事物的不確定性[19]。假設(shè)這種聯(lián)系度用μ表示，A,B兩集合構(gòu)成集對D，即D=(A,B)，若集對D中有N個特征，其中S個特征共屬于A，B集合，P個特征分別屬于兩者，其余F個特征不屬于兩者。具體聯(lián)系度μ可通過數(shù)學(xué)表達(dá)式表述，即其表達(dá)式為

(1)

其中，a=S/N為同一度，b=F/N為差異度，c=P/N為對立度，并且N=S+F+P，同時a,b,c∈[0,1]，a,c相對確定，b相對不確定，a+b+c=1;i為差異標(biāo)記符;j為對立標(biāo)記符。

單純地從同、異、反3個側(cè)面刻畫不確定性顯得過于粗糙，因此，可將聯(lián)系度進(jìn)一步細(xì)化，將式(1)改寫為

μ=a+(b1-+b1+)i+(b2-+b2+)i+…+

(bn-+bn+)i+(c1-+c1+)j+(c2-+

c2+)j+…+(cn-+cn+)j=a+b1-i-+

b1+i++b2-i-+b2+i++…+bn-i-+bn+i++

c1-j-+c1+j++c2-j-+c2+j++…+cn-j-+cn+j+

(2)

式中,i-,j-為i,j的反方向;i+,j+為i,j的正方向。

當(dāng)n=1時得到五元聯(lián)系度公式:

μ=a+(b1-+b1+)i+(c1-+c1+)j=

a+b1-i-+b1+i++c1-j-+c1+j+

(3)

式中，a+b1-+b1++c1-+c1+=1，對于突水水源識別，假設(shè)b1-,b1+為與所屬水源類型相鄰的鄰左區(qū)和鄰右區(qū)系數(shù)，c1-,c1+為與所屬水源類型次相鄰的次鄰左區(qū)和次鄰右區(qū)系數(shù)。

如圖1所示，屬區(qū)、鄰區(qū)和次鄰區(qū)在整個區(qū)間各占1/3，鄰左、次鄰左、鄰右、次鄰右各占當(dāng)1/6，即Y4-Y1=2(Y3-Y2);識別處于所屬水源類型中，則a=1，其他為0;當(dāng)在鄰左區(qū)時，則越靠近屬區(qū)，a越大，b1-越小，反之，越遠(yuǎn)離屬區(qū)，a越小，b1-越大;當(dāng)在鄰右區(qū)時，則越靠近屬區(qū)，a越大，b1+越小，反之，越遠(yuǎn)離屬區(qū)，a越小，b1+越大;當(dāng)在次鄰左區(qū)時，則越靠近屬區(qū)，則a,b1-越大，c1-越小;當(dāng)在次鄰右區(qū)時，越靠近屬區(qū)，則a,b1+越大，c1+越小。

圖1 改進(jìn)集對分析分區(qū)解釋Fig.1 Zoning interpretation map of improved set pair

聯(lián)系度μ的確定，將各類水源各指標(biāo)按上述分5個區(qū)域，設(shè)X∈[Y2,Y3)，X為待測數(shù)據(jù)，Y2，Y3為所屬區(qū)水源類型的下限和上限，上文介紹，當(dāng)識別水源處于所屬區(qū)中，則a=1，其他為0，所以其聯(lián)系度μ=1。當(dāng)X?[Y2,Y3]即處于次鄰左、鄰左、鄰右與次鄰右時的μ[23]為

(4)

1.2 動態(tài)權(quán)重值的確定

動態(tài)權(quán)重是一種將客觀與主觀權(quán)重結(jié)合的一種隨著待評對象不同情況變化的綜合權(quán)重，綜合考慮了主客觀情況，使得權(quán)重的賦值更加合理[24]。客觀權(quán)重采用熵權(quán)法計算，反應(yīng)識別因子的固有性[25]。

令:

(5)

式中，sij為各類水源各指標(biāo)的平均值;vij為各指標(biāo)歸一化值;m為樣本個數(shù);Hi為熵值。

其熵值為

(6)

若式中vij=0，則Hi=0。其熵權(quán)重值為

(7)

主觀權(quán)重通過改進(jìn)層次分析法確定，考慮不同識別因子的實際情況，引入最優(yōu)傳遞矩陣計算主觀權(quán)重，省去一致性校驗[3]。首先對6項識別因子的重要性進(jìn)行比較，建立判斷矩陣A:

其中，如果aij=1表示i比j重要;aij=0表示i和j同樣重要;aij=-1表示j比i重要。

A的最優(yōu)傳遞矩陣為R:

(8)

R的判斷矩陣為D:

式中，

dij=exp(rij)

(9)

利用式(10)，計算各項識別因子的主觀權(quán)重值θi:

(10)

最后確定動態(tài)權(quán)重wi:

(11)

式中，λi為熵權(quán)法計算的客觀權(quán)重;θi為改進(jìn)層次分析計算的主觀權(quán)重。

1.3 集對勢的確定

由式(4)確定聯(lián)系度μ，并組成聯(lián)系度矩陣U，因此可與動態(tài)權(quán)重矩陣W相乘計算綜合聯(lián)系度矩陣Z，其公式為

Z=U·W

(12)

Z反映待評對象中整體的聯(lián)系度，由l，m和n組成。由其可計算聯(lián)系度的集對勢nSPT[26]，其公式為

(13)

2 礦井突水水源識別

2.1 研究區(qū)域概況

葫蘆素井田位于內(nèi)蒙古自治區(qū)東勝煤田呼吉爾特礦區(qū)，面積約為92.76 km2，井田位于鄂爾多斯高原的東北部，東邊與陜西省交界。周圍的中深埋礦區(qū)還包括新街、呼吉爾特、納林河等，如圖2所示。

圖2 研究區(qū)地理位置Fig.2 Location of the study area

圖3 研究區(qū)主要含水層Fig.3 Main aquifers of study area

2.2 水源識別模型建立

根據(jù)葫蘆素井田5種突水水源，將各個指標(biāo)的質(zhì)量濃度繪制成箱型圖(圖4)，圖4中矩形盒上下分別為離子質(zhì)量濃度的上、下四分位數(shù)，中間橫線為中位數(shù)，小方塊為平均值，兩端分別為質(zhì)量濃度最大值和最小值[28]。從箱線圖可以形象直觀地觀察到各個離子在不同含水層中質(zhì)量濃度變化情況。

表1 不同含水層53組水樣水化學(xué)特征
Table 1 Fifty-three groundwater samples from different aquifers

mg/L

續(xù) 表

由圖4對各類水源類型的分析，將圖中的下四分位數(shù)(Q1)和上四分位數(shù)(Q3)設(shè)定為所屬突水水源的限值Y2和Y3。根據(jù)圖1，將Q1-0.5(Q3-Q1)和Q3+0.5(Q3-Q1)設(shè)定為所屬突水水源鄰區(qū)的限值Y1和Y4。

2.3 水源識別模型驗證與分析

突水識別區(qū)間建立后，將建井、開拓、回采3階段取得的已知來源礦井水樣10組(表2)，代入動態(tài)權(quán)-集對分析模型，對模型進(jìn)行驗證、分析。下面以水樣S7為例，詳細(xì)介紹利用動態(tài)權(quán)-集對分析模型進(jìn)行水源類型判識。

圖4 6項識別指標(biāo)的箱型Fig.4 Box-plot of six parameters

mg/L

Xs7=[311.72,0.35,0.62,61.62,430.21,140.09]T

與集對區(qū)間進(jìn)行對比，如將K++Na+=311.72代入第四系突水水源類型中，與K++Na+區(qū)間[10.09,29.20]和[0.54,38.76]對比，其值大于鄰后區(qū)限值38.76，即Xs7(1)∈[38.76,]，因此利用式(4)中次鄰右式計算屬于第四系的聯(lián)系度為μs7(1)=0.06+0.03i++0.90j+，同理可算的各指標(biāo)各突水水源類型的聯(lián)系度矩陣Us7為

表3 識別指標(biāo)動態(tài)權(quán)重
Table 3 Dynamic weights of identification parameters

權(quán)重類型K++Na+Ca2+Mg2+Cl-SO2-4HCO-3客觀權(quán)重0.180.200.100.290.210.01主觀權(quán)重0.190.110.060.150.310.19動態(tài)權(quán)0.200.120.030.250.380.02

利用式(12)將聯(lián)系度矩陣Us7和動態(tài)權(quán)重相乘得到綜合聯(lián)系度矩陣Zs7為

再利用式(13)計算各類突水水源類型的集對勢，nSPT1=0.29，nSPT2=0.68，nSPT3=2.21，nSPT4=7.68和nSPT5=1.23，將集對勢歸一化處理后，并用百分?jǐn)?shù)表示，依次為2.39%，5.59%，18.28%，63.56%和10.18%。根據(jù)集對勢結(jié)果可知直羅組2段含水層的占比為63.56%，可將突水水源判定為直羅組2段。

按照上述方法與步驟，將其余水樣分別代入到動態(tài)集對模型中，其歸一化后的集對勢見表4，并根據(jù)表4繪制出集對勢比例柱狀圖，如圖5所示。

表4 礦井水樣集對勢計算結(jié)果
Table 4 Set pair trend results of ten mine water samples to be tested

水樣序號歸一化集對勢/%第四系白堊系直羅組1段直羅組2段延安組識別結(jié)果實際類型S183.625.684.851.674.18第四系第四系S235.4036.0614.734.319.50白堊系第四系S311.2872.427.273.085.96白堊系白堊系S49.4868.2610.423.468.37白堊系白堊系S50.250.3198.070.291.09直羅組1段直羅組1段S60.050.0499.710.030.17直羅組1段直羅組1段S72.395.5918.2863.5610.18直羅組2段直羅組2段S82.688.0926.1250.6312.49直羅組2段直羅組2段S90.060.030.160.0499.71延安組延安組S100.080.050.330.0699.48延安組延安組

圖5 集對勢比例柱狀Fig.5 Box-plot of six parameters

由集對勢柱狀圖可以明顯識別出所占比例最大的水源類型，即識別突水水源類型。由圖5可知水樣S1第四系集對勢比例最大，為83.03%，其識別水源類型為第四系含水層;水樣S2，水樣S3和水樣S4中白堊系集對勢比例最大，分別為36.06%，72.42%和68.26%，其識別水樣類型為白堊系含水層;水樣S5和水樣S6中直羅組1段集對勢比例最大，分別為98.07%和99.71%，其識別水樣類型為直羅組1段;水樣S7和水樣S8中直羅組2段集對勢比例最大，分別為63.56%和50.63%，其識別水樣類型為直羅組2段含水層;水樣S9和水樣S10中延安組煤系含水層集對勢比例最大，分別為99.71%和99.48%，其識別水樣類型為延安組煤系含水層。

利用動態(tài)權(quán)-集對分析模型的識別結(jié)果與實際的水樣類型進(jìn)行對比，顯然，應(yīng)用動態(tài)權(quán)-集對分析模型識別的10組待測水樣，9組與實際類型相符，僅1組結(jié)果與實際情況不符。水樣S2的水源識別結(jié)果來自于白堊系，但實際來自為第四系，筆者發(fā)現(xiàn)該水樣的識別結(jié)果第四系為35.40%，白堊系為36.06%，兩者相差甚小只有0.66%。因此可用式(13)中l(wèi)，m和n的大小來判斷這兩類水源nSPT的“強弱”。水樣S2的綜合聯(lián)系度矩陣Zs2為

可判斷出第四系和白堊系識別綜合聯(lián)系度中l(wèi)，m和n的比較結(jié)果均為l>n>m即nSPT+(強同勢)。所以其最終的識別結(jié)果仍然來源于白堊系。識別錯誤的原因在于第四系水樣建模數(shù)據(jù)過少，致使模式識別區(qū)間沒有充分考慮到第四系突水水源的特征。

2.4 水源識別模型預(yù)測

2015-04-26，葫蘆素井田21102工作面主回風(fēng)巷施工至3 411 m位置時，距離巷道右側(cè)0.7 m位置中下部出現(xiàn)突水，初始水量為120 m3/h,之后增至180 m3/h，并最終穩(wěn)定至150 m3/h。經(jīng)調(diào)查分析，判斷突水原因為地面“大53”天然氣孔未封閉。為查明氣孔與上下含水層溝通情況，突水過程中連續(xù)采集水樣9組(2015-04-26—05-03)，05-14采集水樣1組，氣孔治理后06-09采集水樣1組，共計11組進(jìn)行水化學(xué)分析(具體數(shù)據(jù)見文獻(xiàn)[27])，利用已建模型，對11個突水水樣進(jìn)行水源預(yù)測，并與實際突水水源對比，結(jié)果見表5。由表5可知，突水初期，突水水源主要是直羅組含水層，先是直羅組1段，然后是直羅組2段;突水中期，水源主要為白堊系含水層;治理過程中的水源為直羅組2段。這均與文獻(xiàn)[27]的研究結(jié)果一致，11組突水水樣識別結(jié)果均與實際情況相符合，并未出現(xiàn)異常。這是由于突水水源主要為直羅組和白堊系含水層，建模過程中直羅組與白堊系樣本量充足，為準(zhǔn)確的突水水源識別奠定了堅實的基礎(chǔ)。

表5 工作面11組突水水樣集對勢計算結(jié)果
Table 5 Set pair trend results of eleven mine water inrush samples from working face

水樣序號取樣日期歸一化集對勢/%第四系白堊系直羅組1段直羅組2段延安組識別結(jié)果實際類型D53-12015-04-265.805.2070.613.1515.24直羅組1段直羅組1段D53-22015-04-277.977.4032.4533.6918.49直羅組2段直羅組2段D53-32015-04-278.045.4330.7437.5018.28直羅組2段直羅組2段D53-42015-04-288.056.3232.2935.0418.29直羅組2段直羅組2段D53-52015-04-298.475.5931.7634.5319.64直羅組2段直羅組2段D53-62015-04-308.324.9931.8834.8819.94直羅組2段直羅組2段D53-72015-05-0111.5228.4623.7617.3918.86白堊系白堊系D53-82015-05-0311.7727.5823.2618.8418.55白堊系白堊系D53-92015-05-030.340.240.570.6898.17延安組延安組D53-102015-05-1410.4725.0822.4021.3620.69白堊系白堊系D53-112015-06-094.124.7322.5828.3540.22直羅組2段直羅組2段

3 結(jié) 論

(1)動態(tài)權(quán)充分考慮了指標(biāo)權(quán)重的客觀固有性和主觀經(jīng)驗性，利用信息熵值計算客觀權(quán)重，使用改進(jìn)層次分析計算主觀權(quán)重，既削弱了人為因素的影響，又考慮了識別指標(biāo)的實際情況，使得權(quán)重賦值合理，有利于突水水源識別。

(3)利用10組礦井水樣對動態(tài)權(quán)-集對分析水源識別模型進(jìn)行驗證，結(jié)果表明，9組水樣的識別結(jié)果與實際情況完全吻合，識別準(zhǔn)確率達(dá)到90%。1組水樣與實際情況不符，為建模時第四系樣本量小，待測水樣超出識別區(qū)間所致 。

(4)利用所建模型對21102工作面11個突水水樣進(jìn)行水源判別，判別結(jié)果均與實際相符合，識別準(zhǔn)確率達(dá)到100%。大量的水質(zhì)樣本為動態(tài)權(quán)-集對分析模型準(zhǔn)確識別突水水源奠定了良好的基礎(chǔ)。