不同精度下地表穩(wěn)定性模型的評(píng)價(jià)

趙曉東， 王順東， 張?zhí)悾?吳鑫俊， 丁 茜

(1.大連大學(xué)建筑工程學(xué)院, 大連 116622； 2.中國(guó)地質(zhì)調(diào)查局南京地質(zhì)調(diào)查中心, 南京 210016)

地質(zhì)災(zāi)害評(píng)價(jià)，國(guó)外亦稱為滑坡敏感性評(píng)價(jià)，中國(guó)對(duì)地質(zhì)災(zāi)害評(píng)價(jià)的劃分包括易發(fā)性、危險(xiǎn)性、易損性及風(fēng)險(xiǎn)評(píng)價(jià)等，其中易發(fā)性評(píng)價(jià)是地質(zhì)災(zāi)害條件和地質(zhì)災(zāi)害發(fā)生的空間概率統(tǒng)計(jì)分析，注重靜態(tài)地質(zhì)災(zāi)害敏感條件分析和發(fā)生的空間概率預(yù)測(cè)，是地質(zhì)災(zāi)害危險(xiǎn)性和風(fēng)險(xiǎn)評(píng)估的基礎(chǔ)。易發(fā)性評(píng)價(jià)的核心內(nèi)容包括地質(zhì)災(zāi)害特征、空間密度、易發(fā)條件和潛在易發(fā)區(qū)域預(yù)測(cè)評(píng)價(jià)，最終形成研究區(qū)地質(zhì)災(zāi)害易發(fā)性分布圖，其主要目的和用途就是為政府制定土地利用規(guī)劃及災(zāi)害預(yù)警提供基礎(chǔ)資料和依據(jù)[1]。

目前滑坡災(zāi)害評(píng)價(jià)的模型眾多，常用的包括SINMAP、SHALSTAB、CAHSM、TRIGRS、SHETRAN、GEOTOP-FS和SUSHI等[2]。這些評(píng)價(jià)模型大致可分為確定性模型、統(tǒng)計(jì)模型、和非參數(shù)模型[3]。確定性模型是以滑坡發(fā)生的物理力學(xué)為基礎(chǔ)，需要邊坡較為詳細(xì)的巖土力學(xué)參數(shù)，結(jié)果會(huì)受到數(shù)據(jù)量和數(shù)據(jù)精度的影響。相對(duì)于其他模型，確定性模型可以較為客觀地反映邊坡的物理變化，其缺點(diǎn)是受到人力、物力等方面的考慮，缺乏較為詳細(xì)的巖土力學(xué)參數(shù)，該方法通常只能適用于單個(gè)滑坡或者小范圍區(qū)域[4]。在常用模型中，侵蝕地表穩(wěn)定性模型(stability index mapping，SINMAP)由于耦合了無限斜坡模型和穩(wěn)態(tài)水文模型，利用由穩(wěn)定狀態(tài)水文模型獲取的地形濕度指數(shù)、由數(shù)字高程模型(digital elevation model，DEM)獲取的坡度、單位寬度匯水面積等數(shù)據(jù)，結(jié)合地質(zhì)考察資料，可在地理信息系統(tǒng)(GIS)空間分析平臺(tái)上建立定量分析模型，獲得地表穩(wěn)定性分級(jí)，實(shí)現(xiàn)對(duì)研究區(qū)域的地表穩(wěn)定性評(píng)價(jià)[5]。

SINMAP模型是Pack等[6]在Montgomery和Dietrich提出的淺層滑坡穩(wěn)定性模型(shallow landsliding stability，SHALSTAB)基礎(chǔ)上的改進(jìn)模型，考慮了土壤物理參數(shù)的不確定性對(duì)輸出結(jié)果的影響，保留了土壤黏聚力及植物根系對(duì)斜坡穩(wěn)定性的加強(qiáng)作用，可模擬不同降雨條件下斜坡穩(wěn)定性評(píng)價(jià)，非常適用于東南沿海地區(qū)臺(tái)風(fēng)暴雨型地質(zhì)災(zāi)害易發(fā)性的評(píng)價(jià)。

