基于DEM的三維微條柱法計算滑移型塌岸安全系數的研究

王小東，馬瀟菲，袁廣祥

(華北水利水電大學， 河南 鄭州 450011)

目前邊坡穩(wěn)定性分析，一般簡化為二維平面應變問題并采用二維極限平衡法來處理，二維條分法是對邊坡實際情況簡化后采用的方法，實際上，任何邊坡失穩(wěn)現象均是在現實的三維空間中發(fā)生的，因此進行三維計算與分析是研究的方向，即將邊坡置于三維空間之中，在傳統(tǒng)的二維條分法的基礎上進行拓展，二維條分法中的土條也相應地拓展為三維條柱法中的條柱，通過分析每個條柱上的受力狀況，根據極限平衡理論，最終計算得到三維安全系數。三維極限平衡法發(fā)展至今，主要的理論模型有三維簡化畢肖普法、三維簡化簡布法、三維斯賓塞法以及三維Morgenstern-Price法等[1-6]。三維極限平衡法計算結果與二維極限平衡法的計算相比，大多偏于安全[7]，也正因此，三維極限平衡法需要進一步研究、發(fā)展和完善。二維極限平衡分析法的計算機實現程序已相對成熟，Simon等[8]基于excel的VBA開發(fā)模式，融合了三種極限平衡計算公式，實現安全系數的計算；一些商業(yè)化的軟件如Geoslope等廣泛應用于邊坡穩(wěn)定性計算，并被廣大的研究者和工程師所接受。在三維極限平衡分析法的計算機程序實現中，Reid等[9]在研究構造火山側翼邊坡三維穩(wěn)定性分析中，采用三維畢肖普法與DEM(Digital Elevation Model，簡稱DEM)相結合，進行邊坡穩(wěn)定性評價；姜清輝等[10]通過多層DEM構造三維地質模型，基于Visual C++可視化開發(fā)平臺，結合OpenGL，開發(fā)了三維邊坡穩(wěn)定性極限平衡分析軟件slope3D。簡化畢肖普法是畢肖普提出的邊坡穩(wěn)定性計算方法，在工程中得到廣泛的應用，其計算方法簡單且計算結果有較高的精度。Hungr于1987年在簡化畢肖普法的基礎上，將該方法在三維上進行了直接的拓展，即三維畢肖普法。其假設條件與二維狀態(tài)時完全一樣[1，11]。此后，許多研究者也做過類似的研究[12-13]，就三維畢肖普法而言，其基本假設和理論基礎均是一致的。

本文在Reid等[9]對構造火山側翼進行三維穩(wěn)定性研究所采用的理論和方法的基礎上，結合高分辨率DEM，探討三維畢肖普法的實現，其優(yōu)勢在于結合了地理信息系統(tǒng)(Geographical Information System，簡稱GIS)軟件平臺，能夠快速實現三維模型的構建，使得三維模型構建與三維安全系數的計算更加的自動化和一體化。

1 三維簡化畢肖普法的基本原理

本研究基于高分辨DEM數據，選擇在二維條分法的應用領域里被工程所接受的畢肖普法進行擴展，即三維畢肖普法，進行三維安全系數的計算。圖1展示在三維狀態(tài)下條柱的全部受力狀況[14]。

基于圖1中條柱的受力狀況，結合Reid等對構造火山側翼進行三維穩(wěn)定性研究所采用的理論和方法，即滿足每個條柱沿豎直方向即Z軸方向的力的平衡，并根據繞X軸的整體力矩平衡求解安全系數，計算三維安全系數的公式為：

圖1 條柱受力狀況Fig.1 Forces acting on individual columns

(1)

其中，R為底滑面抗滑力的力臂；W為條柱重力；Ac為滑動面與條柱相交部分的面積；c為粘聚力；φ為內摩擦角；k為水平地震力加速度系數，其作用點為條柱中點，力臂為e，無地震力作用時，k=0；γu為孔隙壓力比；θ為失穩(wěn)條柱底部的真傾角；α為底滑面沿滑動方向的視傾角。

其中，mα的計算中包括安全系數F，因此，三維畢肖普法的計算也同樣需要迭代運算才能完成，mα采用公式(2)表示：

ma=cosq+tanjsina/F

(2)

重力W的計算與二維不同，即

W=Vγt

(3)

其中，體積V通過計算微條柱的體積得到，通過計算平截頭六面體的體積近似得到，即

V=1/6Δx(S0+4S1+S2)

(4)

該方法由于未考慮x軸和y軸的整體力的平衡，所以根據整體力矩平衡條件求解安全系數時與坐標軸的位置有關，所以該方法比較適合于滑裂面為弧形面的情況。

2 基于DEM的三維畢肖普法實現的方法和步驟

本文提出了基于GIS的三維極限平衡法的實現思路和方法，即參照Reid等的研究[9]，基于所獲取的2.5 m×2.5 m DEM，采用GIS組件開發(fā)模式，完成三維安全系數的計算。

應用三維微條柱法實現三維穩(wěn)定性分析的具體原理及步驟如下：

