庫岸滑坡地質(zhì)災(zāi)害三維演變動態(tài)顯示方法

梁 欣，安 如，張發(fā)明，李 勇，藍秋萍

（1.河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098）

水庫庫岸滑坡是一種與水庫安全運營密切相關(guān)的地質(zhì)災(zāi)害。庫岸滑坡分布于水庫兩側(cè)，穩(wěn)定性易變，一旦坍塌，破壞力巨大，可能造成河道堵塞及道路、橋梁損壞；其伴隨產(chǎn)生的巨型涌浪，影響河道交通安全，甚至影響到對岸居民的生命安全。因此，水庫庫岸滑坡作為一項重大的地質(zhì)工程及水電工程問題，一直是國內(nèi)外學(xué)者廣泛研究的焦點問題。

滑坡預(yù)測分析與建模模擬方面，Petley等對意大利一大型周期性活動滑坡的失穩(wěn)模式進行了預(yù)測分析[1]；Lateltin等對瑞士國家滑坡進行了統(tǒng)一類別劃分，并給出了相應(yīng)的管理辦法[2]；Xie等運用蒙特卡羅方法成功搜索最危險滑動面，建立了滑坡三維穩(wěn)定性分析模型[3]；Tacher等探討了在非均質(zhì)地質(zhì)條件下，滲透率與孔隙水壓力對滑坡穩(wěn)定性的影響[4]；Corominas在GIS技術(shù)支持下，從水位條件影響因子的角度對西班牙Vallcebre滑坡體穩(wěn)定性進行了研究[5]；殷坤龍等運用非連續(xù)變形方法(DDA)對新灘滑坡的動態(tài)演變過程進行了模擬[6]。

滑坡涌浪數(shù)值計算及動畫仿真方面，Geist等成功實現(xiàn)了滑坡災(zāi)害的三維演變過程模擬[7]；Marcello對滑坡涌浪的數(shù)值分析進行了三維研究[8]；潘家錚提出了滑坡入水初始涌浪的計算方法[9]。涌浪模擬幾乎全部建立在理想化的數(shù)值計算基礎(chǔ)上[7-13]，動畫仿真多數(shù)建立在水波模擬的基礎(chǔ)上。隨著研究的深入，涌浪的動畫模擬方法及模擬效果不斷改善[13]。

在以往的滑坡三維演變過程動態(tài)模擬中，大多數(shù)采用VC結(jié)合OpenGL的方法實現(xiàn)，其基本思想是直接圖形法，在數(shù)據(jù)操作中，往往以循環(huán)為基礎(chǔ)展開，大大增加了數(shù)據(jù)操作的復(fù)雜性；在成圖進行局部更新操作時，更體現(xiàn)出其可控制性不高的特點。祝文化等在IDL的支持下，以三峽庫區(qū)榨坊坪滑坡體為例，詳細闡述了滑坡體地面模型和滑動面三維模型的構(gòu)建與可視化顯示方法[14]，但尚缺少動態(tài)集成式的仿真模擬。本文采用IDL面向數(shù)組矩陣運算以及面向?qū)ο笾茍D的基本思想，對三維場景中滑坡滑落、滑坡入水激起涌浪以及涌浪傳播3個過程的動態(tài)集成顯示進行了研究，簡化了數(shù)據(jù)更新操作的難度，提高了對象的可控制性，各過程之間銜接平穩(wěn)，整體穩(wěn)定性高。

1 滑坡演變動態(tài)模擬機理

水庫滑坡失穩(wěn)運動過程是極其復(fù)雜的，因此，只能利用已知的數(shù)據(jù)結(jié)合理論對滑坡滑落全過程進行模擬，而難以做到完全仿真。

在本研究中，可以確定的滑坡要素主要有滑坡中心坐標(biāo)、滑坡沿河長度、滑坡體主滑方向、滑坡體主滑方向長度、滑坡前緣高程、滑坡后緣高程、滑坡估算體積；水位要素主要有水庫蓄水位高度以及滑坡在某蓄水位穩(wěn)定性級別；涌浪要素主要有滑坡入水產(chǎn)生的涌浪高度。

