基于改進AHP法和CRITIC法耦合賦權的松散承壓含水層富水性評價

葛如濤, 陳陸望, 王迎新, 張 杰, 李蕊瑞

(合肥工業大學 資源與環境工程學院,安徽 合肥 230009)

我國華北隱伏型煤田第四系松散承壓含水層覆蓋于煤系地層之上,是煤礦頂板水害防治、地下水資源管理與生態環境保護等研究的熱點。該類含水層以非膠結砂土、砂礫為骨架,砂泥互層明顯,具有承壓性,富水性空間展布差異性較大;此外,受煤礦區長期采動影響,含水層富水性的動態變化明顯[1-2]。目前,對于松散承壓含水層富水性評價與分區,研究方法可分為以下2類:① 以《煤礦防治水細則》為代表,根據水文地質勘探成果,采用抽水試驗得到的單位涌水量作為標準評價指標進行含水層富水性評價與分區[3];② 通過物探手段進行探測,包括直流電法、瞬變電磁法、高密度電法等[4]。物探手段的解譯成果受研究人員經驗和技術水平影響,存在多解性;采用單位涌水量進行富水性分區,由于抽水試驗成本高、過程復雜,抽水試驗鉆孔數量與抽水試驗次數有限,無法反映采動影響下含水層富水性的時空動態變化。

為了提高松散承壓含水層富水性評價與分區的準確性與普適性,文獻[5-6]借助地理信息系統(geographic information system,GIS)的空間信息融合功能,運用沉積控水規律建立多因素復合評價模型;文獻[7]基于層次分析(Analytic Hierarchy Process,AHP)法,借助GIS空間分析功能,開展多元信息融合,對松散承壓含水層的富水性進行分析與評價。采用多因素復合法,關鍵在于如何合理確定各影響因素的權重。采用AHP法或者改進AHP法對主控因素進行賦權[8-9],屬于主觀賦權法,主觀性較強,存在局限性。文獻[10]在研究含水層非均質性問題時采用熵權法對主控因素進行賦權,但該方法在賦權時未考慮專家意見;文獻[11]在研究風化基巖富水性評價方法時,分別使用AHP法、熵權法及耦合兩者的方法來確定影響因素權重,結果表明,耦合AHP法和熵權法確定權重進而預測的富水性類別,與實測數據吻合度明顯高于單獨使用某一賦權法。而熵權法僅根據各影響因素本身變異大小確定其權重,未考慮影響因素之間的沖突性;基于指標相關性的指標權重確定(CRiteria Importance Through Intercriteria Correlation,CRITIC)法兼顧影響因素變異大小和沖突性,在賦權準確程度方面優于熵權法[12]。因此,耦合改進AHP法與CRITIC法的賦權方法會使評價結果更為合理。

本文以淮北煤田松散承壓含水層為研究區,將改進AHP法和CRITIC法進行耦合,來確定第四系松散承壓含水層富水性影響因素的權重,提出一種更科學的、多因素復合的松散承壓含水層富水性評價模型。該模型所需參數均可通過普通地質勘探孔數據獲得,無需進行更多的抽水試驗,降低了抽水鉆孔的施工成本和時間投入;相比于僅使用礦區為數不多的抽水試驗孔數據進行含水層富水性評價與分區,該模型評價準確度更高。此外,由于含水層水位受采動影響,發生永久性或半永久性變化,以致含水層富水性產生相應變化,該富水性評價模型能夠較準確地對研究區富水性進行時空動態評價。

1 研究區地質背景

淮北煤田地處華北板塊,位于華北平原南部、郯廬斷裂帶西側;區內除靈璧、泗縣、濉溪、渦陽等地有新元古界和下古生界基巖出露外,其余絕大部分地區為厚度180～400 m的第四系松散層所覆蓋[13-14];煤田所在區域屬于暖溫帶半濕潤性氣候,年均降水量為700～950 mm,年降雨量最大為1 107 mm;第四系松散層除直接接受大氣降水補給外,還受到地表水體的滲流補給;煤田地表水系發育,沱河、淝河、新汴河、澮河、濉河等河流呈不對稱羽狀展布于淮河北岸,對淺部松散層起到豐蓄枯補的作用[15]?；幢泵禾锼缮幼陨隙驴蓜澐譃?個含水層(組)及3個隔水層(組),第一含水層(簡稱 “一含”)為近地表的潛水孔隙含水層,第二、三、四含水層(簡稱“二含”、“三含”、“四含”)為孔隙承壓含水層[16]。其中,四含直接覆蓋在煤系地層之上,是礦井開采的直接充水水源。含水層多為黏土質砂、粉細砂和中砂交替沉積,一般單層厚度較小,多為0.3～3.5 m。含水層內含有黏土夾層。隔水層多由鈣質黏土與砂質黏土組成,第一、二隔水層厚度較小,普遍小于25 m,隔水性能一般;第三隔水層厚度大,平均厚度達80 m,隔水性能強。

