杜 偉 王佳穎 李玉容 楊國柱
(國網(wǎng)電力空間技術(shù)有限公司, 北京 102209)
近年來我國極端天氣事件時有發(fā)生,受其影響我國多地的滑坡、崩塌和泥石流等地質(zhì)災害事件頻發(fā),由此造成人員傷亡和基礎(chǔ)設(shè)施癱瘓等重大的損失。因此,必須加強對地質(zhì)災害防治的研究,有效保障人民生命財產(chǎn)和國家基礎(chǔ)設(shè)施的安全。為了切實提高國家地質(zhì)災害防治的能力,我國先后頒布了《地質(zhì)災害防治條例》《國務(wù)院關(guān)于加強地質(zhì)災害防治工作的決定》及《全國地質(zhì)災害防治“十三五”規(guī)劃》等法規(guī),計劃建立健全地質(zhì)災害調(diào)查評價、監(jiān)測預警、綜合治理、應急防治四個體系。風險分析包括評估地質(zhì)災害發(fā)生的可能性以及災害發(fā)生的嚴重程度兩個部分,是地質(zhì)災害防治體系中的重要內(nèi)容。[1]
滑坡是一種典型的地質(zhì)災害,具有顯著的不確定性和后果損害性,是典型的風險事件。對滑坡進行風險分析和風險管理的研究是近年來地質(zhì)災害領(lǐng)域的研究熱點和難點。[2]滑坡災害的風險分析可分為區(qū)域和單體尺度。區(qū)域性滑坡風險分析是在大范圍內(nèi)對邊坡進行危險性的時空預測以及對承災體進行易損性的調(diào)查評價,最終繪制危險性區(qū)劃圖,為地區(qū)性的風險管理提供依據(jù);單體滑坡風險分析則側(cè)重于對單一邊坡破壞風險的量化分析,其研究結(jié)果有助于分析滑坡的破壞機理,是在區(qū)域性滑坡風險分析基礎(chǔ)上的細化。[3]
關(guān)于單體滑坡風險的計算,文獻[4-5]提出了不同的計算方法,其中多數(shù)方法須要求解邊坡發(fā)生破壞的概率及災害導致的損失。邊坡的破壞概率與巖土體參數(shù)的不確定性密切相關(guān),可由可靠度理論求得[1]。除了巖土體參數(shù),影響單體滑坡風險分析結(jié)果的重要因素是邊坡的空間幾何形態(tài)、尺寸以及失效后果的估算。通過傳統(tǒng)的測繪方法很難迅速準確地獲得復雜邊坡表面的空間特征,也很難估算潛在的滑坡體積。
基于無人機遙感(UAVRS)獲取邊坡的空間形態(tài)特征是近年來興起的一種新方法。UAVRS系統(tǒng)是指以無人駕駛飛行器為載體,通過掛載相機或激光雷達等遙感傳感器,在無線通信設(shè)備或機載計算機程控系統(tǒng)操控下獲取空間遙感信息的一套系統(tǒng)。UAVRS的優(yōu)點包括設(shè)備成本和使用成本相對較低,使用場景靈活,能夠快速響應和非接觸式等。因此,UAVRS系統(tǒng)已經(jīng)廣泛應用于邊坡的穩(wěn)定性分析、地質(zhì)勘查和變形監(jiān)測等方面。文獻[6]介紹了基于無人機傾斜攝影的公路邊坡三維重建和災害識別方法;文獻[7-8]介紹了使用搭載單鏡頭的小型無人機獲取高陡邊坡的巖體結(jié)構(gòu)面參數(shù)及對露天礦山邊坡建立精細化FLAC3D數(shù)值模型;文獻[9]介紹了通過固定翼無人機監(jiān)測露天采坑地表變形的動態(tài)變化,為滑坡預警預報和應急救援提供數(shù)據(jù)支撐。這些研究表明無人機遙感是獲取復雜坡體表面空間特征的一種新型方法。但是,如何基于無人機遙感進行滑坡災害的定量風險分析仍鮮有報道。
因此,依托無人機遙感平臺,通過激光雷達采集邊坡體的形狀及位置信息,利用ArcMap和Rhino軟件建立邊坡的幾何模型并進行網(wǎng)格劃分;將使用有限差分軟件FLAC3D計算邊坡的安全系數(shù),使用一階可靠度方法計算邊坡的破壞概率;通過對FLAC3D進行二次開發(fā),求解三維復雜邊坡的失效后果。最后進行滑坡災害的定量風險分析。
皖南山區(qū)受地形、氣候、地質(zhì)構(gòu)造和人類工程活動等因素的影響,滑坡災害多發(fā),人民的生命財產(chǎn)安全和公路、鐵路、輸電線路等基礎(chǔ)設(shè)施的平穩(wěn)運行都受到嚴重威脅。因此,選擇皖南山區(qū)池州某縣道旁塔基邊坡作為研究對象,用以研究無人機遙感在滑坡災害風險分析中的應用。該邊坡由碎石土層及下覆的基巖兩部分組成,其土性參數(shù)見表1。
表1 邊坡巖土體參數(shù)Table 1 Soil parameters of slope rock
采用的無人機設(shè)備是大疆經(jīng)緯M600(圖1),其參數(shù)見表2所示。無人機飛行時,搭載數(shù)字綠土LIAIR 220N無人機激光雷達掃描系統(tǒng),用于采集邊坡體的激光點云數(shù)據(jù)。與通過無人機搭載傾斜攝影相機獲取邊坡體的光學影像數(shù)據(jù)相比,對激光點云數(shù)據(jù)處理后可以濾去坡體表面植被的影響,能更真實地反映復雜坡體表面的空間幾何特征。
圖1 大疆經(jīng)緯M600無人機Fig.1 An unmanned aerial vehicle of DJI Matrice 600
表2 無人機參數(shù)Table 2 Technical parameters of UAVs
邊坡的幾何模型包括地表模型及地層模型,結(jié)合二者并進行網(wǎng)格劃分可以建立接近實際情況的坡體精細化數(shù)值模擬網(wǎng)格模型。
許多專家學者已經(jīng)在基于UAVRS采集地表激光點云數(shù)據(jù)并處理得到地表數(shù)字高程模型(DEM)方面做了很多有益的工作[6,8]。基于DEM文件建立邊坡表面(簡稱“坡表”)精細化三維幾何模型,須要用到ArcMap及Rhino軟件。其中,ArcMap軟件主要用于對*.tif格式的DEM文件進行數(shù)據(jù)的光滑性處理,以便消除濾波及影像處理算法誤差引起的DEM畸變,得到相對平滑的坡表幾何形狀,并導出帶有高程信息的點云數(shù)據(jù);Rhino軟件主要用于建立邊坡的實體模型并劃分網(wǎng)格。
2.1.1DEM數(shù)據(jù)的光滑性處理。
在ArcMap軟件中通過“添加數(shù)據(jù)”功能導入*.tif格式的DEM文件;使用“重采樣工具”功能以得到抽稀后的DEM柵格數(shù)據(jù)。圖2為根據(jù)無人機采集的激光點云數(shù)據(jù)得到的邊坡DEM信息。以根據(jù)DEM生成的等高線為例,其重采樣前、后的等高線分別見圖3??梢?抽稀后得到的等高線更明晰平滑。
在上述步驟完成之后,須使用“柵格轉(zhuǎn)點”功能將處理好的DEM柵格數(shù)據(jù)轉(zhuǎn)換成帶有高程的點云數(shù)據(jù),并保存為*.dwg格式的文件。值得注意的是,此時生成的點云數(shù)據(jù)并不含有*.dwg文件可以識別的高程信息,還須要在其“屬性表”里面自行添加高程字段信息。最后,通過“導出至CAD”功能將文件保存為帶有高程信息的*.dwg格式的柵格點云數(shù)據(jù)。
圖2 邊坡的數(shù)字高程信息Fig.2 Digital elevation information of slopes
a—重采樣前的等高線; b—重采樣后的等高線。圖3 重采樣前后的DEM等高線Fig.3 DEM contours before and after resampling
2.1.2邊坡實體模型的建立及網(wǎng)格劃分
在Rhino軟件中導入ArcMap軟件生成的*.dwg格式的柵格點云數(shù)據(jù)文件,使用Resurf插件的“從云點創(chuàng)建單張曲面”功能由點云建立邊坡地表的“曲面”模型。使用“擠出曲面”功能由邊坡的“曲面”模型建立“實體”模型,并通過“修剪”功能去除模型底部多余的部分,最后使用“將平面洞加蓋”功能即可生成如圖4a所示的封閉的邊坡三維實體模型。
a—坡表幾何模型; b—坡體幾何模型。圖4 坡體的精細化幾何模型Fig.4 Refined geometric models of slopes
當邊坡由成層土組成時,須要在模型中建立土層分界面的位置,即建立地層的幾何模型。地層幾何模型與坡表幾何模型的建立方法類似,即:在ArcMap軟件中建立地層界面柵格,再將其導入Rhino軟件中并建立地層界面的“曲面”模型。在此基礎(chǔ)上,還須通過Griddle插件將地表模型和地層模型組合起來形成一個整體,方法是:1)通過“NonManifoldMerge”命令將邊坡實體模型與地層曲面模型相組合。2)使用“testExtendSlits”命令使組合的模型符合Rhino軟件的幾何規(guī)則。3)通過“ExtractSrf”命令去除組合模型的多余部分。
加入地層后的坡體幾何模型見圖4b。
基于無人機遙感獲取的激光點云數(shù)據(jù)建立坡體精細化幾何模型后,為了對邊坡進行數(shù)值模擬分析,還須要對該幾何模型進行網(wǎng)格劃分。采用FLAC3D有限差分軟件進行邊坡的穩(wěn)定性分析,因此,須要生成FLAC3D所需的網(wǎng)格文件*.f3grid。該項工作可在Rhino軟件中實現(xiàn),步驟[10]如下:
1)使用“mesh”命令建立實體模型的表面網(wǎng)格模型。
2)通過“Gsuf”命令將步驟1)表面網(wǎng)格模型重新劃分成符合FLAC3D建模要求的表面網(wǎng)格模型。
3)使用“Gvol”命令將步驟2)表面網(wǎng)格模型轉(zhuǎn)換成符合FLAC3D網(wǎng)格規(guī)則的體網(wǎng)格模型,并將其導出為符合FLAC3D格式的網(wǎng)格文件*.f3grid。如圖5所示,模型的高約147 m,寬度約為191 m,長度約為354 m,基底厚約為41 m。該網(wǎng)格模型包含83 483個單元和67 285個結(jié)點。
圖5 Rhino中含地層信息的邊坡網(wǎng)格模型Fig.5 A grid model of slopes with stratum information in Rhino
風險(R)是指某種特定的危險事件發(fā)生的可能性(破壞概率Pf)與其產(chǎn)生的后果(C)的組合,其基本計算式為:
R=PfC
(1)
可見,在單體滑坡的定量風險分析中,兩個重要步驟是計算邊坡的破壞概率及失效后果,其共同求解基礎(chǔ)是邊坡穩(wěn)定性的三維數(shù)值模擬。
在巖土工程界通常使用安全系數(shù)Fs來定量評價邊坡的穩(wěn)定性,并認為Fs小于1時邊坡失穩(wěn)破壞。一般采用可靠度分析方法來求解邊坡的破壞概率,并取式(2)所示的極限狀態(tài)函數(shù):
Z=g(X)=Fs(X)-1
(2)
其中X=[X1,X2, …,Xn]
式中:X為隨機變量;Fs為邊坡的安全系數(shù),它是變量X的函數(shù)。
邊坡破壞概率的計算式如下:
Pf=1-Φ(β)
(3)
式中:β為邊坡的可靠指標;Φ為標準正態(tài)分布的累積函數(shù)。
研究使用常規(guī)的一階可靠度方法(FORM)求解可靠指標β,其迭代計算式見式(4):
(4)
式中:μX、σX分別為變量X的均值及均方差。
由式(2)~(4)可見,求解邊坡破壞概率的關(guān)鍵環(huán)節(jié)是求解邊坡的安全系數(shù),它可以由FLAC3D自帶的強度折減法求得。對于圖6所示的邊坡模型,在FLAC3D中對模型X方向邊界施加沿X方向的固定約束,對模型Y方向邊界施加沿Y方向的固定約束,固定模型底部邊界的位移。在此基礎(chǔ)上,即可采用solve fos命令求解邊坡的安全系數(shù)。
圖6 FLAC3D三維精細化邊坡模型Fig.6 A 3D refined slope model constructed by FLAC3D
邊坡失效后果的確定是個非常復雜的問題。為簡單起見,參考文獻[5]的方法,用滑坡體的體積V來代表邊坡災害導致的損失C,即令C=V。
文獻[11]指出:滑坡在宏觀上的空間特征就是滑體與滑床相脫離,故可以將二者的分離界面定義為邊坡的滑動面。邊坡產(chǎn)生破壞時滑坡體和滑床均存在變形,但是在滑裂破壞處的變形相對較大。在基于強度折減法的邊坡穩(wěn)定性數(shù)值模擬中,雖然可以由此繪制最大剪切應變增量等值線云來直觀顯示邊坡滑動面的位置,但難以量化表示滑坡體的體積。而從位移的角度來描述這種變形分布時,則可以認為位移等值線最為密集區(qū)域即為剪切帶,并通過臨界位移等值面來表征三維滑坡面。以圖7所示的邊坡位移等值線剖面圖為例,滑動面附近的位移等值線最密集,其位移值較小;滑體部分靠近臨空面,其位移則相對較大。因此,通過定義臨界位移值來確定滑動面位置的同時,可以依據(jù)臨界位移值區(qū)分邊坡的滑體和穩(wěn)定體,即:認為邊坡體中位移大于該臨界位移值的單元屬于邊坡的滑動破壞部分,所有該類單元體積的總和即為邊坡的滑坡體的體積。
基于這一思路,在FLAC3D中進行了二次開發(fā),具體步驟如下:
圖7 邊坡剖面的位移等值線 mFig.7 Contours of displacement for a vertical cross section of slopes
1)載入網(wǎng)格模型后,輸入“l(fā)ist zone gp”,得到所有網(wǎng)格單元編號(id)和該單元上的結(jié)點編號,將網(wǎng)格信息保存為*.txt文件以便后續(xù)調(diào)用。
2)在FLAC3D中讀取步驟1)得到的網(wǎng)格信息,并將其存入數(shù)組中。
3)遍歷步驟2)數(shù)組中每個單元,并對每個單元的所有結(jié)點位移進行判斷。該過程涉及的主要命令有g(shù)p_xdisp、gp_ydisp、gp_zdisp和find_gp。若單元上的每個結(jié)點位移均大于臨界位移值,則將該單元的體積保存至一個數(shù)組,該過程涉及的主要命令是z_volume。
4)對存儲體積的數(shù)組求和,即為滑體體積。
該方法的一個重要步驟是確定邊坡的臨界位移值。參考文獻[11]的建議,根據(jù)特定剖面邊界的剪應變增量極大值點確定滑體剪出口位置,據(jù)此選擇邊坡的臨界位移值為0.2 m。以邊坡的定值法分析為例,求得的對應的滑體體積為9.683×105m3,約占邊坡總體積(7.465×106m3)的13.0%。
基于FLAC3D進行二次開發(fā)求解滑坡體積方法的優(yōu)點是能方便地求解滑動面以上、坡表以下所有單元的體積。當坡表的幾何形態(tài)由無人機遙感精確獲取后,即可準確地確定任意形狀滑坡體的體積,在此基礎(chǔ)上可以方便地進行滑坡災害的風險分析。
對單體滑坡災害進行概率分析時,須要明確研究對象穩(wěn)定性的影響因素及該因素所服從的概率分布。圖6所示的邊坡由碎石土層及基巖組成,影響邊坡穩(wěn)定性的主要變量是上層土的強度參數(shù)值,記上層土的黏聚力及內(nèi)摩擦角分別為c及φ,則基本變量X=[c,φ]。設(shè)基本變量服從正態(tài)分布,變異系數(shù)為0.1~0.4。在此條件下,邊坡的破壞概率、失效后果及風險值與強度參數(shù)變異系數(shù)的關(guān)系分別見圖8~10。圖中,顏色由藍至紅表示表示數(shù)值逐漸增加。
由圖8可知:隨著黏聚力和內(nèi)摩擦角變異系數(shù)的增大,邊坡的破壞概率逐漸增加。其中,黏聚力變異系數(shù)對邊坡破壞概率的影響較小,而內(nèi)摩擦角變異系數(shù)對邊坡破壞概率影響較大。內(nèi)摩擦角變異系數(shù)的增加會導致邊坡破壞概率的大幅增加;內(nèi)摩擦角變異系數(shù)較大時,黏聚力變異系數(shù)的變化對邊坡破壞概率的影響程度減弱。
圖8 邊坡破壞概率與強度參數(shù)變異系數(shù)的關(guān)系Fig.8 Relations between failure probability of slopes and variation coefficients of strength parameters
表3給出了邊坡的性能等級與破壞概率的關(guān)系。[12]結(jié)合圖8及表3可知:當內(nèi)摩擦角變異系數(shù)為0.4時,邊坡達到了危險級別,此時須對該邊坡開展進一步監(jiān)測。在黏聚力變異系數(shù)為0.1而內(nèi)摩擦角變異系數(shù)從0.1增長到0.4的過程中,邊坡的性能從平均以上級別降至危險級別,說明在工程實踐中要重視對邊坡內(nèi)摩擦角變異系數(shù)的取值。
表3 邊坡的預期性能等級[12]Table 3 Expected performance grades of slopes
由圖9可知:邊坡失效后果近似隨著黏聚力變異系數(shù)的減小及內(nèi)摩擦角變異系數(shù)的增大而增大,直至趨于穩(wěn)定。此外,黏聚力變異性較小時,會一定程度上減少內(nèi)摩擦角變異性對失效后果的影響。而對于上述危險級別的邊坡(內(nèi)摩擦角變異系數(shù)為0.4),其發(fā)生破壞時還可能會造成更加嚴重的后果。
由圖10可知:隨著黏聚力和內(nèi)摩擦角變異系數(shù)的增大,滑坡災害風險均隨之增加。圖10總體上與圖8較為類似,說明對于該例而言,滑坡災害風險主要由破壞概率控制。
圖9 邊坡失效后果與強度參數(shù)變異系數(shù)的關(guān)系Fig.9 Relations between failure consequences of slopes and variation coefficients of strength parameters
上述算例分析表明,將無人機遙感與邊坡穩(wěn)定的數(shù)值模擬及可靠度分析方法相結(jié)合,可以定量地確定滑坡災害的風險。若研究對象有多個邊坡,還可以基于滑坡災害的量化風險分析結(jié)果確定各個邊坡的危險等級,并對所有邊坡進行重要性排序。
圖10 滑坡災害風險與強度參數(shù)變異系數(shù)的關(guān)系 104Fig.10 Relations between landslide risks and variation coefficients of strength parameters
提出了一種基于無人機遙感的滑坡災害風險分析方法,該方法將無人機遙感、邊坡穩(wěn)定性的三維數(shù)值模擬及可靠度分析方法相結(jié)合,具有低成本、靈活快速且結(jié)果較為準確的優(yōu)點,為線路設(shè)施巡檢和滑坡災害風險的快速響應提供了新的思路。主要結(jié)論如下:
1)基于無人機遙感獲取的坡表激光點云數(shù)據(jù),結(jié)合ArcMap及Rhino軟件可以建立坡體精細化的幾何模型及數(shù)值計算模型。
2)在FLAC3D中進行了二次開發(fā),提出了根據(jù)臨界位移值識別滑坡體并計算滑坡體體積的方法。
3)量化分析了強度參數(shù)變異系數(shù)與滑坡災害風險的關(guān)系,得到了強度參數(shù)變異系數(shù)與滑坡災害風險的關(guān)系云。