本研究中，需要建立三維場景、滑坡和河道水面3個模型，并利用數(shù)組運算功能獲取指定區(qū)域的數(shù)據(jù)，在面向?qū)ο笾茍D思想的指導(dǎo)下，根據(jù)場景中的坐標(biāo)等屬性特點將3個模型依次添加進三維虛擬現(xiàn)實空間中。建立三維場景后，通過實時更新3個模型的屬性來實現(xiàn)滑坡滑落、激起涌浪、涌浪傳播等動態(tài)過程。滑坡三維演變過程的模擬設(shè)計思路如圖1所示，從左側(cè)啟動（場景重置）開始。

圖1 水庫滑坡滑落過程總體設(shè)計思路

2 滑坡演變動態(tài)過程場景的構(gòu)建

IDL對象圖形系統(tǒng)[15]是一個內(nèi)建的對象類庫，通過選擇適當(dāng)?shù)膬?nèi)建對象、構(gòu)造層次結(jié)構(gòu)來創(chuàng)建圖形場景。圖形對象是功能性的封裝，即每一個對象都包含自己獨立的屬性和方法，圖形創(chuàng)建之后，數(shù)據(jù)及屬性駐留內(nèi)存之內(nèi)，便于修改和反復(fù)使用。圖2為本文對象層次結(jié)構(gòu),從下到上依次加入到上一級容器中。

圖2 對象層次結(jié)構(gòu)

2.1 三維地形對象的構(gòu)建

地形數(shù)據(jù)的采集和輸入是生成DEM模型的關(guān)鍵和基礎(chǔ)。本文中地形數(shù)據(jù)的采集是根據(jù)地形等高線進行，通過二維等高線的矢量化，經(jīng)過ArcGIS的處理獲得 DEM數(shù)據(jù)。紋理映射是建立逼真三維地形場景的重要手段。地形紋理數(shù)據(jù)采用與DEM相同區(qū)域的高分辨率2.5 m SPOT5影像數(shù)據(jù)。

呂秋靈等采用層次細節(jié)方法和多分辨率模型方法，對三維地形的真實性研究做了闡述[16]。在實體建模中，IDL支持利用子塊作為三維單元實體進行整個地形的三維顯示，本文通過對DEM數(shù)據(jù)進行數(shù)據(jù)平滑計算后，運用GRID表面模型建模，建立三維地形。

2.2 河流水面與滑坡體對象的構(gòu)建

河流水面層的繪制比較簡單，只需要根據(jù)DEM數(shù)據(jù)和水位現(xiàn)狀信息將等高面提取出來，然后再進行紋理映射即可。

滑坡體的繪制則要根據(jù)具體的滑坡中心點坐標(biāo)、滑坡沿河長度、滑坡主滑方向長度、滑坡主滑方向4個要素來進行。由于滑坡實測數(shù)據(jù)的難得性，因此，在對滑坡4個地理要素（如圖3所示）進行分析后，將滑坡形狀在外部應(yīng)用程序（如PhotoShop、Windows自帶繪圖工具等）中繪制出來（圖3右下中間紅色區(qū)域為滑坡部分，背景色設(shè)置為黑色），另存為圖片，作為滑坡形狀文件輸入應(yīng)用程序。在應(yīng)用程序中，將滑坡中心點坐標(biāo)、滑坡主滑方向、滑坡主滑方向長度、滑坡沿河長度等作為基本參數(shù)輸入。在讀取滑坡形狀數(shù)據(jù)的同時，利用數(shù)組運算方法進行滑坡形狀重構(gòu)、滑坡在視圖中的無損旋轉(zhuǎn)等，再根據(jù)滑坡中心坐標(biāo)，獲取滑坡區(qū)域所對應(yīng)的DEM數(shù)據(jù)。

圖3 滑坡地理要素示意圖

2.3 光照效果的添加

人類對光線強弱變化的反應(yīng)，要比顏色變化靈敏度更高，添加光照在很大程度上增強了模型的立體感。因此，在三維場景繪制的最后添加光照。圖4所示為三維場景構(gòu)建結(jié)果。

圖4 三維場景構(gòu)建示意圖

