基于GIS的邊坡三維穩(wěn)定性計(jì)算

楊 偉，丁伯陽，潘曉東，朱益軍，陶海冰

（1.浙江工業(yè)大學(xué) 巖土工程研究所，浙江 杭州 310032；2.浙江省交通規(guī)劃設(shè)計(jì)研究院，浙江 杭州 310014；3.杭州師范大學(xué) 理學(xué)院，浙江 杭州 310036）

基于GIS的邊坡三維穩(wěn)定性計(jì)算

楊 偉1，丁伯陽1，潘曉東1，朱益軍2，陶海冰3

（1.浙江工業(yè)大學(xué) 巖土工程研究所，浙江 杭州 310032；2.浙江省交通規(guī)劃設(shè)計(jì)研究院，浙江 杭州 310014；3.杭州師范大學(xué) 理學(xué)院，浙江 杭州 310036）

根據(jù)高速公路高邊坡的鉆孔資料，利用Arc Map建立進(jìn)行三維邊坡穩(wěn)定性分析的TIN數(shù)據(jù)層和柵格數(shù)據(jù)層，將柵格數(shù)據(jù)模型的各層對應(yīng)柵格塊視為一個(gè)滑坡單元體，通過讀取相應(yīng)柵格值進(jìn)行穩(wěn)定性計(jì)算.針對順層滑坡的特點(diǎn)，采用ArcGIS集成的VBA開發(fā)環(huán)境，提出了一種三維折線法分析邊坡的穩(wěn)定性.該方法通過對柵格數(shù)據(jù)的重采樣，可以計(jì)算任意滑動(dòng)方向三維邊坡的安全系數(shù)，通過選擇最小安全系數(shù)來確定邊坡的主滑動(dòng)方向.結(jié)合某高速公路大型巖質(zhì)滑坡的工程實(shí)例，通過比較三維折線法計(jì)算結(jié)果與某典型剖面的二維不平衡推力法計(jì)算結(jié)果，討論該方法對工程滑坡治理的指導(dǎo)意義.

三維邊坡；邊坡穩(wěn)定；GIS；滑動(dòng)方向；折線法

GIS（地理信息系統(tǒng)）軟件在土木工程中的應(yīng)用越來越廣泛，在公路邊坡這種涉及大量空間信息存儲、處理和分析的工程實(shí)踐中，更加體現(xiàn)出其優(yōu)越性.目前，滑坡穩(wěn)定性分析常用二維的極限平衡法，對某一可能滑體，采用計(jì)算幾個(gè)不同剖面的穩(wěn)定系數(shù)來綜合評價(jià).雖然出現(xiàn)了很多三維極限平衡法，多以柱體單元為計(jì)算單元，由于計(jì)算方法繁瑣，在二維GIS中較難實(shí)現(xiàn)，通常只借助GIS作為空間信息存儲和處理的工具，通過數(shù)據(jù)的轉(zhuǎn)換實(shí)現(xiàn)在其他軟件中的三維穩(wěn)定性計(jì)算.對于邊坡三維問題，謝漠文結(jié)合GIS柵格數(shù)據(jù)和4個(gè)邊坡穩(wěn)定極限平衡模型開發(fā)了一個(gè)用于邊坡三維穩(wěn)定分析的GIS擴(kuò)展模塊，提供了一種便利的數(shù)據(jù)處理方法和分析方法［1］.

利用ArcGIS的VBA開發(fā)環(huán)境，采用工程中廣泛應(yīng)用的折線法，進(jìn)行滑坡三維穩(wěn)定性計(jì)算，可以方便地確定滑坡體的三維穩(wěn)定系數(shù)，并對滑坡體的假定的可能滑動(dòng)方向進(jìn)行試算，找出最有可能的滑動(dòng)方向.

1 滑坡體模型的建立

