SHALSTAB模型在浙南林溪流域滑坡預(yù)測(cè)中的應(yīng)用

孫 強(qiáng),張?zhí)?伍劍波,王赫生,朱延輝,韓 帥

(中國地質(zhì)調(diào)查局南京地質(zhì)調(diào)查中心,江蘇 南京 210016)

浙南林溪流域所屬的浙閩丘陵山地區(qū)是我國群發(fā)性滑坡重點(diǎn)防治區(qū)之一,多年飽受臺(tái)風(fēng)暴雨襲擾。該區(qū)域地形起伏大,構(gòu)造發(fā)育,巖體破碎、風(fēng)化強(qiáng)烈,頻發(fā)的極端降雨與人類工程活動(dòng)誘發(fā)了大量的地質(zhì)災(zāi)害。盡管滑坡形成和觸發(fā)的因素是眾所周知的,但其發(fā)生與誘發(fā)因素之間的映射關(guān)系仍然是一個(gè)巨大的挑戰(zhàn),也一直是地質(zhì)災(zāi)害領(lǐng)域研究的難點(diǎn)和熱點(diǎn)[1]?；乱装l(fā)性評(píng)價(jià)為地質(zhì)災(zāi)害管理提供重要信息,是提升防災(zāi)減災(zāi)效果的重要措施之一[2]。

國內(nèi)外對(duì)滑坡預(yù)測(cè)進(jìn)行了大量的研究,理論方法主要有統(tǒng)計(jì)分析模型、確定性模型等。統(tǒng)計(jì)分析模型一般較少考慮滑坡發(fā)生的內(nèi)因,注重對(duì)影響滑坡的各種因子與滑坡分布之間函數(shù)關(guān)系進(jìn)行分析,由于統(tǒng)計(jì)分析模型是基于數(shù)據(jù)進(jìn)行預(yù)測(cè), 因此在一個(gè)地區(qū)所建立的模型一般不適用于其他地區(qū)[3]。確定性模型以物理力學(xué)為基礎(chǔ),其各項(xiàng)參數(shù)及計(jì)算過程有實(shí)際的物理意義。在進(jìn)行區(qū)域斜坡穩(wěn)定性的計(jì)算,由于土體的不均勻性,輸入?yún)?shù)通常是以點(diǎn)代面,模型的準(zhǔn)確性受限于輸入?yún)?shù)的數(shù)據(jù)量和數(shù)據(jù)精度,輸入的參數(shù)越多,精度越高,模型與實(shí)際情況的吻合度越高,但實(shí)際計(jì)算過程中,對(duì)每個(gè)計(jì)算的單元進(jìn)行采樣測(cè)試是難以實(shí)現(xiàn)的,也是不經(jīng)濟(jì)的。由于缺乏詳細(xì)的土體力學(xué)參數(shù),該方法通常僅適用于單個(gè)滑坡或者小范圍區(qū)域。但近幾年,確定性模型已經(jīng)由單一的邊坡穩(wěn)定性模型,逐漸發(fā)展為能進(jìn)行區(qū)域斜坡穩(wěn)定性分析的耦合模型。

目前的確定性模型有很多種,如CHASM、SHALSTAB、SINMAP、TRIGRS、SHETRAN、GEOTOP FS和SUSHI等。選擇合適的模型對(duì)正確認(rèn)識(shí)和預(yù)測(cè)滑坡有著重要的意義。BATHURST J C等[4]強(qiáng)調(diào)用水文空間分布的物理模型來預(yù)測(cè)滑坡。將水文模型與滑坡穩(wěn)定性進(jìn)行耦合是目前被廣泛認(rèn)可的降雨型滑坡預(yù)測(cè)方法[5]。在這些模型中,淺層滑坡穩(wěn)定性模型(SHALSTAB)耦合了穩(wěn)態(tài)水文模型與無限邊坡穩(wěn)定模型,模型中的假設(shè)條件和降雨因素的考量較符合浙南林溪流域的滑坡破壞模式。