3 滑坡失穩(wěn)運動過程數(shù)值計算

3.1 滑坡運動速度計算

國內(nèi)外對于滑坡運動的研究有很多[17]，對滑坡運動速度的計算所采用的方法也很多，其中，應(yīng)用最多的是把滑坡體作為斜面的運動質(zhì)點，進而進行滑坡運動速度的計算。根據(jù)牛頓運動定律，滑坡運動方程為：

式中，為滑坡的質(zhì)量；g為重力加速度；α為滑面傾角；f為動摩擦系數(shù)，f=tanφ；φ為動摩擦角；u為孔隙水壓力；c為其他阻力，如粘聚力等。經(jīng)推導(dǎo)：

式中，v為滑速；v0為初速；s為滑坡（質(zhì)心）沿滑坡的滑距；H為滑坡（質(zhì)心）的垂直降距；L為滑坡（質(zhì)心）水平運動距離。

在實際運用中，H的數(shù)值近似采用滑坡中心的垂直降距，孔隙水壓力和其他阻力等實測數(shù)據(jù)是很難得到的，在應(yīng)用程序計算中，通常忽略u和c的影響，而改用f來綜合反映，因此，提高了對f值的精確度要求。

3.2 滑坡入水激起涌浪高度的計算

潘家錚于1980年提出了初始浪高的計算方法[9]：當(dāng)岸坡發(fā)生水平運動時，初始涌浪可表示為

當(dāng)岸坡發(fā)生垂直運動時，初始涌浪可表示為

式中，ξ0為激起的涌浪初始高度（單位：m）；h為水庫平均深度（單位：m）；v為岸坡水平運動速度（單位：m/s）；v″為岸坡垂直運動速度（單位：m/s）；g為重力加速度（單位：m/s2）。

f函數(shù)可表示如下：

3.3 涌浪傳播數(shù)值計算

貝塞爾方程是在圓柱坐標(biāo)（球坐標(biāo)）下使用分離變量法求解拉普拉斯方程和亥姆霍茲方程時得到的，因此貝塞爾方程在波動問題以及各種涉及有勢場的問題中占有重要地位。本文涌浪傳播過程中的數(shù)值計算采用第一類貝塞爾方程獲得。第一類α階貝塞爾方程Jα（x）是貝塞爾方程當(dāng)α為整數(shù)或α非負時的解。

當(dāng)α為整數(shù)時，貝塞爾的一種積分式為：

另一種積分表達式為：

根據(jù)貝塞爾方程可以繪制出曲線、二維平面圖以及涌浪傳播的三維靜態(tài)數(shù)值模擬圖，如圖5所示。

圖5 貝塞爾方程曲線、二維貝塞爾數(shù)值平面、三維貝塞爾靜態(tài)波動面（左圖曲線100個點，中圖512×512，右圖為三維波動面數(shù)值靜態(tài)圖）

4 滑坡穩(wěn)定性與模擬應(yīng)用分析

4.1 滑坡穩(wěn)定性分析

本文所選擇模擬對象為瀾滄江流域較大滑坡體，位于小黑江右岸，滑坡力學(xué)性質(zhì)為牽引式，滑坡體物質(zhì)組成以粘土和粉土夾碎塊石為主，估算體積達到5 800 000 m3，滑坡深度約為20 m，滑坡沿河長度為640 m，主滑方向長度為320 m，屬于古滑坡類型，河谷地貌類型屬于U型谷。雖然天然狀況下處于穩(wěn)定狀態(tài)，但在蓄水過后，大部分將被淹沒，而且距離大壩較近，在蓄水位從812 m驟降到765 m時，滑坡極易發(fā)生，入水涌浪高度經(jīng)估算達到40.35 m?；乱坏┌l(fā)生，到達大壩水位高度為1.28 m，涌浪沖擊波將影響到大壩的安全運營。因此，在配合滑坡數(shù)值分析的情況下，本文利用已有數(shù)據(jù)對滑坡滑落過程進行三維可視化模擬，將二維數(shù)據(jù)直觀表現(xiàn)出來，同時也可配合管理人員進行滑坡體的監(jiān)測與管理。

