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

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

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

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

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

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

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

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

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

1 滑坡體模型的建立

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

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

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

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

圖2 滑坡體Fig.2 A landslide mass

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

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

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

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

其中：ψi－1＝cos（θxzi－1－θxzi）－sin（θxzi－1－θxzi）tgφi／Fs3D 為傳遞系數(shù)；θA為方位角 ，θxz為滑面與X軸的夾角；Wi為各滑塊重力；Ai為各滑面面積.可以計算出X軸滑動方向的整體的穩(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.式中各項的值可從柵格數(shù)據(jù)中取值得到.

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

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

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

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

4 滑動方向的確定

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

5 工程實例

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

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

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

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

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

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

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

6 結(jié) 論

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

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

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

［3］吳秀芹.ArcGIS 9地理信息系統(tǒng)應(yīng)用與實踐［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ù)字地形圖上滑坡動態(tài)方向穩(wěn)定性系數(shù)的計算與滑坡滑動方向的確定［J］.巖石力學(xué)與工程學(xué)報，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

浙江省科技廳基金資助項目（2009C33047）

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

（

劉 巖）