國內(nèi)許多學(xué)者使用SHALSTAB模型進(jìn)行區(qū)域的滑坡預(yù)測(cè)。同霄等[6]對(duì)淺層黃土滑坡的預(yù)測(cè)研究顯示淺層滑坡分布于雨水匯積的區(qū)域,與實(shí)際情況較吻合?？党萚7]利用SHALSTAB模型進(jìn)行黃土滑坡預(yù)測(cè),模型使用了近70個(gè)勘察點(diǎn)的土體力學(xué)參數(shù),計(jì)算結(jié)果顯示正確率為70. 23%。王佳佳等[8]依據(jù)三峽庫區(qū)二期和三期滑坡勘察資料,統(tǒng)計(jì)了3類滑坡滑體的重度和抗剪強(qiáng)度參數(shù),基于SHALSTAB模型對(duì)巴東新城進(jìn)行了滑坡災(zāi)害預(yù)測(cè)。陳悅麗等[9]通過Monte Carlo法隨機(jī)生成1 000組土體黏聚力和內(nèi)摩擦角的值,計(jì)算斜坡失穩(wěn)的概率,模擬結(jié)果顯示當(dāng)滑坡危險(xiǎn)性閾值取10%時(shí),有90%以上的歷史滑坡點(diǎn)落在危險(xiǎn)區(qū)域內(nèi)。由于模型中需要的土體力學(xué)參數(shù)通常需要實(shí)驗(yàn)分析獲取,如果研究區(qū)域面積較大,考慮人力、物力和時(shí)間等的成本,并且試驗(yàn)測(cè)試本身存在誤差,土體力學(xué)參數(shù)具有高的空間變異性[10]。目前的研究證實(shí),對(duì)模型輸入的土體力學(xué)參數(shù)進(jìn)行合理的取值是非常重要的[11]。

土體的礦物組成對(duì)抗剪強(qiáng)度和水理性質(zhì)有著重要的控制作用[12-13],而土體的礦物組成與其母巖巖性是相關(guān)的。對(duì)各種地層巖性分布區(qū)進(jìn)行分析測(cè)試進(jìn)而獲得該區(qū)域土體力學(xué)參數(shù),使模型更符合實(shí)際情況,同時(shí)滿足經(jīng)濟(jì)和可行性的需求。在考慮到土體力學(xué)參數(shù)分布不均,基于同種巖性的土體力學(xué)參數(shù)的概率分布特征進(jìn)行判斷斜坡的穩(wěn)定性更加合理有效。與王佳佳等[8]、陳悅麗等[9]方法不同,本文模型采用的巖土力學(xué)參數(shù)為大樣本實(shí)測(cè)數(shù)據(jù),在充分分析不同巖性的土體力學(xué)參數(shù)概率分布特征的基礎(chǔ)上,通過SHALSTAB模型進(jìn)行了林溪流域滑坡預(yù)測(cè)。

1 研究區(qū)概況

研究區(qū)位于浙江省溫州市中部的林溪流域,橫跨瑞安市與甌海區(qū)(圖1(a))。林溪流域面積約61 km2,是飛云江下游的一條支流,干流長度8 km左右,流域內(nèi)溪流縱橫,在瑞安市湖嶺鎮(zhèn)南部上街村匯入飛云江。

流域最高海拔位于北側(cè)瑞安市境內(nèi)大巖頭，海拔為905.4 m,最低海拔位于流域南部的上街村,海拔12.1 m,東側(cè)的五峰尖為瑞安市和甌海區(qū)的界山(圖1(b))。流域內(nèi)山高坡陡、溝窄谷深,上游的支溝水系為源發(fā)性溪流,溝谷比降一般100‰～300‰,水流湍急,每年汛期河水位暴漲暴落,常引起溪流塌岸。在林源村一帶,林溪被截?cái)?修筑有林溪水庫,受水庫調(diào)蓄作用,下游溪流潺湲,谷地開闊。據(jù)多年統(tǒng)計(jì)數(shù)據(jù),流域全年平均降雨量為1 200 mm左右,降雨集中,7～9月降水量占近全年降水量的70%以上。

圖1 研究區(qū)地理位置(a)和地勢(shì)圖(b)Fig. 1 Geographical location(a) and topography(b)of the study area

林溪流域?qū)偃A南褶皺系,浙東南褶皺帶泰順—祖州斷拗的中部。燕山期地質(zhì)構(gòu)造-火山-巖漿活動(dòng)強(qiáng)烈,斷裂十分發(fā)育,大面積分布火山巖。流域內(nèi)分布的主要地層巖性有燕山期花崗巖，上侏羅統(tǒng)高塢組(J3g)凝灰?guī)r,下白堊統(tǒng)館頭組(K1g)凝灰?guī)r、粉砂巖等。受構(gòu)造影響,區(qū)域巖體風(fēng)化強(qiáng)烈,松散層厚度大,為滑坡發(fā)育提供了基礎(chǔ)地質(zhì)條件。隨著城鎮(zhèn)化的發(fā)展,區(qū)內(nèi)居民建房、公路建設(shè)等大量工程活動(dòng)產(chǎn)生邊坡開挖,并在近年來極端降雨頻發(fā)的氣象條件下,誘發(fā)了大量滑坡。2016年開展的1∶1萬地質(zhì)災(zāi)害調(diào)查成果顯示該區(qū)自1997年以來,共發(fā)育滑坡34處,密度達(dá)0.56處/km2。