2 富水性評價模型的建立

2.1 影響因素確定

影響松散承壓含水層富水性的因素是多元的。本文通過收集淮北煤田青東煤礦、朱仙莊煤礦、祁南煤礦、祁東煤礦相關水文地質資料,分析上述4個煤礦松散承壓含水層抽水試驗60個有效鉆孔的數據,發現其富水性主要受補給情況、含水層厚度、巖性及砂泥分布情況控制。其中,巖性為定性指標,故先按類別對其量化,然后用級配系數來定量表示巖性特征。最終選取隔水系數R、水頭系數G、含水層厚度S、最厚砂層厚度M、級配系數D、砂泥互層系數P共6個影響因素。

(1) 隔水系數R。松散承壓含水層富水性與地表水體、大氣降水直接或間接補給有關,含水層埋深越淺,接受的補給越充分。隔水層的存在使得隔水層下伏含水層接受的補給受到抑制,其抑制程度與上覆隔水層厚度相關。將松散層某一層位含水層的上覆隔水層累加厚度與該層位含水層底板埋深的比值定義為隔水系數R,即

(1)

其中:r為上覆隔水層總厚度;H為該含水層底板埋深。R越大,該含水層接受補給越困難。

(2)水頭系數G。承壓水頭的變化能夠反映含水層承壓水位的變化,將某一層位含水層承壓水頭T與該含水層頂板埋深H0的比值定義為水頭系數G,即

(2)

G越大,說明單位厚度的含水層給水能力越強,與地表的水力聯系越強。

(3) 含水層厚度S。S是影響含水層富水性強弱的重要因素,地下水賦存情況與S密切相關,在其他影響因素相差不大的情況下,S越大,含水層富水性越強。

(4) 最厚砂層厚度M。M是指某一含水層內厚度最大的一層砂礫層的厚度,只要沒有黏土夾層,可以是不同粒徑砂礫層累加厚度。在S一定的條件下,M越大,說明層內黏土層越薄、層數越少,則含水層儲水空間越大,層內水力聯系越強,富水性越強。在S及其他影響因素相差不大的情況下,M與富水性強弱相關。

(5) 級配系數D。在其他影響因素一定的情況下,松散含水層土體粒徑組合不同,其富水性不同。根據抽水試驗成果可知,卵礫石層較砂礫層富水性強,粒徑大的砂礫層較粒徑小的富水性強。將土體粒徑組合對含水層富水性的貢獻大小進行量化,見表1所列,進而構建表征含水層富水性強弱的粒徑組合系數,即級配系數D,D值越大,富水性越強。

D計算公式為:

D=∑dγcγ/S

(3)

其中:dγ為粒徑類別的量化值;cγ為不同粒徑土層的厚度;γ為土體粒徑類別。

表1 松散層按土體粒徑級別賦值結果

(6) 砂泥互層系數P。P是在砂泥比基礎上進行改進的評價指標。以往研究采用的砂泥比僅簡單計算含水層內砂礫層與黏土層累計厚度的比值,而未考慮垂向砂泥互層分布對含水層富水性的影響。鉆探成果顯示,淮北煤田內各級松散承壓含水層在垂向上一般都為多層結構,具體表現為砂礫層與黏土層交互沉積,故將這樣的一組上覆砂礫層及其下伏黏土層定義為一個砂泥互層。在含水層厚度相近時,互層層數越多,單個互層內黏土比率越大,意味著層內砂礫層之間的水力聯系越差,富水性相應較弱;反之,富水性較強。

朱仙莊煤礦8105-四含檢2孔與青東煤礦2015-水1孔2個抽水鉆孔揭露的四含厚度與砂泥比均相差無幾。前者有2個互層,作為隔水層的黏土層將砂層分割開,砂層之間水力聯系受阻;而后者僅有1個互層,為單一砂層含水層。2個鉆孔抽水試驗得到的單位涌水量q分別為2.8×10-5、0.01 L/(s·m),驗證了上述觀點。依據各個互層中黏土厚度比率對富水效果的貢獻大小將其量化,構建表征砂泥互層影響富水性強弱的砂泥互層系數P,P值越大,層內水力聯系越弱,富水性越弱。P計算公式為:

P=∑pεlε/S

(4)

其中:pε為互層內黏土厚度所占比率的賦值;lε為各個互層的厚度;ε為各互層序號。

