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

(1.大連大學建筑工程學院, 大連 116622； 2.中國地質調查局南京地質調查中心, 南京 210016)

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

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

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

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

1 SINMAP模型及應用

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

(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為土壤導水系數(shù),m2/d；R為地下水側向穩(wěn)態(tài)補給,mm/d；a為單位匯水面積,m2；h為地下水高度,m；q為流量,m3/d。

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

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

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

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

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

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

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

2.1 研究區(qū)域概況

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

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

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

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

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

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

巖土物理參數(shù)來自土體采樣點取樣的室內(nèi)土工試驗，土體抗剪強度指標黏聚力和內(nèi)摩擦角采用快剪實驗得到。風化層厚度z來自于對鉆孔土樣的分析。

工程地質巖組由1∶5 萬區(qū)域地質圖地層巖性歸類獲得，如圖3所示玉壺流域內(nèi)主要分布有六種巖組，Hif(軟硬相見塊狀晶、玻屑熔結凝灰?guī)r夾凝灰質粉砂巖)、Hs(軟硬相間的層狀沉凝灰?guī)r、凝灰質粉砂巖夾流紋晶屑玻屑凝灰?guī)r)、Pb(堅硬塊狀玄武巖等基性噴出巖不易崩滑巖組)、Qg[堅硬塊狀花崗(斑)巖等酸性侵入巖、潛火山巖]、Rr[堅硬塊狀流紋(斑)巖等酸性巖]、Scf(軟硬相間凝灰質粉砂巖、沉凝灰?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模型實現(xiàn)

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

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

為研究DEM精度及降雨條件對地質災害評價結果的影響，區(qū)分了兩種工況，劃分標準見表3。降雨數(shù)據(jù)來自對玉壺流域內(nèi)雨量站數(shù)據(jù)的統(tǒng)計，2016年超強臺風“尼伯特”經(jīng)過溫州帶來了短時間內(nèi)的強降雨，文成縣日最大降雨量達到400 mm，而2017年期間沒有較大臺風過境，日最大降雨量為166.5 mm。圖5為反距離權重法插值得到的研究區(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ā)性制圖柵格累計圖Fig.6 Susceptibility cartographic raster cumulative map

3 結果對比分析

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

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

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

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

同一精度下，工況B(降雨量較大)的總體不穩(wěn)定區(qū)域面積最大，這也符合推測，降雨量會影響T/R范圍進而影響輸出結果，但無條件不穩(wěn)定區(qū)域相同，推測原因可能是由于該區(qū)域汛期降雨量大，日最大降雨量都在100 mm以上，地下水側向補給容易達到飽和狀態(tài)，超過界限后降雨量的增大對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)計Table 4 Statistics of SINMAP model output

4 結論

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

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

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