孫 琳,任娜娜,李云安,胡樂健
(中國地質(zhì)大學(xué)(武漢)工程學(xué)院,湖北 武漢 430074)
巖溶即喀斯特(KARST)。地下水和地表水對(duì)可溶性巖石的破壞和改造作用,這種作用及其所產(chǎn)生的地貌現(xiàn)象和水文地質(zhì)現(xiàn)象總稱為巖溶[1]。巖溶本身是一種常見的地質(zhì)現(xiàn)象,但是隨著我國經(jīng)濟(jì)的發(fā)展,基礎(chǔ)工程的建設(shè),巖溶塌陷已經(jīng)成為我國主要地質(zhì)災(zāi)害之一。據(jù)記載,人類活動(dòng)造成的巖溶塌陷存在著逐年遞增的趨勢(shì)。在巖溶地區(qū),巖溶塌陷不明顯,但在人類工程的影響下,存在巖溶塌陷的可能性;因此進(jìn)行巖溶發(fā)育強(qiáng)弱預(yù)測(cè)、并對(duì)重點(diǎn)區(qū)域進(jìn)行巖溶塌陷的排查具有重要的現(xiàn)實(shí)意義。
地質(zhì)災(zāi)害預(yù)測(cè)是地質(zhì)災(zāi)害評(píng)估中重要的一環(huán),一直廣受關(guān)注。國內(nèi)學(xué)者對(duì)巖溶塌陷不確定性預(yù)測(cè)也開始了一些研究。雷明堂等[2]率先將GIS技術(shù)應(yīng)用在巖溶塌陷不確定性評(píng)價(jià)中。隨后一些知識(shí)驅(qū)動(dòng)模型應(yīng)用在巖溶塌陷預(yù)測(cè)中,其中權(quán)重的定量化取值最終依靠的是專家的理論和經(jīng)驗(yàn)。主要方法有層次分析法(AHP)[3-5]、模糊層次分析法[6]、模糊綜合評(píng)判法[7-8]、灰色聚類法[9]等。但是基于知識(shí)驅(qū)動(dòng)的模型受人為主觀性影響較大,不同情況下必將損失掉空間數(shù)據(jù)本身深層的信息。不同學(xué)者又將數(shù)據(jù)驅(qū)動(dòng)模型運(yùn)用到巖溶塌陷預(yù)測(cè)當(dāng)中,其中通過計(jì)算來確定各影響因素的權(quán)重。數(shù)據(jù)驅(qū)動(dòng)模型主要包括人工神經(jīng)網(wǎng)絡(luò)(ANN)[10]、邏輯回歸[11]、證據(jù)權(quán)法[12]等。這些模型預(yù)測(cè)方法在滑坡、泥石流災(zāi)害區(qū)劃預(yù)測(cè)中已經(jīng)進(jìn)行大量的應(yīng)用。但是這些模型在巖溶塌陷區(qū)劃中只有少量的應(yīng)用,基本呈現(xiàn)出單點(diǎn)開花,個(gè)別學(xué)者應(yīng)用的現(xiàn)象,尤其是數(shù)據(jù)驅(qū)動(dòng)模型迄今為止只有個(gè)別學(xué)者研究使用過,不具備說服力。因此,運(yùn)用證據(jù)權(quán)模型對(duì)巖溶路基發(fā)育強(qiáng)弱進(jìn)行預(yù)測(cè),依托GIS技術(shù)并且有著深厚的統(tǒng)計(jì)學(xué)理論基礎(chǔ),這將為巖溶研究者提供一個(gè)更科學(xué)、直觀的研究背景。
本文選取丹霞互通樞紐巖溶區(qū)域,通過計(jì)算巖溶發(fā)育的先驗(yàn)概率,客觀解釋了各影響因素與巖溶的關(guān)系,選取條件獨(dú)立的證據(jù),利用證據(jù)權(quán)模型結(jié)合GIS技術(shù)確定巖溶發(fā)育強(qiáng)弱的后驗(yàn)概率,從而得出在丹霞互通樞紐建設(shè)影響下巖溶塌陷易發(fā)分區(qū)。最后進(jìn)行預(yù)測(cè)概率的ROC精度評(píng)價(jià)分析。
證據(jù)權(quán)法(the weight of evidence,wofe)是基于貝葉斯準(zhǔn)則的一種離散的多元統(tǒng)計(jì)方法?;谪惾~斯準(zhǔn)則計(jì)算出各證據(jù)層因子的權(quán)重,有效避免了人的主觀性。證據(jù)權(quán)法最開始應(yīng)用于醫(yī)學(xué)領(lǐng)域,在20世紀(jì)80年代BONHAM-CARTER[13]和AGTERBERG等人[14]將證據(jù)權(quán)法應(yīng)用于地質(zhì)礦產(chǎn)領(lǐng)域,現(xiàn)已基本成熟。隨后又有學(xué)者將其應(yīng)用于滑坡、泥石流等地質(zhì)災(zāi)害領(lǐng)域。證據(jù)權(quán)首先進(jìn)行先驗(yàn)概率的計(jì)算,就是用已發(fā)現(xiàn)的巖溶點(diǎn)進(jìn)行計(jì)算;然后進(jìn)行各個(gè)證據(jù)因子權(quán)重的計(jì)算;之后還要進(jìn)行證據(jù)因子的優(yōu)選;最后運(yùn)用GIS將各權(quán)重因子疊加得出巖溶發(fā)育的后驗(yàn)概率圖并進(jìn)行證據(jù)因子獨(dú)立性檢驗(yàn)。
證據(jù)權(quán)模型原理如下。假設(shè)研究區(qū)域柵格數(shù)為N[T],且每個(gè)柵格最多只有一個(gè)巖溶點(diǎn),巖溶點(diǎn)的單元數(shù)為N[K],所以任意單元發(fā)生巖溶的先驗(yàn)概率公式為:
(1)
(2)
當(dāng)Wj+>0,而Wj-<0,表示該證據(jù)有利于巖溶形成;反之,不利于巖溶形成,因此這時(shí)該證據(jù)與巖溶形成無關(guān)。且Wj與證據(jù)層關(guān)系如下所示:
(4)
證據(jù)因子優(yōu)選,在計(jì)算每個(gè)證據(jù)層的權(quán)重后,需要將證據(jù)層二值化處理。但是一個(gè)證據(jù)層可能存在多級(jí)、多個(gè)權(quán)重,就需要確定最優(yōu)切值(cut of fvaule)。最優(yōu)切值可以通過對(duì)比值C(contrast)來確定。C=Wj+-Wj-表示一組空間證據(jù)與一組預(yù)測(cè)目標(biāo)之間的空間關(guān)聯(lián)度量。C值大于0為正相關(guān),且C值越大,關(guān)聯(lián)越大;反之亦然。但是,在一些特定情況下,訓(xùn)練點(diǎn)較少,沒有明確的最大C值,可以采用學(xué)生化檢驗(yàn)[15]。學(xué)生化檢驗(yàn)(the studentizedC)stud(C)是一種顯著性檢驗(yàn)。
(5)
式中:s2(Wj+)和s2(Wj-)表示為Wj+和Wj-的方差。計(jì)算方式如下:
(6)
在顯著性檢驗(yàn)之后,取顯著性水平α=0.05,stud(C)>1.96,認(rèn)為95%與目標(biāo)預(yù)測(cè)相關(guān),表明優(yōu)于其他切值。當(dāng)存在多個(gè)stud(C)大于1.96,可選取N{Bj∩K}最大者[12]。
最后可計(jì)算后驗(yàn)概率。為了計(jì)算簡(jiǎn)便,引入定義的幾率O{K},其計(jì)算如下:
(8)
經(jīng)過推導(dǎo),可以得出:
(10)
本文利用ArcGIS Desktop 10.2軟件等功能和Excel 2013,成功地對(duì)位于廣東省韶關(guān)市仁化縣丹霞樞紐互通區(qū)域進(jìn)行了公路路基巖溶發(fā)育的后驗(yàn)概率計(jì)算,并給出巖溶發(fā)育分區(qū)圖和預(yù)測(cè)精度評(píng)價(jià)。
擬建丹霞樞紐互通位于廣東省韶關(guān)市仁化縣,路線經(jīng)打鐵沖村與新建村?;ネ▍^(qū)設(shè)置有A、B、C、D、E、F、G、H共8個(gè)匝道,總計(jì)長度為6 258.7 m(本文只選取其中交匯的最主要一部分)?;ネ▍^(qū)以橋梁及路基形式通過,線路上部為沖洪積平原,地形平緩,地勢(shì)相對(duì)平坦,以水田、魚塘為主,有一條水溝穿過;丘陵低山,山體高大、山坡較陡,植被發(fā)育,場(chǎng)區(qū)有粵贛高速、水泥路及在建的韶贛鐵路與外界相通,交通條件較好。場(chǎng)區(qū)地貌上屬于低山丘陵單元區(qū)(圖1)。
圖1 研究區(qū)范圍及其概況Fig.1 The overview of the research area
2.2.1證據(jù)因子數(shù)據(jù)準(zhǔn)備
2.2.2地層巖性
基于研究區(qū)的地層地質(zhì)圖和鉆孔數(shù)據(jù),且?guī)r溶發(fā)育在可溶性巖石中,所以地層圖中忽略第四系粉質(zhì)黏土。地層分別是上白堊統(tǒng)南雄組(K2n)含礫砂巖和泥盆系上統(tǒng)錫礦山組(D3x)角礫狀灰?guī)r。表1中是通過計(jì)算確定證據(jù)因子和巖溶的關(guān)系,其中錫礦山組角礫狀灰?guī)r地層的stud(C)>1.96,可將該地層作為本證據(jù)圖層的二值化處理切值。
表1 地層巖性權(quán)重及對(duì)比度
2.2.3地下水
巖溶發(fā)育的一大基本條件就是可溶性水。本文根據(jù)原始數(shù)據(jù)計(jì)算出地下水位距離地面的距離,并做出地下水距地面距離分區(qū)圖。根據(jù)實(shí)際情況,地下水距離地面最小為1.25 m,最大為18.75 m。地下水距地面距離均勻劃分,可將地下水距地面的距離分為5個(gè)區(qū)域,分別是1.25~4.75 m、4.75~8.25 m、8.25~11.75 m、11.75~15.25 m、15.25~18.75 m。由表2可知,區(qū)域1.25~4.75 m可作為本證據(jù)圖層的最優(yōu)切值。且該條件符合理論經(jīng)驗(yàn)認(rèn)知,地下水水位高的地方巖溶發(fā)育強(qiáng)烈。
表2 地下水距地面距離權(quán)重及對(duì)比度
2.2.4頂板厚度
頂板厚度是巖溶發(fā)育頂部距巖層頂面的距離。頂板厚度能夠表明巖溶在可溶性巖層中發(fā)育的位置。頂板厚度對(duì)巖溶發(fā)育有一定的影響,其中頂板越厚,水循環(huán)越差,巖溶可能就不發(fā)育。根據(jù)實(shí)際情況,將頂板厚度分為6個(gè)區(qū)域,分別是:0~3 m、3~6 m、6~9 m、9~12 m、12~18 m、18~35 m。由表3可知,證據(jù)優(yōu)選后,區(qū)域0~3 m可作為本證據(jù)圖層的最優(yōu)切值。并且符合理論經(jīng)驗(yàn)認(rèn)知,在巖層分界面,巖溶發(fā)育強(qiáng)烈。
表3 頂板厚度權(quán)重及對(duì)比度
2.2.5斷層
斷層可以影響地下水流向,對(duì)巖溶發(fā)育有一定的影響。在研究區(qū)有一條斷層破碎帶WF1,物探成果發(fā)現(xiàn)成帶狀分布,向下延伸深度大,近北東走向,視寬度約25 m,影響寬度約60 m,構(gòu)造垂向延伸深度較大,構(gòu)造水平展布長度較大。根據(jù)距斷層距離分區(qū),分別是:0~100 m、100~200 m、200~300 m、300~400 m、<200 m,<300 m。由表4可知,證據(jù)優(yōu)選后,區(qū)域0~100 m可作為本證據(jù)圖層的最優(yōu)切值。
2.2.6證據(jù)因子優(yōu)選
每個(gè)證據(jù)因子幾乎都有多級(jí)分區(qū)。這就需要將每個(gè)證據(jù)圖層進(jìn)行二值化處理。二值化處理之前按照方法確定最優(yōu)切值。經(jīng)計(jì)算(表1~表4),選取的最優(yōu)切值為:地層巖性中灰?guī)r(D3x)、地下水距地面距離1.25~4.75 m、頂板厚度0~3 m、距斷層距離<100 m。
表4 斷層緩沖區(qū)權(quán)重及對(duì)比度
2.3.1條件獨(dú)立性檢驗(yàn)
(11)
式中:m——各證據(jù)圖層存在狀態(tài)的組合個(gè)數(shù);
Pi——第i種組合狀態(tài)下的后驗(yàn)概率;
N[i]——第i種組合狀態(tài)下的柵格數(shù)。
2.3.2危險(xiǎn)區(qū)劃分
通過ArcGIS軟件生成巖溶發(fā)育的后驗(yàn)概率圖,巖溶發(fā)育的地方容易產(chǎn)生巖溶塌陷,尤其是在公路路基施工設(shè)備等增加地面荷載的情況下,所以在施工過程中應(yīng)該提高警惕。在本文中巖溶發(fā)育后驗(yàn)概率圖也可叫做巖溶塌陷風(fēng)險(xiǎn)后驗(yàn)概率圖。后驗(yàn)概率圖在分級(jí)過程中需要確定分級(jí)閾值,通常采用CAPP(Cumulative Area—Posterior Probability)曲線的拐點(diǎn)確定閾值[17-18]。如圖2所示,考慮到先驗(yàn)概率為0.015,則優(yōu)先選擇離先驗(yàn)概率最近的拐點(diǎn),閾值分別為0.007 6、0.014 2、0.022 0。則將后驗(yàn)概率分四級(jí),分別對(duì)應(yīng)穩(wěn)定、較穩(wěn)定、較危險(xiǎn)、危險(xiǎn)(圖3)。
圖2 后驗(yàn)概率面積頻率曲線Fig.2 The cumulative area posterior probability curve
圖3 巖溶塌陷風(fēng)險(xiǎn)后驗(yàn)概率圖Fig.3 The posteriori probability map showing the risk of karst collapse
本文采用ROC(receiver operating characteristic curve)曲線[19-20]進(jìn)行巖溶預(yù)測(cè)準(zhǔn)確性評(píng)價(jià)。ROC預(yù)測(cè)曲線橫坐標(biāo)假陽性為預(yù)測(cè)巖溶區(qū)域(不包括真實(shí)巖溶區(qū)域)百分比,縱坐標(biāo)真陽性為實(shí)際巖溶區(qū)域百分比。ROC曲線下面積(AUC,area under curve)為預(yù)測(cè)概率(圖4),預(yù)測(cè)概率為66.2%。
圖4 ROC預(yù)測(cè)曲線Fig.4 The ROC curve
在公路路基巖溶發(fā)育區(qū),引入常用在礦產(chǎn)資源預(yù)測(cè)中的證據(jù)權(quán)模型,定量的分析計(jì)算每個(gè)證據(jù)因子的權(quán)重,并據(jù)此得出線性公路路基巖溶發(fā)育強(qiáng)度分區(qū)圖,由于人類工程活動(dòng)易引發(fā)巖溶塌陷,公路路基巖溶發(fā)育強(qiáng)度分區(qū)圖也可叫做人類工程作用下公路路基巖溶塌陷發(fā)生危險(xiǎn)分區(qū)圖。危險(xiǎn)分區(qū)圖劃分為四個(gè)區(qū),分別為穩(wěn)定(0.001 8~0.007 6)、較穩(wěn)定(0.007 6~0.014 2)、較危險(xiǎn)(0.014 2~0.022 0)、危險(xiǎn)(0.021 9~0.039 3)。其中穩(wěn)定區(qū)占27.1%、較穩(wěn)定區(qū)占35.4%、較危險(xiǎn)區(qū)占32.9%、危險(xiǎn)區(qū)占4.6%,可根據(jù)后驗(yàn)概率圖,注意危險(xiǎn)區(qū)施工。根據(jù)ROC的精度評(píng)價(jià),預(yù)測(cè)成功率為66.2%(圖4)。巖溶最發(fā)育的地方基本上都在危險(xiǎn)區(qū),并結(jié)合預(yù)測(cè)成功率可知該分區(qū)圖基本符合實(shí)際結(jié)果。
(1)需要指出的是,巖溶發(fā)育強(qiáng)度預(yù)測(cè)和巖溶塌陷預(yù)測(cè)略有不同,一個(gè)是以研究區(qū)內(nèi)已知巖溶發(fā)育點(diǎn)為訓(xùn)練點(diǎn),一個(gè)是以研究區(qū)內(nèi)已知巖溶塌陷為訓(xùn)練點(diǎn)。在本研究區(qū)內(nèi)巖溶塌陷點(diǎn)基本不存在,所以預(yù)測(cè)巖溶發(fā)育強(qiáng)烈區(qū)在人類工程作用下為巖溶塌陷危險(xiǎn)區(qū)。并且本文將巖溶發(fā)育強(qiáng)烈區(qū)作為公路路基施工過程中容易塌陷的區(qū)域,需注意做好防范措施。如果巖溶塌陷的訓(xùn)練點(diǎn)足夠多,可以直接預(yù)測(cè)巖溶塌陷危險(xiǎn)區(qū)。影響巖溶發(fā)育的因素有很多,地表水流域,地下水變化等,但是由于數(shù)據(jù)缺少,所以無法進(jìn)行計(jì)算。并且采用的是巖溶鉆孔數(shù)據(jù),鉆孔數(shù)據(jù)具有很強(qiáng)的目的性,隨機(jī)性不夠強(qiáng)。如有條件,可進(jìn)行巖溶塌陷危險(xiǎn)區(qū)預(yù)測(cè)。
(2)在進(jìn)行證據(jù)因子優(yōu)選時(shí),每個(gè)證據(jù)圖層幾乎都有多級(jí)分區(qū),證據(jù)因子優(yōu)選后,每個(gè)證據(jù)圖層只有兩種狀態(tài),這種情況就忽略了證據(jù)層中其他分區(qū)的狀態(tài),這就造成了最終結(jié)果的系統(tǒng)誤差,這是證據(jù)權(quán)模型的弊端。為了解決這種弊端,成秋明等[21]提出模糊證據(jù)權(quán)模型,這樣就可以最大限度利用證據(jù)圖層中的信息。在計(jì)算證據(jù)條件獨(dú)立性檢驗(yàn)時(shí),運(yùn)用的omnibus test檢驗(yàn)方法是一種近似的檢驗(yàn),大致符合條件獨(dú)立性。