pε賦值見表2所列。

表2 單個互層內黏土厚度比率級別賦值結果

淮北煤田青東、朱仙莊、祁南、祁東4個煤礦松散承壓含水層60個抽水試驗有效鉆孔及其層位見表3所列,對應60個鉆孔的6個影響因素取值見表4所列。

表3 淮北煤田抽水試驗60個鉆孔名稱及其層位

2.2 權重確定

2.2.1 改進AHP法的賦權方法

傳統的AHP法在構造判斷矩陣時,各指標之間相對重要程度的判斷會受到該領域專家研究方向、工作經驗等影響,具有一定的主觀性;此外,該方法固化了標度,不具有擴展性和傳遞性,導致矩陣表述與實際情況相差很大[17]。改進AHP法的賦權方法如下:

(1) 采用極大值法和極小值法對影響因素進行歸一化處理,消除影響因素單位不同對數據比較造成的影響。極大值法適用于與評價結果成正比的影響因素數據歸一化,極小值法適用于與評價結果成反比的影響因素數據歸一化。極大與極小值法計算公式分別為:

(5)

(6)

其中:Xi為第i個影響因素數據歸一化后的量化值;i為影響因素序號;maxxi、minxi分別為第i個影響因素數據歸一化前的最大值和最小值。

(2) 運用SPSS軟件中的CORREL函數分別對各影響因素和評價標準值(鉆孔的單位涌水量q)進行相關性分析[18],計算出各影響因素與評價標準值的相關系數,相關系數越大,表明該因素與評價標準值擬合得越好,則該因素相對重要程度越高。采用期望標度法,根據重要程度的傳遞性法則將各因素兩兩比較排序,依次可以計算出判斷矩陣未知元素的值。

期望標度取值見表5所列。

表5 期望標度取值

(3) 計算各影響因素的樣本標準差s(i)(i=1,2,…,n)。樣本標準差越大,說明該因素內部變化形式越豐富,其對于綜合評價的影響越大。因此,可采用各影響因素的樣本標準差來表征各影響因素對富水性的影響程度,從而構建判斷矩陣B1-2。B1-2中元素計算公式為:

(7)

其中:s(i)、s(j)分別為影響因素i和影響因素j的樣本標準差;smax、smin分別為影響因素樣本標準差的最大值和最小值;bmin為相對重要程度參數,bmin=min{9,int[smax/smin+0.5]},int為取整函數。

(4) 由于多階判斷矩陣較為復雜,矩陣中會出現某些元素前后矛盾的現象,有必要檢驗判斷矩陣的一致性,保證輸出結果及最終權重向量的準確性。一致性指標(consistency index,CI)IC檢驗公式為:

IC=(λmax-n)/(n-1)

(8)

IR=IC/RC

(9)

其中:λmax為判斷矩陣的最大特征值;RC為隨機一致性比率(consistency rate,CR);IR為隨機一致性指標(random index,RI)。當CR小于0.1時,可認為判斷矩陣滿足一致性要求,反之,需要調整判斷矩陣。對于六階矩陣,對應的隨機一致性指標RI[19]為1.24。

(5) 求出各個判斷矩陣的最大特征值及其對應的特征向量,標準化處理后可得所求的權重向量W1-1、W1-2。若各判斷矩陣符合一致性檢驗要求,則綜合評價矩陣是可行的。故可取各權重向量組成元素ei的平均值,構建結合矩陣B1-1、B1-2的權重向量W1為:

(10)

改進AHP法具體實現過程如下:

(1) 由2.1節可知,R、P與富水性近似呈負相關性,G、S、M、D與富水性近似呈正相關性。采用對應的歸一化方法((5)式、(6)式)對表4影響因素數據進行歸一化處理。

(2) 分別對各影響因素和評價標準值進行相關性分析,得出各影響因素與評價標準值之間的相關系數,6個影響因素按相關系數從大到小排序依次為:R(0.645)、M(0.492)、P(0.386)、S(0.374)、D(0.344)、G(0.271)。采用期望標度法構建判斷矩陣B1-1為:

B1-1=

(3) 計算得到影響因素樣本標準差,6個影響因素按其樣本標準差從大到小排序依次為:P(0.283)、M(0.253)、S(0.236)、D(0.216)、R(0.188)、G(0.174)。經計算,bmin=min{9,int[2.118]}=2。根據(7)式構建判斷矩陣B1-2為:

B1-2=

(4) 判斷矩陣B1-1、B1-2均為六階矩陣,經計算其最大特征值λmax分別為6.000 1、6.002 2;通過(8)式、(9)式可得CR1-1為2.0×10-5,小于0.1;CR1-2為3.8×10-4,小于0.1。因此,判斷矩陣B1-1、B1-2均滿足一致性檢驗要求。