滑坡體的DEM模型和各土層的TIN數(shù)據(jù)、柵格數(shù)據(jù)，可以根據(jù)鉆孔資料，通過合適的插值方法得到.各土層的上下關(guān)系通過柵格法進(jìn)行空間拓?fù)潢P(guān)系的處理［2］.ArcGIS提供的空間插值方法有反距離權(quán)重插值、樣條法插值和克里格插值.無論選擇哪種插值方法，樣本點(diǎn)越多，樣本點(diǎn)分布越廣，插值結(jié)果越接近實(shí)際值［3］.而在實(shí)際工程中，鉆孔點(diǎn)的分布可能達(dá)不到計(jì)算精度的要求，這時(shí)可以引入虛擬鉆孔來提高精度.

順層滑坡的滑動(dòng)面，通?？梢酝ㄟ^鉆孔資料得到破碎滑動(dòng)帶來確定，或分別對各地質(zhì)層層底作為潛在滑面，分別做穩(wěn)定計(jì)算而得到.圖1為邊坡DEM模型，圖2為以破碎帶為滑動(dòng)面得到的滑坡體.

圖1 邊坡DEM模型Fig.1 A slope digital elevation model

2 滑坡體三維穩(wěn)定計(jì)算

圖2 滑坡體Fig.2 A landslide mass

對滑坡體進(jìn)行柱狀單元?jiǎng)澐?，由于滑坡體基于柵格數(shù)據(jù)，事實(shí)上可以看成是以柵格大小為面積，以地面到滑面距離為高的柱狀體.與滑坡體穩(wěn)定性分析相關(guān)的數(shù)據(jù)，如：各土層底面高度，水位線高度，滑動(dòng)面傾角，滑動(dòng)面傾向等，存儲在相應(yīng)的柵格數(shù)據(jù)中［4］.在ArcGIS的VBA開發(fā)環(huán)境中，可以在滑坡區(qū)域內(nèi)提取相應(yīng)的柵格數(shù)據(jù)進(jìn)行計(jì)算.

假定滑體的主滑動(dòng)方向?yàn)閄軸，Y軸垂直滑動(dòng)方向.每一個(gè)柵格單元都與X軸垂直.柵格柱體間的條間力可分解為水平向的剪力Q、沿X軸方向的推力Ei和Y軸方向的擠壓力Eh，對于每一個(gè)滑體的受力情況如圖3所示.

圖3 柵格柱狀體受力圖Fig.3 Force on a cell column

采用折線法計(jì)算，忽略滑塊水平方向的剪力的作用，假定Q＝0，滑塊只傳遞主滑動(dòng)方向的推力，不考慮擠壓力的影響，則基于柵格數(shù)據(jù)的滑體向主滑動(dòng)方向的推力［5－6］為

其中：ψi－1＝cos（θxzi－1－θxzi）－sin（θxzi－1－θxzi）tgφi／Fs3D 為傳遞系數(shù)；θA為方位角 ，θxz為滑面與X軸的夾角；Wi為各滑塊重力；Ai為各滑面面積.可以計(jì)算出X軸滑動(dòng)方向的整體的穩(wěn)定系數(shù)為

其中：Ri＝（W cosθi·sinφi＋c Aicosφi）sinθA；Ti＝W sinθi·sinθA.當(dāng)所有滑體的φi值相同，則：Ri＝（W cosθi·tgφi＋c Ai）sinθA；Ti＝W sinθi·sinθA.式中各項(xiàng)的值可從柵格數(shù)據(jù)中取值得到.

3 計(jì)算機(jī)程序?qū)崿F(xiàn)

ArcGIS下的VBA開發(fā)環(huán)境中，提供了眾多柵格數(shù)據(jù)的接口和方法.對于基于同一坐標(biāo)的土層?xùn)鸥駭?shù)據(jù)集合，可以直接讀取柵格值進(jìn)入計(jì)算程序，不必進(jìn)行數(shù)據(jù)轉(zhuǎn)換就在Arc Map中實(shí)現(xiàn)了邊坡三維穩(wěn)定性計(jì)算.與邊坡計(jì)算有關(guān)的各土層的柵格數(shù)據(jù)可以從基于鉆孔數(shù)據(jù)而來的TIN數(shù)據(jù)得到.這樣得到的柵格數(shù)據(jù)層在邊界處會有柵格不重合的現(xiàn)象，即在同一坐標(biāo)點(diǎn)上，有可能只有地面高程數(shù)據(jù)，而沒有滑動(dòng)面高程數(shù)據(jù).在計(jì)算程序中，邊界上無法同時(shí)取到地面數(shù)據(jù)和滑動(dòng)面數(shù)據(jù)的點(diǎn)不讀入數(shù)據(jù).為了計(jì)算簡便，對土層參數(shù)取加權(quán)平均值，作為整個(gè)柵格柱體的參數(shù)，將柵格柱體作為單一土層考慮.