2 SHALSTAB計(jì)算原理

MONTGOMERY D R和DIETRICH W E[14]于1994年提出了SHALSTAB模型,此模型基于莫爾-庫侖破壞準(zhǔn)則,結(jié)合了無限邊坡模型和穩(wěn)態(tài)水文模型。

2.1 無限邊坡模型

由于研究區(qū)內(nèi)的滑坡規(guī)模小,滑面埋深淺,其滑體長度遠(yuǎn)大于滑體厚度,無限邊坡穩(wěn)定性模型是非常有效的。該模型主要考慮巖土體的下滑力和抗滑力,忽略了邊緣效應(yīng)。無限邊坡穩(wěn)定性的模型是基于莫爾-庫侖破壞準(zhǔn)則,其修正公式為

ρs·g·z·sinθ·cosθ=c+(ρs·g·z·cos2θ-

ρw·g·h·cos2θ)·tanφ,

(1)

式中:c是土體的黏聚力,Pa;φ是土的內(nèi)摩擦角,°;ρs是濕土的密度,kg·m-3;z是土體厚度,m;θ是斜坡坡度,°;ρw是水的密度,kg·m-3;h是地下水位,m;g是重力加速度,m·s-2。

2.2 穩(wěn)態(tài)水文模型

滑坡的穩(wěn)定條件與地下水位直接相關(guān)。因此,建立水文模型來估算地下水位是十分必要的。TOPMODEL模型是最常見描述邊坡穩(wěn)態(tài)潛流的水文地質(zhì)概念模型[15],其表達(dá)式是達(dá)西定律方程的變形方程式

q·a=b·ks·h·sinθ·cos,

(2)

式中:q是有效降雨強(qiáng)度,m·d-1;a是積水面積,m2;b是排水區(qū)寬度,m;ks是土體飽和滲透系數(shù)(假設(shè)沿土層深度是均勻的),m·d-1;h是地下水位,m;g是重力加速度,m·s-2。

根據(jù)公式(1)和公式(2)可以建立無條件穩(wěn)定、無條件不穩(wěn)定及基于q/T函數(shù)關(guān)系的SHALSTAB方程組。

無條件不穩(wěn)定,即土體無地下水時(shí),斜坡仍然無法保持穩(wěn)定。

(3)

另外一種極端條件發(fā)生在h/z為1(完全飽和)時(shí),土體抗剪強(qiáng)度較大,斜坡仍然保持穩(wěn)定。

(4)

SHALSTAB最終方程為

(5)

SHALSTAB模型需要輸入的參數(shù)為c,φ,ρs和z,其他變量(a,b和θ)從數(shù)值高程模型(DEM)提取。SHALSTAB主要是根據(jù)有效降雨強(qiáng)度q和土壤的導(dǎo)水系數(shù)T之間的比值(q/T)與土壤強(qiáng)度在地形上所能產(chǎn)生的最大飽和狀態(tài)進(jìn)行斜坡穩(wěn)定性的判斷。將斜坡分為7個(gè)穩(wěn)定等級(jí),除了無條件穩(wěn)定和無條件不穩(wěn)定2個(gè)等級(jí)外,其他5個(gè)等級(jí)是q/T的一個(gè)函數(shù)(按照倍數(shù)關(guān)系)，表1列出了其他5個(gè)等級(jí)的q/T斷點(diǎn)值。

表1 q/T及l(fā)og(q/T)換算表

3 輸入?yún)?shù)

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

模型所需的坡度、積水面積等參數(shù)均由數(shù)值高程模型(DEM)獲取。DEM分辨率越高,地圖越能清晰呈現(xiàn)真實(shí)的流域面積和坡度,但是沒有一個(gè)“完美”的DEM分辨率可以適用所有地區(qū)的邊坡穩(wěn)定性分析[16]。WU S M等[17]發(fā)現(xiàn)當(dāng)其他參數(shù)精度有限時(shí),單純的提高DEM分辨率對(duì)評(píng)價(jià)結(jié)果影響很小。考慮其他輸入?yún)?shù)的取樣密度,模型中使用的DEM來源于浙江省地理信息與測(cè)繪局(浙江溫州地區(qū)2013年測(cè)繪成果),精度為5 m。區(qū)域內(nèi)的滑坡幾何特征均為實(shí)測(cè)。通過劃定滑坡邊界形成滑坡分布圖,并用于檢驗(yàn)?zāi)Ｐ托Ч?。在本研究?所有的滑坡均為發(fā)生過或者有明顯變形跡象的滑坡。