(1)確定滑移型塌岸邊界。通過野外調查，結合遙感解譯，確定涉水滑坡堆積體范圍，以矢量多邊形的形式在二維空間中表達，如圖2(a)所示。

(2)確定失穩(wěn)滑動的方向?；菩退妒Х€(wěn)滑動方向是進行塌岸體三維條柱劃分的基礎。在塌岸體詳細調查的基礎上，對塌岸可能的失穩(wěn)方向作出判斷。

(3)三維微條柱劃分。以DEM的空間分辨率為參考，對所要進行穩(wěn)定性計算的三維岸坡進行微條柱劃分，即沿滑動方向和與滑動方向垂直的方向生成一系列垂直相交的直線，覆蓋整個塌岸體在二維平面的投影范圍，并與塌岸體邊界線進行相交運算，最終完成對整個塌岸體的條柱劃分，如圖2(b)所示。

圖2 條柱劃分過程Fig.2 The division process of individual columns

(4)滑動面的選擇與確定。首先確定滑移形塌岸體岸坡性質，即判斷塌岸岸段土體為粘性土還是無粘性土，這樣便可確定塌岸體可能的失穩(wěn)滑動模式，為選擇試算滑動面提供依據。

(5)確定滑動角?；瑒用娲_定之后，滑動角也基本能夠確定，對于滑動面為平面的情況，如圖3(a)所示的二維剖面，每一個條柱的滑動角即為滑動面的傾角；對于滑動面為球形或橢球形時，如圖3(b)所示的二維剖面，每一個條柱的滑動角是不同的，即沿著過條柱底部中心點的切平面與水平面的夾角；對于組合滑動面的情況，如圖3(c)所示的二維剖面，即不同的滑動面區(qū)域采用相應的滑動角；當然，除了圖3給出的這三種情況外，滑動面還有很多種情況，但均可以通過平面與平面、平面與球面的組合方式得到。

圖3 滑坡面形狀Fig.3 The shape of Landslide surface

(6)土層劃分與參數載入。由于考慮到水庫蓄水對滑移型岸坡的影響，因此，即使針對均質岸坡而言，也被蓄水位線分割為水上和水下兩個部分，因此要在豎直方向上對土體進行分層，針對不同土層采用不同的物理力學參數，按土層頂部和底部高程確定不同土層采用不同的參數，通過鍵盤輸入或文件導入的方式載入程序，參與運算。

通過(1)—(6)步可完成對塌岸體三維地質模型的構建，下面則針對極限平衡法進行計算和分析。

(7)依據三維簡化畢肖普法計算安全系數。為了計算安全系數，首先按照以上所述(1) —(6)步構建三維塌岸體地質模型。

(8)結果輸出與表現。采用三維畢肖普法計算安全系數以及對應的最危險滑動面，并在三維空間上進行展示。

3 算法實現與實例分析

本研究以Visual Studio為可視化開發(fā)平臺， 嵌入ArcGIS Engine(簡稱AE)地理信息系統(tǒng)二次開發(fā)組件，以高分辨率DEM數據為基礎，基于三維微條柱法實現的基本方法和步驟，以三維簡化畢肖普法為基礎，完成了塌岸體三維地質模型的構建和三維安全系數的計算。以金沙江溪洛渡庫區(qū)老木溝滑移型塌岸體為例，進行岸坡穩(wěn)定性評價，過程如下：

(1)基于二維畢肖普法搜索得到臨界滑動面

沿滑動方向繪制主剖面線，參照瑞典條分法實現過程中構建試算滑動面的方法，構建試算圓弧滑動面，采用畢肖普法獲得最小安全系數對應的最危險滑動面的位置，并記下圓心坐標和半徑，如圖4所示。這一過程和瑞典條分法的實現過程一致。

圖4 二維畢肖普法獲得臨界滑動面Fig.4 Critical slip surface obtained by 2D Bishop method

(2)基于AE組件二次開發(fā)構建球形滑動面

得到二維狀態(tài)下的臨界滑動面所在的圓弧的圓心和半徑之后，就可通過該圓心和半徑構造三維球形滑動面，其過程在AE二次開發(fā)中實現時，最終的目標是獲得沿與主剖面線垂直方向按一定距離間隔平移后得到的一組剖面線在球形滑動面上投影點的三維坐標，如圖5所示。

圖5顯示了在AE地圖容器中構造球面上點的過程，圖5(a)中的間隔距離即為過該剖面線所在球體截面圓與過主剖面線截面圓之間的垂直距離，將其投影到過主剖面線所在的截面圓上，如圖5(a)所示，可得到一系列同心圓，這些圓的半徑可由球半徑和與過主剖面所在截面圓的距離通過勾股定理計算得到，如圖5(b)所示。經過這一過程之后，便構造出一組同心圓弧，將其離散成點之后，計算圓弧滑動面上點的高程，就得到了該點的三維坐標，最終得到一組球面上的三維坐標點，插值即可生成三維球形滑動面。

圖5 三維球形滑動面的構造過程示意圖Fig.5 Construction process of 3D spherical sliding surface