求解安全系數(shù)時(shí)，按以下步驟：1）首先假定安全系數(shù)K＝1，讀取第一行第一列的柵格數(shù)據(jù)，計(jì)算得到R1，T1；2）依次讀取下一列的柵格數(shù)據(jù)，計(jì)算對應(yīng)的傳遞系數(shù)ψi，得到的R′i＝Ri－1ψi＋Ri，T′i＝Ti－1ψi＋Ti繼續(xù)往下傳遞；3）逐行、逐列計(jì)算，按式（2）得到新的安全系數(shù)K′；4）當(dāng)｜K′－K｜＞0.01，則K＝K′，重新進(jìn)行計(jì)算，直到K值滿足要求.計(jì)算框架圖如圖4所示.

圖4 穩(wěn)定系數(shù)計(jì)算框圖Fig.4 Calculating diagram of stability cofficent of landslide

4 滑動(dòng)方向的確定

確定滑坡滑動(dòng)方向的傳統(tǒng)方法是根據(jù)地表和滑動(dòng)面形態(tài)，通過經(jīng)驗(yàn)方法估計(jì)確定，缺乏定量基礎(chǔ)，基于柵格數(shù)據(jù)的邊坡穩(wěn)定性計(jì)算公式中引入了傾角、傾向等參數(shù)，則為滑坡體滑動(dòng)方向的確定提供了條件.選取不同的滑動(dòng)方向計(jì)算出對應(yīng)的邊坡的安全系數(shù)，則安全系數(shù)最小的某個(gè)滑動(dòng)方向，為最可能的滑動(dòng)方向［7］.利用ArcGIS可以較方便的確定主滑動(dòng)方向.在ArcGIS下，柵格數(shù)據(jù)可以通過Rotate工具重采樣得到新的柵格數(shù)據(jù)，提供的重采樣方法有三種：最鄰近采樣（NEAREST），雙線性（BILINEAR）和三次卷積采樣（CUBIC）.最鄰近采樣精度較低，對旋轉(zhuǎn)后得到的柵格數(shù)據(jù)進(jìn)行坡度分析和坡向分析誤差較大.為提高內(nèi)插精度，選用三次卷積采樣，取待計(jì)算點(diǎn)周圍相鄰的16個(gè)點(diǎn)，可先在某一方向上內(nèi)插，如先在X方向上，每四個(gè)值依次內(nèi)插四次，再根據(jù)四次的計(jì)算結(jié)果在Y方上內(nèi)插，最終得到內(nèi)插結(jié)果.對旋轉(zhuǎn)后的滑動(dòng)面數(shù)據(jù)進(jìn)行坡度分析和坡向分析，得到旋轉(zhuǎn)后的坡度角和方位角的柵格數(shù)據(jù)層.此時(shí)得到的柵格單元與原坐標(biāo)系X軸垂直，可直接套用上式計(jì)算X軸滑動(dòng)方向的安全系數(shù)，不必進(jìn)行坐標(biāo)變換.利用上面提到的折線法計(jì)算程序，假定某一方向?yàn)橹骰瑒?dòng)方向，只需通過柵格數(shù)據(jù)的旋轉(zhuǎn)，即可以計(jì)算出不同主滑動(dòng)方向的安全系數(shù).

5 工程實(shí)例