3.2 土體厚度

土體厚度的獲取方法有很多種,一般通過動(dòng)態(tài)圓錐貫入儀測(cè)定土體的厚度,也可以通過滑坡體的觀測(cè)估計(jì)土體深度[18]。本文土體的厚度數(shù)據(jù)源于區(qū)域地質(zhì)災(zāi)害調(diào)查,調(diào)查精度為1∶1萬。調(diào)查嚴(yán)格按照2014年中國地質(zhì)調(diào)查局頒布的《DD2014—08地質(zhì)災(zāi)害調(diào)查評(píng)價(jià)技術(shù)要求(1∶50 000)》[19],調(diào)查路線間距一般為300～500 m,調(diào)查點(diǎn)密度3.4處/km2。研究區(qū)內(nèi)人類工程活動(dòng)強(qiáng)烈,有大量的人工邊坡揭露了土體的厚度。通過統(tǒng)計(jì)有土體厚度調(diào)查數(shù)據(jù)的118個(gè)調(diào)查點(diǎn),進(jìn)行回歸克里格插值形成土體厚度的柵格面。結(jié)果顯示區(qū)域內(nèi)南向坡土體較厚,土體深度最大8.60 m,平均深度為4.01 m左右,隨著坡度的升高,土體厚度逐漸變薄(圖2)。

圖2 土體厚度與坡度(a)、坡向(b)關(guān)系圖Fig. 2 Relation between soil thickness and slope gradient(a), slope aspect(b)

3.3 土體參數(shù)

研究區(qū)主要出露3種巖性,分別為粉砂巖、凝灰?guī)r以及花崗巖。沿調(diào)查路線均勻布點(diǎn)共采集了208組樣品(粉砂巖30組、凝灰?guī)r52組以及花崗巖126組),進(jìn)行土工試驗(yàn)后,得到了土體的抗剪強(qiáng)度參數(shù)、重度、飽和滲透系數(shù)和粒徑分布。由于研究區(qū)滑坡的厚度小,發(fā)生破壞速度快,土體的抗剪強(qiáng)度參數(shù)(φ和c)通過原狀土樣的直剪試驗(yàn)測(cè)定。3種巖性的抗剪強(qiáng)度參數(shù)分布如圖3所示:

圖3 不同巖性區(qū)域土體的抗剪強(qiáng)度參數(shù)特征Fig. 3 Parameter characteristics of shear strength of different lithological soil masses

粉砂巖地區(qū)的土體c值較高,但φ值較小;花崗巖地區(qū)的土體c值較低,φ值較大,且樣品間的差異明顯(表2)。土體的c值變異系數(shù)54.0%～69.2%,而φ值則變化較小,變異系數(shù)13.4%～21.9%。

表2 不同巖性區(qū)域土體的c、φ值特征

多個(gè)研究[20-21]表明土體的黏聚力和內(nèi)摩擦角呈正態(tài)分布。研究區(qū)的土體抗剪強(qiáng)度參數(shù)的分布也符合正態(tài)分布,根據(jù)實(shí)測(cè)值擬合的正態(tài)分布曲線參數(shù)如表3所示。

表3 不同巖性區(qū)域土體抗剪強(qiáng)度的正態(tài)分布參數(shù)

模型采用單側(cè)置信區(qū)間下限值(置信水平為0.90)進(jìn)行預(yù)測(cè)模擬。對(duì)研究區(qū)土體重度及滲透系數(shù)的考慮與土體厚度相同,通過回歸克里格插值形成連續(xù)的柵格面。

4 討論

4.1 穩(wěn)定性結(jié)果及評(píng)價(jià)分析

將公式(3)、公式(4)、公式(5)計(jì)算結(jié)果進(jìn)行疊加，計(jì)算各單元的穩(wěn)定性指數(shù),按照表1進(jìn)行重分類后得到穩(wěn)定指數(shù)分級(jí)圖(圖4)。

圖4 穩(wěn)定性指數(shù)分級(jí)圖Fig. 4 Stability index map