圖6以老木溝堆積體為例，展示了通過AE二次開發(fā)構建三維滑動面的過程，如圖所示，第一步：由二維畢肖普法計算得到主剖面上二維臨界圓弧滑動面；第二步：由二維滑動面構造三維球形滑動面的一組三維點坐標，三維球形滑動面在平面上的投影為一橢圓；第三步，由三維點坐標通過插值生成三維球形滑面DEM。

圖6 構造三維滑動面并生成滑動面DEM表面Fig.6 Constructing 3D sliding surface and generating DEM surface of sliding surface

(3)基于三維畢肖普法計算得到三維安全系數

三維滑動面構建完成之后即可按照三維畢肖普法的計算公式，對岸坡所在區(qū)域的多邊形進行微條柱劃分，如圖7所示，將塌岸多邊形沿滑動方向進行正方形網格剖分，分別投影在原始坡面DEM和三維滑動面DEM之上，以便于進行每個條柱上的計算，根據公式(1)進行三維安全系數的計算，與二維畢肖普法一樣，三維畢肖普法的實現同樣是一個迭代的過程，即采用公式(1)和公式(2)完成這一過程。

圖7 網格與坡面和滑動面的關系Fig.7 The relationship between mesh and original slope DEM and spherical sliding surface DEM

對于式(1)—式(4)中所包含θ，α和V等參數，在AE組件開發(fā)中可采用下列方法確定：

真傾角θ的確定：真傾角θ與真傾向有關，是指條柱底部與三維球形滑動面相交的面的中心所在的切平面與水平面的夾角，從幾何意義上講，是真傾向與其水平面上投影線的夾角，真傾向可通過對三維滑動面DEM進行坡向分析得到，如圖8中坡向圖所示。真傾角則可通過對三維滑動面DEM進行坡度分析得到，如圖8中坡度圖所示，在二次開發(fā)進行三維畢肖普法計算時，基于條柱底面中心點位置坐標和坡度圖獲得該條柱底面的真傾角θ。

圖8 三維球形滑動面坡向圖和坡度圖Fig.8 Aspect map and slope map of 3D spherical sliding surface

視傾角α的確定：視傾角α與真傾向和視傾向有關，對某個條柱底面而言，存在多個視傾角，但對于三維畢肖普法計算安全系數而言，視傾角α是指沿岸坡主滑方向的視傾角，即將主滑方向作為視傾向，通過視傾向與真傾向之間的夾角ω，如圖9所示，可通過公式tanθ=tanα·cosω，計算得到視傾角α。

圖9 真傾角與視傾角之間的關系Fig.9 The relationship between true dip angle and apparent dip angle

最終，以金沙江溪洛渡庫區(qū)老木溝滑移型塌岸體為例進行三維安全系數的計算，其在天然狀態(tài)下的抗剪強度指標為：c=25 kPa，φ=40°，γ=20 KN/m3，以25 m間隔構建條柱，通過三維畢肖普法計算得到的三維安全系數為1.48，而采用二維畢肖普法計算得到安全系數為1.04。由此可見，三維畢肖普的計算結果偏于安全。引起這種差異的原因較多，其中一個原因是由滑動面的形狀造成的，球形滑動面在現實中很少見到，圖7也明確顯示出球形滑動面在平面的投影范圍沒有完全覆蓋實際調查中所確定的塌岸體的范圍，因此，在理論上是一種可行的方法，在實際中是一種近似的方法。

4 結論

本文通過AE二次開發(fā)的方式，將三維畢肖普法引入GIS軟件平臺，利用GIS平臺的數據處理與空間分析功能，耦合極限平衡理論計算邊坡安全系數的方法，實現了三維圓弧滑動面的構造和三維安全系數的計算，主要結論如下：

1)以高分辨率DEM為結合點，耦合地理信息系統(tǒng)分析方法與三維極限平衡法，進行安全系數的計算。高分辨率DEM提供了詳細的地表信息，十分適合邊坡的三維模型構建，本文將三維畢肖普法引入GIS軟件平臺，實現了三維安全系數的計算，為單體穩(wěn)定性評價與區(qū)域穩(wěn)定性評價的結合提供了新的思路和方法。

2)提高了邊坡三維模型構建與三維安全系數的計算的自動化、一體化程度。AE二次開發(fā)的方式提高了三維地質模型構建的效率，即通過圖切剖面線方式構建球形滑動面，提供了邊坡三維模型構建的自動化程度，并將GIS中坡度、坡向等空間分析方法與三維極限平衡法相結合，實現了三維安全系數計算的一體化。

3)三維滑動面的構造需要進一步的優(yōu)化。應用本論文方法，以金沙江溪洛渡庫區(qū)老木溝滑移型塌岸體為例，進行三維穩(wěn)定系數的計算，與二維計算結果相比，結果偏于安全。主要的原因在于將滑動面視為球形滑動面，是一種理想的狀態(tài)，實際的滑動面大多是非圓弧的，因此需要在圓弧滑動面的基礎上進行優(yōu)化，可在程序中加入控制性結構面或控制性節(jié)點，使其與實際情況更加接近。