杭金衢K103高邊坡是一個(gè)大型的巖質(zhì)滑坡，該滑體體積約為160萬m3，屬于工程施工誘發(fā)的、潛伏型路塹滑坡.根據(jù)地質(zhì)分析和定量評價(jià)，該滑坡治理前長期處于緩慢蠕滑狀態(tài)，經(jīng)歷公路開挖后4年多緩慢蠕動(dòng)變形，2005年春季，滑坡體后緣及兩側(cè)均發(fā)現(xiàn)多條裂縫，坡腳擋墻出現(xiàn)拉裂、擋墻沉降縫錯(cuò)位、墻面與坡面局部鼓起等現(xiàn)象，相應(yīng)路段高速公路路面也出現(xiàn)了不同程度的裂縫開裂與路面鼓包現(xiàn)象.滑坡方向大致為正東向，與公路延伸方向近于正交.滑坡體的厚度一般在15～40 m之間，局部大于40 m.滑坡巖體破碎，風(fēng)化強(qiáng)烈，呈碎塊石夾泥狀，節(jié)理面風(fēng)化后大多呈泥狀，潛在滑面主要沿破碎帶發(fā)展，滑面分布次生夾泥，強(qiáng)度低.根據(jù)現(xiàn)場地質(zhì)勘察結(jié)論，該滑坡滑面前緩后陡，目前還處于蠕滑變形階段，按照規(guī)范要求，對其取穩(wěn)定系數(shù)K＝0.98，利用Bishop法對該滑面進(jìn)行二維的反演分析，得到滑面參數(shù)的加權(quán)平均值為：c＝20.5，φ＝14.8.

利用高程數(shù)據(jù)生成地面TIN數(shù)據(jù)，再將其轉(zhuǎn)換為柵格數(shù)據(jù).利用鉆孔中的軟弱夾層位置數(shù)據(jù)與滑坡裂隙編輯數(shù)據(jù)生成滑面TIN，也將其轉(zhuǎn)換為柵格數(shù)據(jù).對滑動(dòng)面進(jìn)行坡度分析和坡向分析，得到滑動(dòng)面傾角柵格圖和滑動(dòng)面方位角柵格圖，見圖5.利用VBA開發(fā)的折線法計(jì)算程序計(jì)算安全系數(shù)，當(dāng)c＝20.5，φ＝14.8時(shí)，正東方向的安全系數(shù)為：K＝1.152.對傾向進(jìn)行統(tǒng)計(jì)分析，多分布在60°～120°范圍內(nèi)，在此范圍內(nèi)假定不同主滑動(dòng)方向所得到安全系數(shù)如表1所示.

圖5 坡度角和方位角Fig.5 Slope angle and aspect angle

表1 安全系數(shù)表Table 1 Slope Stability coefficient on different direction

由此推斷，最危險(xiǎn)滑動(dòng)方向在傾向?yàn)?10°，相對應(yīng)的安全系數(shù)為1.096.如果在計(jì)算程序中考慮地下水的作用，安全系數(shù)可能會進(jìn)一步降低.該結(jié)果與實(shí)際邊坡的情況較吻合，有一定的參考價(jià)值.

用不平衡推力法對滑坡進(jìn)行穩(wěn)定性計(jì)算，該法在水利部門、鐵路部門廣泛應(yīng)用.選取滑坡上某一典型剖面，用二維不平衡推力法計(jì)算安全系數(shù)，通過調(diào)整安全系數(shù)的大小，最終使滑坡出口處的剩余下滑推力為0，此時(shí)得到的安全系數(shù)為滑坡的穩(wěn)定系數(shù).同時(shí)計(jì)算地下水位變化時(shí)，滑坡的穩(wěn)定系數(shù)變化情況［8］，見表2.

表2 地下水位變化時(shí)滑坡安全系數(shù)表Table 2 Slope stability coefficient on different groundwater level

可見，在完全疏干的情況下，該典型剖面的安全系數(shù)取值為1.26，與計(jì)算程序得到的滑動(dòng)方向正東10°范圍內(nèi)的結(jié)果比較接近，但并非最小值.而最危險(xiǎn)滑動(dòng)方向的確定，對緩慢發(fā)展的工程滑坡治理具有一定的指導(dǎo)意義.

6 結(jié) 論

