王小東，馬瀟菲，袁廣祥

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

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

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

1 三維簡(jiǎn)化畢肖普法的基本原理

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

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

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

(1)

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

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

ma=cosq+tanjsina/F

(2)

重力W的計(jì)算與二維不同，即

W=Vγt

(3)

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

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

(4)

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

2 基于DEM的三維畢肖普法實(shí)現(xiàn)的方法和步驟

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

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

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

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

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

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

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

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

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

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

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

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

(8)結(jié)果輸出與表現(xiàn)。采用三維畢肖普法計(jì)算安全系數(shù)以及對(duì)應(yīng)的最危險(xiǎn)滑動(dòng)面，并在三維空間上進(jìn)行展示。

3 算法實(shí)現(xiàn)與實(shí)例分析

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

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

沿滑動(dòng)方向繪制主剖面線，參照瑞典條分法實(shí)現(xiàn)過程中構(gòu)建試算滑動(dòng)面的方法，構(gòu)建試算圓弧滑動(dòng)面，采用畢肖普法獲得最小安全系數(shù)對(duì)應(yīng)的最危險(xiǎn)滑動(dòng)面的位置，并記下圓心坐標(biāo)和半徑，如圖4所示。這一過程和瑞典條分法的實(shí)現(xiàn)過程一致。

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

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

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

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

圖5 三維球形滑動(dòng)面的構(gòu)造過程示意圖Fig.5 Construction process of 3D spherical sliding surface

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

圖6 構(gòu)造三維滑動(dòng)面并生成滑動(dòng)面DEM表面Fig.6 Constructing 3D sliding surface and generating DEM surface of sliding surface

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

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

圖7 網(wǎng)格與坡面和滑動(dòng)面的關(guān)系Fig.7 The relationship between mesh and original slope DEM and spherical sliding surface DEM

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

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

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

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

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

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

4 結(jié)論

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

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

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

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