SINMAP模型的輸入?yún)?shù)分為兩類，一類為由DEM轉(zhuǎn)化而來的地形數(shù)據(jù)，另一類是巖土物理力學(xué)參數(shù)。DEM精度會(huì)對(duì)第一類參數(shù)產(chǎn)生直接影響，隨著精度的降低坡度會(huì)更加平緩；DEM分辨率越高，地圖越能清晰呈現(xiàn)真實(shí)的流域面積和坡度。但DEM精度越高，其獲取的成本就越高。要做到技術(shù)參數(shù)可行，經(jīng)濟(jì)成本合理就必須考察不同精度DEM對(duì)易發(fā)性評(píng)價(jià)模型的影響程度。為了探明DEM精度對(duì)應(yīng)用SINMAP模型進(jìn)行地質(zhì)災(zāi)害易發(fā)性評(píng)價(jià)的影響，本文選用免費(fèi)全球30 m DEM、日本ALOS衛(wèi)星的12.5 m和課題組購(gòu)買的5 m三種DEM精度，兩種日最大降雨量工況，對(duì)SINMAP模型的地質(zhì)災(zāi)害易發(fā)性評(píng)價(jià)結(jié)果進(jìn)行了對(duì)比分析研究。

1 SINMAP模型及應(yīng)用

SINMAP是由Pack和Tarboton等開發(fā)出的基于無限斜坡模型和穩(wěn)態(tài)水文模型的侵蝕地表穩(wěn)定性模型。在該模型中，假設(shè)了滑動(dòng)面平行于地表且忽略其邊緣作用，采用地表土層穩(wěn)定的抗滑力與滑動(dòng)力之比為安全系數(shù)，其表達(dá)式為

(1)

(2)

(3)

(4)

(5)

式中：Fs為安全系數(shù)；C為黏聚力因子；w為地形濕度指數(shù)；z為土壤厚度，m；Cr為植物根系產(chǎn)生的黏聚力，N/m2；Cs為土體自身黏聚力，N/m2；θ為地形坡度，(°)；ρs為土壤密度,kg/m3；ρw為水的密度,kg/m3；r為水與土壤密度之比；g為重力加速度,m/s2；φ為內(nèi)摩擦角,(°)；T為土壤導(dǎo)水系數(shù),m2/d；R為地下水側(cè)向穩(wěn)態(tài)補(bǔ)給,mm/d；a為單位匯水面積,m2；h為地下水高度,m；q為流量,m3/d。

穩(wěn)定性指數(shù)(SI)的定義來自安全系數(shù)，模型假設(shè)研究區(qū)域內(nèi)巖土物理參數(shù)呈均勻分布，對(duì)于某一地面點(diǎn)來說，當(dāng)C、tanφ取最小值，R/T取最大值時(shí)，此時(shí)安全系數(shù)Fs最小，該點(diǎn)穩(wěn)定性處于最差狀態(tài)。當(dāng)C、tanφ取最大值，R/T取最小值時(shí)，此時(shí)安全系數(shù)Fs為最大，該點(diǎn)穩(wěn)定性處于最佳狀態(tài)。當(dāng)Fs，min>1時(shí)，該點(diǎn)無條件穩(wěn)定，此時(shí)穩(wěn)定性指數(shù)SI=FSmin；當(dāng)Fs,min<1，F(xiàn)s,max>1時(shí)，該點(diǎn)存在發(fā)生不穩(wěn)定的可能，此時(shí)穩(wěn)定性指數(shù)SI定義為該地面點(diǎn)處于穩(wěn)定狀態(tài)的概率，即SI=Prob(FS>1)；當(dāng)Fs,max<1時(shí)，該點(diǎn)處于不穩(wěn)定狀態(tài)，此時(shí)穩(wěn)定性指數(shù)SI=Prob(Fs>1)=0，根據(jù)SI的大小將邊坡的穩(wěn)定性劃分為5個(gè)等級(jí)，SINMAP輸出結(jié)果分類見表1，關(guān)于SINMAP的詳細(xì)理論參見文獻(xiàn)[6]。