水位驟降時，滑坡體內(nèi)部孔隙水向坡體外滲流，導(dǎo)致土體強度參數(shù)降低，穩(wěn)定系數(shù)降低。Fredlund[18]和張文杰[19]在對邊坡穩(wěn)定性影響的研究中，分別由飽和滲透系數(shù)、土-水特征曲線和極限平衡法得到了水位水文條件對邊坡穩(wěn)定性的影響規(guī)律。用摩根斯坦-普萊斯法計算該滑坡體穩(wěn)定系數(shù)為0.93，畢肖普法計算穩(wěn)定系數(shù)為0.91。由此可以得出，在水電站不斷營運的過程中，在水位不斷變化的情況下，滑坡極易失穩(wěn)滑落。圖6為水位由812 m驟降到765 m時的局部穩(wěn)定性示意圖，綠色區(qū)域為不穩(wěn)定區(qū)域。

圖6 水位由812 m驟降至765 m局部穩(wěn)定性示意圖

4.2 模擬應(yīng)用

在實際滑坡入水過程中，滑坡入水產(chǎn)生連續(xù)不斷的對水面的沖擊力而激起涌浪，因此在考慮到運行效率以及真實性的前提下，本文采用將涌浪激發(fā)點設(shè)置在滑坡主滑方向線上，采用10次觸發(fā)涌浪的方法制造連續(xù)涌浪，以此創(chuàng)造出涌浪不斷被激起的起伏效果，提高模擬的真實性。圖7為滑坡動態(tài)演變部分截圖。

圖7 滑坡三維演變動態(tài)局部截圖

4.3 模擬分析

在滑坡數(shù)據(jù)量選擇時，既要將滑坡整體區(qū)域全部截取出來，又要考慮到場景構(gòu)建及滑坡數(shù)據(jù)運算的效率，因此本文選擇的模擬數(shù)據(jù)量網(wǎng)格數(shù)量為1 024×1 024，實地距離范圍為2 560 m×2 560 m。電腦硬件配置為Inter(R) Pentium(R) 2.00 GHz，內(nèi)存為1.00 GB?；逻^程發(fā)生時，每一次的滑坡繪制時間及水面繪制時間如圖8（橫軸為統(tǒng)計點個數(shù)Number，縱軸為時間Second）。

圖8 滑坡繪制時間及水面繪制時間統(tǒng)計結(jié)果

通過圖8可得到，滑坡體單次繪制平均時間為0.029 s，標(biāo)準(zhǔn)偏差為0.057，水面數(shù)據(jù)單次繪制平均時間為0.17 s，標(biāo)準(zhǔn)偏差為0.024 s。從整體的運行視覺效果觀察（見圖7），整個動態(tài)過程能夠平穩(wěn)銜接起來，并且在有一定的時間間隔的情況下，滑坡滑落及涌浪的傳播2個連續(xù)不斷的過程能夠更真實地表現(xiàn)出來。在將整個對象按照圖2規(guī)則劃分為多個模型的基礎(chǔ)上,使用數(shù)組矩陣運算，簡化了模型數(shù)據(jù)更新的難度，數(shù)組元素的調(diào)用更直觀、靈活[15]。圖7的模擬結(jié)果顯示，整個模擬過程運行穩(wěn)定性高；圖8的時間統(tǒng)計結(jié)果說明了各子過程之間能夠平穩(wěn)過渡。結(jié)果表明，整個滑坡動態(tài)演變過程之間銜接平穩(wěn)，各對象數(shù)據(jù)操作簡便，對象的可控制性高，整體穩(wěn)定性好。

5 結(jié) 語

本文以瀾滄江糯扎渡水電站某大型滑坡為例，通過數(shù)值分析及現(xiàn)實計算，成功實現(xiàn)了滑坡失穩(wěn)滑落、滑坡入水激起涌浪、涌浪在河道傳播三者的無縫結(jié)合，較好地完成了滑坡滑落全過程的動態(tài)模擬。此次研究為后期運用DDA方法進行滑坡三維演變的動態(tài)顯示提供了有力的支持。面相對象制圖思想的成功應(yīng)用，滿足了“應(yīng)對變化，提高復(fù)用”的要求，降低了數(shù)據(jù)處理與制圖具體算法之間的耦合度，不僅使程序本身具有更高層次的抽象度，而且具有高度的易擴充性和安全性。本文研究表明，程序文本層次結(jié)構(gòu)簡單分明，理論與實際操作耦合度高，方法應(yīng)用靈活性高。