(5) 求出各判斷矩陣對應最大特征值的特征向量,標準化處理后得到所求權重向量。B1-1對應的影響因素權重W1-1中6個影響因素的權重分別為:R,0.349 4;G,0.071 9;S,0.121 6;M,0.205 5;D,0.093 5;P,0.158 1。B1-2對應的影響因素權重W1-2中6個影響因素的權重分別為:R,0.168 2;G,0.111 7;S,0.227 6;M,0.214 9;D,0.161 3;P,0.116 2。

B1-1、B1-2均滿足一致性檢驗要求,故可用(10)式求得最終R、G、S、M、D、P6個影響因素權重,進而得到權重向量W1,即W1=[0.258 8 0.091 8 0.174 6 0.210 2 0.127 4 0.137 2]。

2.2.2CRITIC賦權方法

CRITIC法是文獻[20]提出的一種基于客觀條件的賦權方法。該方法主要根據各指標的信息量和相關性對其賦權,分別用對比強度和指標沖突性來反映各指標的信息量和相關性。對比強度以各指標內部變異大小來衡量,可以通過標準差來反映變異大小;標準差越大,該指標反映的信息量越大,其權重相應較大。指標沖突性反映指標之間的相關性特征,通過計算各指標相互之間的相關系數來表示;兩個指標的相關系數越大,其正相關性越強,則這兩個指標沖突性較低。CRITIC法不僅考慮影響因素變異程度對權重的影響,還考慮各因素之間的沖突性特征,而熵權法僅考慮到前者,因此CRITIC法在綜合評價效果上要優于熵權法[12]。因素兩兩之間的Pearson相關系數ρ、各因素的標準差σ計算公式分別為:

(11)

(12)

設Ej為第j個影響因素所包含的信息量,Ej越大,該因素相對重要程度越高,由此可以計算出第j個因素的客觀權重wj,進而可求出權重向量W2=[w1w2…wn]。Ej、wj計算公式為:

(13)

(14)

CRITIC法具體實現過程如下:

(1) 與改進AHP法相同,首先用(5)式、(6)式對影響因素的數據進行歸一化處理,消除由于影響因素單位不同造成的影響,然后用(12)式計算得到6個影響因素內部的標準差分別為:P,0.283 2;M,0.253 3;S,0.236 2;D,0.216 4;G,0.187 6;R,0.173 6。

(2) 利用SPSS軟件,通過(11)式計算得到影響因素兩兩之間的Pearson相關系數,見表6所列。

(3) 根據(13)式、(14)式可得R、G、S、M、D、P6個影響因素權重,進而得到權重向量W2,即

W2=[0.171 3 0.187 2 0.157 3

0.125 9 0.171 7 0.186 5]。

表6 影響因素之間的Pearson相關系數

2.2.3 耦合賦權

為使建立的耦合賦權模型兼顧該領域專家的經驗和實測數據的客觀性特征,在改進AHP法賦權與CRITIC法賦權求出主、客觀權重向量的基礎上,基于總偏差最小化的原則,采用博弈論對改進AHP法和CRITIC法進行耦合,從而得到耦合賦權模型。

(1) 采用2種方法對各影響因素進行賦權,得到2組不同的權重數據。2個評價結果的隨機性組合為:

(15)

其中:W為所有可能出現的權重數據組合向量;α1、α2為線性組合系數;W1、W2為采用某種賦權方法所獲得的權重數據組合向量。

(2) 基于博弈論的權重集結本質,對線性組合系數α1、α2優化處理,據此建立決策模型。

(3) 根據矩陣的微分性質,計算求得(15)式的最優化一階導數條件矩陣,即

(16)

(17)

(18)

(5) 獲取最優化系數后,建立最優化權重耦合模型為:

(19)

2.2.4 含水層富水性評價模型的建立

多因素復合的松散承壓含水層富水性評價流程如圖1所示。

圖1 松散承壓含水層富水性評價的流程圖

基于改進AHP法和CRITIC法耦合得到最終的權重向量W*,根據線性加權的方法建立淮北煤田多因素復合的松散承壓含水層富水性評價模型為:

0.170 5f(S)+0.190 1f(M)+

0.138 0f(D)+0.148 9f(P)

(20)

其中:V為富水性綜合指數,取值范圍為0～1;Ai為各影響因素數據歸一化后的數值;f(R)、f(G)、f(S)、f(M) 、f(D)、f(P)分別為R、G、S、M、D、P6個影響因素數據歸一化后的數值。