表1 SINMAP模型輸出值分類Table 1 Output value classification of SINMAP model

中外學(xué)者對(duì)滑坡穩(wěn)定性的預(yù)測(cè)大多應(yīng)用基于ArcGIS軟件的統(tǒng)計(jì)模型，復(fù)合了確定性模型、信息量、邏輯回歸、潛勢(shì)度等多種模型及方法[7-10]。將SINMAP模型應(yīng)用于地質(zhì)災(zāi)害易發(fā)性評(píng)價(jià)中國(guó)也早有先例，朱遠(yuǎn)樂等[11]將依據(jù)層次分析法得到的滑坡預(yù)測(cè)易發(fā)性評(píng)價(jià)柵格圖與應(yīng)用SINMAP模型集成分析得到尾礦庫(kù)地表穩(wěn)定性指數(shù)分布圖。姚志永[12]對(duì)比了SINMAP和TRIGRS(transient rainfall infiltration and gridbased regional slope-stability model)兩種模型在湖南省辰溪縣暴雨型淺層滑坡易發(fā)性評(píng)價(jià)中的差異，得出在該區(qū)域SINMAP模型更為準(zhǔn)確。周偉[13]研究了在白龍江流域，邏輯回歸模型在整體上優(yōu)于SINMAP模型但針對(duì)單個(gè)滑坡而言后者能更好反映坡體穩(wěn)定性。莊建琦等[14]研究發(fā)現(xiàn)在降雨強(qiáng)度為30 mm的條件下SINMAP模型與其建立的黃土地區(qū)淺層滑坡發(fā)育危險(xiǎn)性評(píng)價(jià)模型結(jié)果最吻合。

中外學(xué)者將穩(wěn)定性模型引入地質(zhì)災(zāi)害評(píng)價(jià)非常廣泛，Rabonza等[15]在菲律賓中部使用SINMAP軟件，將地形、土壤強(qiáng)度、水文參數(shù)和DTM數(shù)據(jù)輸入模型，計(jì)算相應(yīng)的安全系數(shù)，利用SINMAP生成的滑坡圖與2002—2014年高分辨率衛(wèi)星圖像得到的滑坡清查結(jié)果高度一致。Michel等[16]對(duì)巴西南部山區(qū)的達(dá)河流域應(yīng)用SHALSTAB和SINMAP分別進(jìn)行滑坡易發(fā)性制圖，得出該流域SINMAP的模擬結(jié)果穩(wěn)定區(qū)域預(yù)測(cè)面積過大，而SHALSTAB模型模擬結(jié)果良好。

前文敘述過SINMAP模型是在SHALSTAB模型基礎(chǔ)上改進(jìn)而成，SHALSTAB模型的提出者Dietrich等[17]及研究者Ramos等曾指出：模型的模擬精度會(huì)隨著輸入?yún)?shù)精度的提高而提高，地形數(shù)據(jù)分辨率以及土壤物理參數(shù)對(duì)模型模擬效果的影響尤為顯著，同樣推測(cè)DEM精度也會(huì)對(duì)SINMAP模型的輸出結(jié)果產(chǎn)生很大影響。