基于GIS柵格層數(shù)據(jù)模型能將各層對應(yīng)柵格塊視為一個(gè)滑坡單元體，極大的方便了邊坡三維穩(wěn)定性的計(jì)算.在ArcGIS環(huán)境下，結(jié)合VBA二次開發(fā)技術(shù)，能直接讀取ArcGIS中柵格數(shù)據(jù)的值進(jìn)行穩(wěn)定性計(jì)算，避免了復(fù)雜的數(shù)據(jù)轉(zhuǎn)換工作.基于柵格數(shù)據(jù)的三維折線法邊坡穩(wěn)定性計(jì)算公式，引入了傾角、傾向等參數(shù)，對定量確定滑動(dòng)方向提供了條件，對不同滑動(dòng)方向計(jì)算所得到的安全系數(shù)最小的方向?yàn)樽羁赡艿幕瑒?dòng)方向.不同滑動(dòng)方向的安全系數(shù)相差最大的可到30%，可見滑動(dòng)方向的選取具有一定的工程意義.邊坡三維穩(wěn)定性計(jì)算模型對不同土層的參數(shù)取用的是加權(quán)平均值，相當(dāng)于作為單一土層考慮，且未考慮地下水作用，有待改進(jìn).確定滑動(dòng)方向時(shí)需要對較大范圍進(jìn)行試算，計(jì)算量較大.可以進(jìn)一步探究邊坡的傾向與安全系數(shù)之間的解析式，來確定可能的滑動(dòng)方向，這是下一步要繼續(xù)深入的工作.

［1］謝漠文，江崎哲郎，周國云.基于邊坡單元的三維滑坡災(zāi)害評價(jià)的 GIS方法［J］.巖石力學(xué)與工程學(xué)報(bào)，2003，22（6）：969-976.

［2］陶海冰，蔣蓓.基于數(shù)據(jù)的柵格法模擬三維土層研究［J］.杭州師范大學(xué)學(xué)報(bào)，2007，7（6）：469-473.

［3］吳秀芹.ArcGIS 9地理信息系統(tǒng)應(yīng)用與實(shí)踐［M］.北京：清華大學(xué)出版社，2007.

［4］謝漠文，蔡美峰，江崎哲郎.基于GIS邊坡穩(wěn)定三維極限平衡方法的開發(fā)及應(yīng)用［J］.巖土力學(xué)，2006，27（1）：117-122.

［5］陳斗勇.折線形滑面邊坡安全系數(shù)的一種新的求算法［J］.水文地質(zhì)工程地質(zhì)，1991，18（6）：46-47.

［6］殷宗澤.土工原理［M］.北京：中國水利水電出版社，2007.

［7］李先華，管群，楊武年，等.數(shù)字地形圖上滑坡動(dòng)態(tài)方向穩(wěn)定性系數(shù)的計(jì)算與滑坡滑動(dòng)方向的確定［J］.巖石力學(xué)與工程學(xué)報(bào)，1999，18（3）：308-311.

［8］方鵬飛，朱益軍，朱向榮.杭金衢高速公路K103滑坡治理［J］.工程勘察，2009（1）：31-35.

Three-dimensional slope stability analysis based on GIS

YANG Wei1，DING Bo-yang1，PAN Xiao-dong1，ZHU Yi-jun2，TAO Hai-bin3
（1.Geotechnical Engineering Institute，Zhejiang University of Technology，Hangzhou 310032，China；2.Institute of Communications，Planning，Design ＆Research，Hangzhou 310014，China；3.College of Science，Hangzhou Nomal University，Hangzhou 310036，China）

According to the borehole date of Highway，the TIN data layer and grid data layer for three-dimensional slope stability analysis were carried out by use of Arc Map.The value of each gird cell of different layers are picked up as slide block，it makes convenience in the slope stability analysis.Considering the property of bedding landslide，the three-dimensional slope stability analysis program based on VBA of ArcGIS is put forward in this paper.The program can calculate the Stability factor of any directions by resampling of raster date，thus the real sliding direction can be determined by the minimumValue of the Stability factor.Taking one of the large rocky landslide as an example，the results of the three-dimensional with the two-dimensional method of a typical profile of the slope are compared.

three-dimensional slope；slope stability analysis；GIS；sliding direction；line method

TU457

A

1006-4303（2012）01-0092-04

2010-10-12

浙江省科技廳基金資助項(xiàng)目（2009C33047）

楊 偉（1982—），男，浙江金華人，碩士研究生，主要從事巖土工程數(shù)字化研究，E-mail：airball150336＠126.com.

（

劉 巖）