[1]Petley D N, Mantovani F, Bulmer M H, et al. The Use of Surface Monitoring Data for the Interpretation of Landslide Movement Patterns[J]. Geomorphology, 2005，66(1): 133-147

[2]Lateltin O, Haemmig C, Raetzo H, et al. Landslide Risk Management in Switzerland[J]. Landslides, 2005，2(4): 313-320

[3]Xie M, Esaki T, Zhou G Y. Geographic Informationsystemsbased Three-dimensional Critical Slope Stability Analysis and Landslide Hazard Assessment[J].Journal of Geotechnical and Geoenvironmental Engineering, 2003, 129(12): 110-118

[4]Tacher L, Bonnard C, Laloui L, et al. Modelling the Behaviour of a Large Landslide with Respect to Hydrogeological and Geomechanical Parameter Heterogeneity[J]. Landslides, 2005，2(1): 3-14

[5]Corominas J,Santacana N. Stability Analysis of the Vallcebre Translational Slide, Eastern Pyrenees (Spain) by Means of a GIS[J].Natural Hazards,2003，30(3): 473-485

[6]殷坤龍, 姜清輝, 汪洋.新灘滑坡運動全過程的非連續(xù)變形分析與仿真模擬[J].巖石力學(xué)與工程學(xué)報,2002，21(7): 959-962

[7]Geist E L, Jakob M, Wieczorek G F,et al. Preliminary Assessment of Landslide-induced Wave Hazards, Tidal Inlet, Glacier Bay National Park, Alaska[R].Department of the Interior and US Geological Survey,2003

[8]Risio M, Girolamo P, Bellotti G, et al. Three-dimensionalExperiments on Landslide Generated Waves at a Sloping Coast[J].Coastal Engineering, 2009，56(5):659-671

[9]潘家錚.建筑物的抗滑穩(wěn)定和滑坡分析[M].北京:水利出版社,1980

[10]Fritz H M, Hager W H, Minor H E. Landslide Generated Impulse Waves 2: Hydrodynamic Impact Craters[J]. Experiments in Fluids, 2003，35(6):520-532

[11]周劍華.水庫滑坡涌浪災(zāi)害的數(shù)值研究[J].長江科學(xué)院院報,2003，20(2):7-9

[12]李未,王如云,張長寬.滑坡涌浪的產(chǎn)生與傳播波形分析與計算[J].水科學(xué)進展, 2004， 15(1): 45-49

[13]寧德志.快速多極子邊界元方法在完全非線性水波問題中的應(yīng)用[D].大連:大連理工大學(xué),2005

[14]祝文化, 田金華, 池秀文，等.IDL支持下的滑坡可視化方法研究[J].武漢理工大學(xué)學(xué)報,2004，26(8): 55-56

[15]韓培友.IDL可視化分析與應(yīng)用[M].西安:西北工業(yè)大學(xué)出版社,2006

[16]呂秋靈, 張俊霞. 三維地形可視化及其實時顯示方法[J]. 河海大學(xué)學(xué)報: 自然科學(xué)版, 2002，30(4): 85-87

[17]汪洋, 殷坤龍.水庫庫岸滑坡速度及其涌浪災(zāi)害研究[D].武漢:中國地質(zhì)大學(xué), 2005

[18]Fredlund D G, Xing A, Fredlund M D, et al. The Relationship of the Unsaturated Soil Shear Strength Function to the Soil Water Characteristic Curve[J]. Canadian Geotechnical Journal,1995，32:440- 448

[19]張文杰,詹良通,凌道盛,等,水位升降對庫區(qū)非飽和土質(zhì)岸坡穩(wěn)定性的影響[J].浙江大學(xué)學(xué)報:工學(xué)版, 2006，40(8):1 365-1 370