聚類分析法能夠從樣本數據入手,在保證組間差異最大化、組內差異最小化的基礎上恰當有效地進行分類。利用SPSS軟件中的系統聚類分析法,對60個鉆孔的富水性綜合指數分級分類。由于淮北煤田松散層未發現極強富水性的抽水鉆孔資料,將計算得到的富水性綜合指數進行3級劃分,結果如圖2所示。

圖2 淮北煤田松散承壓含水層富水性綜合指數聚類分析結果

通過歐氏距離7.5的點縱穿樹狀圖畫同型線,Ⅰ區對應弱富水性,Ⅱ區對應中等富水性,Ⅲ區對應強富水性。Ⅰ區、Ⅱ區與Ⅲ區之間的中斷值見表7所列。

表7 淮北煤田松散承壓含水層富水性評價分區中斷值

表7中,遵循《煤礦防治水細則》[3]中的含水層富水性4級劃分原則,選取研究區松散承壓含水層抽水試驗鉆孔對應的富水性綜合指數最大值作為劃分強富水(Ⅲ區)和極強富水(Ⅳ區)的中斷值。

3 富水性評價模型的驗證

為了驗證本文建立的多因素復合的松散承壓含水層富水性評價模型的準確性,收集淮北煤田錢營孜煤礦、任樓煤礦、孫疃煤礦及臥龍湖煤礦相關數據,整理得到松散層有效抽水試驗鉆孔11個,采用11個鉆孔的數據進行模型驗證。首先根據(20)式計算11個鉆孔數據的富水性綜合指數,并根據中斷值的大小進行分類,然后將該模型得到的含水層富水性與由鉆孔抽水試驗得到的富水性進行對比,驗證結果見表8所列。11個鉆孔數據中,8個弱富水層位均評價正確;對于3個中等富水層位,2個評價正確;評價不正確的鉆孔為一含抽水試驗孔。由于一含為近地表潛水或半承壓含水層,其性質與承壓含水層有所差別,因此,本文模型驗證正確率達90.91%。表8中,單位涌水量q的單位為L/(s·m)。

表8 單位涌水量q富水性與模型評價結果對照

4 工程應用

4.1 資料收集與處理

選擇朱仙莊煤礦北部采區,對其松散承壓含水層四含進行富水性評價與分區。將63個普通地質鉆孔柱狀圖與相關水文資料,按照不同影響因素整理出各自數據并歸一化處理后,通過多因素復合的淮北煤田松散承壓含水層富水性評價模型,計算出對應各個鉆孔的富水性綜合指數V。

4.2 富水性評價與動態預測

朱仙莊煤礦北部采區受采動影響,四含水位明顯降低,其水位埋深及變動情況見表9所列。由表9可知,從2017年12月至2020年7月,四含水位平均下降36.75 m。

依據該含水層富水性綜合指數V與中斷值,得到朱仙莊煤礦北部采區四含富水性分區,如圖3所示。對比圖3a、圖3b可以看出,朱仙莊煤礦北部采區受采動影響,中等富水區明顯減小。根據抽水鉆孔I-I-6在1964年7月的抽水試驗,單位涌水量q為0.136 L/(s·m),該數據顯示該區域為中等富水,根據2017年水位進行評價,該區域仍為中等富水,但根據2020年水位進行評價,該區域已變為弱富水。

表9 朱仙莊煤礦北部采區四含水位埋深及變動情況

圖3 朱仙莊煤礦北部采區不同水位埋深下的四含富水性分區結果

5 結 論

(1) 選取隔水系數R、水頭系數G、含水層厚度S、最厚砂層厚度M、級配系數D、砂泥互層系數P共6個影響因素作為評價指標,基于博弈論對改進AHP法和CRITIC法進行耦合,建立多因素復合的松散承壓含水層富水性評價模型,建立的耦合賦權模型具有兼顧該領域專家經驗和實測數據的客觀性特征。

(2) 對該多因素復合的松散承壓含水層富水性評價模型進行驗證,評價結果正確率高達90.91%,所選評價指標均可通過普通地質勘探孔數據得到,無需過多開展現場抽水試驗與相關物探工作,解決了由于抽水孔數量少導致分區結果不理想的問題,也降低了抽水孔的施工成本,避免了運用物探手段進行富水性評價時,解譯成果受研究人員的經驗和技術水平影響、有主觀局限性、存在多解性等問題。

(3) 松散承壓含水層的水位等水文地質條件受采動影響,會發生永久性或半永久性變化,以致含水層富水性也相應發生改變。通過該多因素復合的松散承壓含水層富水性評價模型,可以對研究區富水性開展時空動態評價與預測。