如何選取量化指標(biāo)進(jìn)行模擬結(jié)果的對(duì)比，不同學(xué)者使用了不同的方法，大多數(shù)研究人員應(yīng)用野外獲得的土壤物理參數(shù)如內(nèi)摩擦角、抗剪強(qiáng)度、飽和密度，但這些參數(shù)在一定區(qū)域內(nèi)變化較大，而且?guī)в姓`差，Seefelder等[18]提出了一種用于山區(qū)滑坡問題分區(qū)的方法，該方法有效且最適合于水文盆地，他認(rèn)為可以用流域內(nèi)工程地質(zhì)巖組對(duì)區(qū)域進(jìn)行劃分工況。Regiane Mara Sbroglia應(yīng)用這種劃分方法對(duì)巴西某山區(qū)地帶進(jìn)行了研究，模擬結(jié)果良好[19]。

2 研究區(qū)域數(shù)據(jù)整理

2.1 研究區(qū)域概況

研究數(shù)據(jù)范圍處于溫州市飛云江流域，以玉壺流域重點(diǎn)研究，研究區(qū)域位置見圖1。溫州陸域面積11 784 km2，其中山地面積占78.2%，地勢(shì)從西南向東北呈現(xiàn)梯形傾斜，有洞宮、括蒼、雁蕩山脈，泰順的白云尖，海拔1 611 m，為全市最高峰。屬亞熱帶季風(fēng)氣候區(qū)，年降水量1 500～1 900 mm，8—10月是臺(tái)汛期，多年平均臺(tái)風(fēng)登陸和受臺(tái)風(fēng)影響3.5次/年。玉壺流域位于溫州市東南地區(qū)，隸屬于文成縣，面積103.9 km2。

圖1 研究區(qū)域位置Fig.1 Location of study area

受地質(zhì)構(gòu)造和第四紀(jì)晚期快速海退的控制，山區(qū)河流深切、相對(duì)高差大、坡面陡峭，坡度多大于30°，坡面穩(wěn)定性差；組成山體的巖石以侏羅系白堊系火山巖為主，山坡表面和坡腳殘積層一般1～5 m；山區(qū)植被覆蓋率高，河谷和緩坡多被開墾成梯田或?yàn)榫用窠ㄔO(shè)用地，山區(qū)人口分散。域內(nèi)山地多，山體穩(wěn)定性差，人為工程活動(dòng)強(qiáng)，臺(tái)風(fēng)暴雨頻繁，是地質(zhì)災(zāi)害高易發(fā)區(qū)和多發(fā)區(qū)，是浙江省突發(fā)性地質(zhì)災(zāi)害最嚴(yán)重的市。據(jù)不完全統(tǒng)計(jì)，已查明的地質(zhì)災(zāi)害隱患點(diǎn)近2 000處，約占全省的20%以上；受威脅人口30 007人，占全省的21.5%。

2.2 地形數(shù)據(jù)

5 m數(shù)字高程數(shù)據(jù)為課題組購(gòu)入數(shù)據(jù)，使用ArcGIS10.2軟件從DEM當(dāng)中提取山脊線，繪制流域分布圖，依據(jù)中華人民共和國(guó)水利行業(yè)標(biāo)準(zhǔn)控制流域面積并命名[20]，最終選取飛云江上游玉壺流域作為研究區(qū)域。12.5 m DEM數(shù)據(jù)來源于日本ALOS衛(wèi)星，30 m DEM取自NASA全球30 m高程數(shù)據(jù)，因其相鄰兩種精度比值大小相近，便于比較，故均取其原始大小計(jì)算。利用ArcGIS軟件Hillshade工具從DEM數(shù)據(jù)當(dāng)中提取出玉壺流域山體陰影(圖2)。

2.3 巖土物理力學(xué)參數(shù)

巖土物理參數(shù)來自土體采樣點(diǎn)取樣的室內(nèi)土工試驗(yàn)，土體抗剪強(qiáng)度指標(biāo)黏聚力和內(nèi)摩擦角采用快剪實(shí)驗(yàn)得到。風(fēng)化層厚度z來自于對(duì)鉆孔土樣的分析。

工程地質(zhì)巖組由1∶5 萬(wàn)區(qū)域地質(zhì)圖地層巖性歸類獲得，如圖3所示玉壺流域內(nèi)主要分布有六種巖組，Hif(軟硬相見塊狀晶、玻屑熔結(jié)凝灰?guī)r夾凝灰質(zhì)粉砂巖)、Hs(軟硬相間的層狀沉凝灰?guī)r、凝灰質(zhì)粉砂巖夾流紋晶屑玻屑凝灰?guī)r)、Pb(堅(jiān)硬塊狀玄武巖等基性噴出巖不易崩滑巖組)、Qg[堅(jiān)硬塊狀花崗(斑)巖等酸性侵入巖、潛火山巖]、Rr[堅(jiān)硬塊狀流紋(斑)巖等酸性巖]、Scf(軟硬相間凝灰質(zhì)粉砂巖、沉凝灰?guī)r、粉砂巖、砂礫巖、砂巖)，研究區(qū)域內(nèi)巖組分布見圖3，具體占比見表2。

圖2 玉壺流域山體陰影圖Fig.2 Hillshade map of yuhu basin

圖3 巖組分布Fig.3 Distribution of rock group

表2 巖土參數(shù)分布

2.4 SINMAP模型實(shí)現(xiàn)

SINMAP模型的實(shí)現(xiàn)環(huán)境為ArcGIS10.2版本，軟件采用柵格單元計(jì)算，將采樣點(diǎn)巖土參數(shù)數(shù)據(jù)整理成覆蓋整個(gè)玉壺流域的柵格數(shù)據(jù)(黏聚力c、內(nèi)摩擦角φ、風(fēng)化層厚度z、密度ρ)。C、R/T、tanφ為可變參數(shù)，根據(jù)研究區(qū)域內(nèi)采樣點(diǎn)巖土力學(xué)參數(shù)分別指定C、tanφ的上限與下限，R/T參數(shù)由室內(nèi)土工試驗(yàn)及參考其他學(xué)者研究確定取值范圍，并假定它們?cè)谥付ǚ秶鷥?nèi)均勻分布[12-13]。

地形坡度由DEM在ArcGIS中使用坡度計(jì)算獲得。單位匯水面積的取值與每個(gè)柵格內(nèi)模擬水流流向有關(guān)，流向的確定采用水文分析當(dāng)中常用的D8算法，流向與柵格邊平行取邊長(zhǎng)，與對(duì)角線平行則取柵格對(duì)角線長(zhǎng)度。匯水面積A在數(shù)值上等于流量與柵格面積的乘積，其中流量由DEM在ArcGIS中使用自行開發(fā)的計(jì)算工具獲取。可以理解為單位匯水面積越大，徑流強(qiáng)度越大，對(duì)坡體沖刷能力越強(qiáng)，圖4為單位匯水面積示意圖，SHALSTAB模型提出者Dietrich繪制，有修改[21]。

為研究DEM精度及降雨條件對(duì)地質(zhì)災(zāi)害評(píng)價(jià)結(jié)果的影響，區(qū)分了兩種工況，劃分標(biāo)準(zhǔn)見表3。降雨數(shù)據(jù)來自對(duì)玉壺流域內(nèi)雨量站數(shù)據(jù)的統(tǒng)計(jì)，2016年超強(qiáng)臺(tái)風(fēng)“尼伯特”經(jīng)過溫州帶來了短時(shí)間內(nèi)的強(qiáng)降雨，文成縣日最大降雨量達(dá)到400 mm，而2017年期間沒有較大臺(tái)風(fēng)過境，日最大降雨量為166.5 mm。圖5為反距離權(quán)重法插值得到的研究區(qū)域內(nèi)降雨量分布圖。

圖4 單位匯水面積示意圖Fig.4 Schematic diagram of unit catchment area

表3 工況劃分Table 3 Working condition of division

圖5 降雨量分布圖Fig.5 Rainfall distribution map

圖6 易發(fā)性制圖柵格累計(jì)圖Fig.6 Susceptibility cartographic raster cumulative map

3 結(jié)果對(duì)比分析

圖6為各種工況下的易發(fā)性制圖，按照模型提出者的分級(jí)標(biāo)準(zhǔn)進(jìn)行分類，通過統(tǒng)計(jì)易發(fā)性制圖(圖7)的數(shù)據(jù)得到SINMAP模型輸出值統(tǒng)計(jì)(表4)和柵格累計(jì)(圖6)，表4中列出了不同DEM精度與降雨量組合后模型的輸出值分類結(jié)果。

圖7 易發(fā)性制圖Fig.7 Susceptibility map

可以看出三種精度結(jié)果差異較大，主要是穩(wěn)定和無條件不穩(wěn)定區(qū)域的差別，降雨量較小的A工況下，5 m精度的無條件不穩(wěn)定區(qū)域?yàn)?2.1%，穩(wěn)定區(qū)域?yàn)?1.6%；12.5 m精度的無條件不穩(wěn)定區(qū)域?yàn)?.5%，穩(wěn)定區(qū)域?yàn)?4.1%；30 m精度的無條件不穩(wěn)定區(qū)域?yàn)?.5%，穩(wěn)定區(qū)域?yàn)?4.2%；DEM精度較高的5 m條件下識(shí)別出的不穩(wěn)定區(qū)域接近12.5 m精度的兩倍，30 m精度的3倍，足見DEM精度對(duì)模型輸出值的影響非常大，主要是影響無條件不穩(wěn)定區(qū)域的大小。相同降雨量不同DEM精度間的對(duì)比很明顯表現(xiàn)出一條規(guī)律：精度越高，模型識(shí)別出的無條件不穩(wěn)定區(qū)域越多。

對(duì)比工況A下三種DEM精度輸出結(jié)果：5 m精度下無條件不穩(wěn)定區(qū)域?yàn)?2.1%，12.5 m精度的無條件不穩(wěn)定區(qū)域?yàn)?.5%，30 m精度的無條件不穩(wěn)定區(qū)域?yàn)?.5%，像元大小同樣增大2.5倍，由5 m增大到12.5 m時(shí)無條件不穩(wěn)定區(qū)域減小量為5.6%，而12.5 m增大到30 m時(shí)無條件不穩(wěn)定區(qū)域面積減小量?jī)H為2%，這說明精度降低后模型識(shí)別不穩(wěn)定區(qū)域能力的影響迅速減弱。

同一精度下，工況B(降雨量較大)的總體不穩(wěn)定區(qū)域面積最大，這也符合推測(cè)，降雨量會(huì)影響T/R范圍進(jìn)而影響輸出結(jié)果，但無條件不穩(wěn)定區(qū)域相同，推測(cè)原因可能是由于該區(qū)域汛期降雨量大，日最大降雨量都在100 mm以上，地下水側(cè)向補(bǔ)給容易達(dá)到飽和狀態(tài)，超過界限后降雨量的增大對(duì)SINMAP模型不穩(wěn)定區(qū)域的增加影響較小。相差三倍的日最大降雨量數(shù)據(jù)在5、12.5、30 m三種精度下模型穩(wěn)定區(qū)域面積分別相差2%、3.6%、4.5%，無條件不穩(wěn)定區(qū)域面積相同。

表4 SINMAP模型輸出統(tǒng)計(jì)Table 4 Statistics of SINMAP model output

4 結(jié)論

(1)DEM精度對(duì)SINMAP模型下地質(zhì)災(zāi)害易發(fā)性評(píng)價(jià)結(jié)果影響非常大，特別是對(duì)極端情況的影響。隨著DEM精度的提高，模型輸出的無條件不穩(wěn)定區(qū)域擴(kuò)大，無條件穩(wěn)定區(qū)域減小，中間過渡地帶變化并不明顯。在降雨量相同的情況下，DEM精度較高的5 m條件下識(shí)別出的不穩(wěn)定區(qū)域接近12.5 m精度的兩倍，30 m精度的3倍。

(2)DEM精度對(duì)SINMAP模型下地質(zhì)災(zāi)害易發(fā)性評(píng)價(jià)結(jié)果影響是非線性的，高精度轉(zhuǎn)為中精度，中精度轉(zhuǎn)為低精度，模型輸出無條件不穩(wěn)定區(qū)域下降速率降低很快。當(dāng)DEM精度降低到一定范圍后，精度的降低對(duì)模型輸出的地質(zhì)災(zāi)害易發(fā)性制圖影響很小。

(3)日最大降雨量超過100 mm情況下降雨量對(duì)模型輸出值的影響不大，相差3倍的日最大降雨量數(shù)據(jù)在5、12.5、30 m三種精度下模型穩(wěn)定區(qū)域面積僅相差2%、3.6%、4.5%，無條件不穩(wěn)定區(qū)域面積相同。