模型結(jié)果顯示,研究區(qū)不穩(wěn)定的區(qū)域較連續(xù),主要位于斜坡的下部及深切河谷的兩側(cè)且坡度陡峭的地區(qū)。無條件穩(wěn)定的區(qū)域占研究區(qū)總面積的68.10%,主要為山體頂部,土體厚度薄及地形平坦的區(qū)域。

4.2 模型質(zhì)量評(píng)價(jià)

模型質(zhì)量通過評(píng)估實(shí)際發(fā)生的滑坡與預(yù)測(cè)滑坡之間的空間重合來實(shí)現(xiàn)的,SORBINO G等[22]提出捕獲率(SI)和誤判率(EI)兩項(xiàng)量化指數(shù)能有效檢驗(yàn)?zāi)Ｐ唾|(zhì)量(圖5)。SI指的是模型預(yù)測(cè)滑坡區(qū)域與實(shí)際發(fā)生滑坡重合區(qū)域占實(shí)際發(fā)生滑坡區(qū)域的百分比,而EI是模型預(yù)測(cè)滑坡區(qū)域占實(shí)際未發(fā)生滑坡區(qū)域的百分比。

(12)

(13)

SORBINO G等[22]發(fā)現(xiàn)模型中輸入不同的參數(shù),預(yù)測(cè)結(jié)果的SI值與EI值均發(fā)生變化,但兩者之間呈線型相關(guān), SI/EI比值具有更好的評(píng)價(jià)模型效果,比值越高說明模型的預(yù)測(cè)效果越好。用2016年開展的林溪流域1∶1萬地質(zhì)災(zāi)害調(diào)查獲得的實(shí)際滑坡數(shù)據(jù)作為驗(yàn)證,不同log(q/T)水平為不穩(wěn)定區(qū)判別標(biāo)準(zhǔn)的SI值與EI值如表4所示。

表4 不同log(q/T)水平的量化指數(shù)

隨著log(q/T)水平的提高,模型預(yù)測(cè)的滑坡面積逐漸擴(kuò)大,SI值逐漸增大, EI值同時(shí)增大, SI/EI值由4.9減小到2.2。單純提高log(q/T)水平,滑坡與預(yù)測(cè)滑坡之間的空間重合度增大,同時(shí)導(dǎo)致未發(fā)生滑坡的區(qū)域內(nèi),預(yù)期滑坡的面積增大,誤判率升高。log(q/T)≤-3.1時(shí),量化指數(shù)的SI/EI值最大,作為預(yù)測(cè)滑坡判別標(biāo)準(zhǔn)是較恰當(dāng)?shù)?滑坡預(yù)測(cè)與實(shí)際情況有較高的重合度,保持較低的誤判率。以此為標(biāo)準(zhǔn)時(shí),滑坡預(yù)測(cè)的捕獲率為62.50%,誤判率為17.79%。

5 結(jié)論

(1)研究區(qū)土體的抗剪強(qiáng)度參數(shù)符合正態(tài)分布,c值表現(xiàn)出較大的空間差異,變異系數(shù)54.0%～69.2%,而φ值變化較小,變異系數(shù)13.4%～21.9%。3種巖性分布區(qū)土體的抗剪強(qiáng)度參數(shù)不同,粉砂巖地區(qū)的土體c值較高,φ值較小;花崗巖地區(qū)的土體c值較低,φ值較大。

(2)模型預(yù)測(cè)的不穩(wěn)定區(qū)域較連續(xù),主要位于斜坡的下部及深切河谷的兩側(cè),且坡度陡峭的地區(qū),地形坡度以及土體厚度對(duì)模擬結(jié)果影響較大。無條件穩(wěn)定的區(qū)域占研究區(qū)總面積的68.10%,主要為山體頂部、土體厚度薄及地形平坦的區(qū)域。

(3)量化指數(shù)分析顯示隨著log(q/T)水平的提高,模型預(yù)測(cè)的滑坡區(qū)域逐漸擴(kuò)大,SI值逐漸增大, EI值同時(shí)增大, SI/EI值由4.9減小到2.2。單純的提高log(q/T)水平,隨著滑坡與預(yù)測(cè)滑坡之間的空間重合度增大,會(huì)導(dǎo)致預(yù)期滑坡的面積增大,誤判率升高。

(4)以log(q/T)≤-3.1作為預(yù)測(cè)滑坡的判別標(biāo)準(zhǔn),量化指數(shù)的SI/EI值最大,預(yù)測(cè)效果較好,真實(shí)的捕獲率為62.50%,誤判率較低,為17.79%。