基于無人機和Rockfall Analyst的崩塌落石特征分析與運動學(xué)模擬——以察雅縣崩塌落石為例

劉福臻，李旭德，王軍朝，劉建康，陳 龍，張佳佳

(1.西南石油大學(xué) 土木工程與測繪學(xué)院，四川 成都 610500; 2.中國地質(zhì)調(diào)查局地質(zhì)災(zāi)害防治技術(shù)中心，四川 成都 611734)

崩塌落石作為常見的地質(zhì)災(zāi)害，是指斜坡上的部分危巖體受到力的作用后發(fā)生滾落的現(xiàn)象，是邊坡勘察防護中的重點[1]。引發(fā)崩塌落石的條件簡單，容易崩落，截至2018年昌都地區(qū)就共有2 192處地質(zhì)災(zāi)害及隱患點，崩塌共計647處，占總數(shù)的29.52%，其中小型崩塌落石共計403處[2]。目前針對崩塌落石災(zāi)害，在形成機制、運動學(xué)模擬、易發(fā)性等方面的研究成果較多，調(diào)查手段也越來越豐富。

董秀軍等[3]在2006年就提出將三維激光掃描技術(shù)應(yīng)用于高陡邊坡地質(zhì)調(diào)查測量；王鳳艷等[4]也早在2008年將近景攝影測量技術(shù)用在巖體結(jié)構(gòu)面調(diào)查。但是二者對調(diào)查環(huán)境的要求較高，復(fù)雜環(huán)境下工作難以開展。隨著小型無人機技術(shù)的快速發(fā)展，憑借其高精度、靈活機動、實用性強，能克服危險環(huán)境影響，可提供多角度、高分辨率影像等獨特優(yōu)點，近幾年已逐漸應(yīng)用于各種災(zāi)害地質(zhì)調(diào)查的研究[5-13]。

影響滾石運動軌跡和規(guī)律的因素多且情況復(fù)雜[14]。在國內(nèi)已有眾多學(xué)者基于數(shù)值模擬軟件研究影響滾石運動的參數(shù)和運動規(guī)律，常用的是二維模擬軟件Rockfall和三維模擬軟件Rockyfor3D。Rockfall是一個基于統(tǒng)計的軟件程序，計算時將塊石的崩落過程當(dāng)作一個質(zhì)點，在人為規(guī)定好的剖面線上讓其墜落、碰撞、滾落，應(yīng)用極為廣泛[15-18]。但其崩塌路徑是斜坡上的特定剖面，具有一定的人為控制性。Rockyfor3D作為三維模擬軟件，在森林地區(qū)和鐵路沿線的落石災(zāi)害評估上都有應(yīng)用[19-22]。該軟件在模擬計算時會考慮單個塊石的形狀和尺寸大小，通常塊石形狀精確數(shù)據(jù)的獲取具有一定的難度[23]。

Rockfall Analyst是基于地理信息系統(tǒng)，將分布式建模和物理落石過程相結(jié)合的三維模擬軟件。模擬時采用集中質(zhì)量法，使用矢量和柵格數(shù)據(jù)，能在ArcGIS平臺有效呈現(xiàn)落石軌跡、飛行彈跳高度以及速度等的分布，在山區(qū)落石災(zāi)害預(yù)測以及風(fēng)險評價方面都有應(yīng)用[24-25]。相比Rockyfor3D和Rockfall，該軟件不需要統(tǒng)計落石大小、形狀的參數(shù)，也無需人為規(guī)定落石運動的剖面方向，整體考慮的巖體參數(shù)較少，又減少了人為干擾因素，因此更適用于本文研究區(qū)環(huán)境下崩塌落石的運動學(xué)模擬。

察雅縣為G349國道起始路段，沿線多為高山峽谷地段，每逢雨季沿線崩塌落石災(zāi)害頻發(fā)，而且隨機性大、物源位置分布高、環(huán)境復(fù)雜，傳統(tǒng)的人工手持羅盤等工具的接觸式調(diào)查無法完成，由于公路上車流量較多，采用近景攝影測量和三維激光掃描技術(shù)方法不僅極為不便，而且也無法全面調(diào)查評估災(zāi)害特征，不能滿足調(diào)查需求。復(fù)雜的地形環(huán)境為災(zāi)害調(diào)查評估和防治帶來了較大的困難，嚴(yán)重威脅著過往車輛安全，調(diào)查分析研究區(qū)崩塌落石災(zāi)害特征及威脅范圍，可為該處落石災(zāi)害防治以及沿線災(zāi)害評估提供技術(shù)借鑒。因此，本文借助現(xiàn)場調(diào)查和無人機三維傾斜攝影測量技術(shù)，分析了研究區(qū)崩塌落石物源區(qū)巖體特征、堆積區(qū)特征以及巖體穩(wěn)定性。獲取了研究區(qū)高精度DEM，結(jié)合基于ArcGIS的三維模擬軟件Rockfall Analyst，對物源區(qū)塊石進行了運動學(xué)模擬。

1 研究區(qū)概況

昌都市位于青藏高原東部，地處橫斷山脈以東高山深谷區(qū)，海拔3 000～6 378 m，平均海拔3 500 m，地勢北高南低。在大地構(gòu)造上， 區(qū)域上處于巨型青藏滇緬構(gòu)造， 淺層地殼結(jié)構(gòu)特征為板片拼合結(jié)構(gòu)[2]。區(qū)域基本地震加速度為0.2g，地震烈度為VII度，50年內(nèi)超越該烈度值的可能性為10%[26]。察雅縣位于昌都市南東側(cè)，縣城屬于寬谷向峽谷過度的地貌，主河麥曲從縣城穿過。縣城南側(cè)斜坡基巖出露少，整體為滑坡群，縣域內(nèi)地層以侏羅系汪布組和三疊系奪蓋拉組為主(圖1)。

圖1 崩塌點周邊地質(zhì)圖Fig.1 Geological map of the surrounding of the rockfall

本文研究區(qū)位于縣城北西側(cè)的一處巖質(zhì)斜坡，坡頂高程3 187 m，坡腳處高程3 094 m。該斜坡基巖以紫紅色粉砂質(zhì)泥巖為主，夾少量淺灰綠色粉砂巖，差異性風(fēng)化嚴(yán)重。崩塌落石點所在坡體巖體破碎，結(jié)構(gòu)面發(fā)育，曾發(fā)生過多次小型崩塌落石，坡腳形成堆積區(qū)，目前仍存在多處危險源，物源分布位置高，地勢陡峭，最低一處物源與坡腳公路的高差為15 m，傳統(tǒng)的人工調(diào)查方式、近景攝影測量和三維激光掃描技術(shù)方法均無法滿足調(diào)查需求。高位物源塊石一旦誘發(fā)崩落，對坡腳G349國道、鄉(xiāng)道和果園產(chǎn)生著極大的威脅，如圖2所示。

圖2 崩塌落石位置概況Fig.2 Location of collapse and rockfall

2 研究方法

2.1 無人機三維傾斜攝影測量

相比傳統(tǒng)的攝影測量，無人機傾斜攝影可以從多個不同角度拍攝采集數(shù)據(jù)，作業(yè)機動靈活，受地形環(huán)境影響小，能夠完成傳統(tǒng)的人工手段無法進行的高陡邊坡危巖體調(diào)查工作，構(gòu)建的三維模型可以真實反映地物特征，實現(xiàn)非接觸式測量。無人機傾斜攝影測量的災(zāi)害調(diào)查流程主要由調(diào)查區(qū)資料準(zhǔn)備、野外拍攝和內(nèi)業(yè)影像數(shù)據(jù)處理和分析結(jié)果現(xiàn)場驗證四部分組成。流程圖如圖3所示。調(diào)查區(qū)資料包括遙感影像、地形資料、歷史災(zāi)害調(diào)查資料等，野外拍攝包括布像控點和飛行拍攝，內(nèi)業(yè)影像處理主要為影像三維重建。

圖3 調(diào)查流程Fig.3 Investigation process

2.2 Rockfall Analyst運動學(xué)模擬

災(zāi)害模擬分析和評估災(zāi)害危險性對于災(zāi)害防治具有實在意義，常用的分析評估方法為物理力學(xué)模型分析法[27]。力學(xué)數(shù)值模擬模型在牛頓運動定律和碰撞理論的基礎(chǔ)上，通過大量試驗分析總結(jié)出影響落石的特征參數(shù)，利用運動學(xué)和動力學(xué)完成滾石運動軌跡、速度等的模擬[28,29]。本文采用基于ArcGIS的三維數(shù)值模擬軟件Rockfall Analyst來模擬研究區(qū)崩塌落石的動力學(xué)。該模型采用了集中質(zhì)量法，落石動力學(xué)不受塊石的質(zhì)量和形狀的影響，分析時只考慮震源位置和初始方向的內(nèi)在隨機性。巖石和邊坡之間的相互作用由兩個恢復(fù)系數(shù)(RN、RT)和一個摩擦系數(shù)特定描述，在模擬過程中，這3個系數(shù)是不變的[30]。

該軟件可以模擬落石飛行、碰撞反彈和滾動這3個主要的運動過程。落石運動的路徑、高度和速度通過拋物線方程來模擬。位移和速度如式(1)(2)所示：

(1)

(2)

其中，g為重力加速度(取9.8m/s2)，X0，Y0，Z0是滾石在三維空間中的初始位置，Vx0，Vy0，Vz0是滾石在x、y和z方向的初始速度。

碰撞和反彈階段首先是要確定巖石飛行結(jié)束時的撞擊點，即飛行滾石和柵格表面的交點，然后由質(zhì)點碰撞反彈模型(圖4)表示彈跳速度矢量。彈跳速度矢量根據(jù)恢復(fù)系數(shù)(RN、RT)計算。RN、RT可以在野外現(xiàn)場用高速攝影機搭配解析軟件測得[31]，數(shù)學(xué)表達式如式(3)，彈跳速度矢量由式(4)定義：

(3)

(4)

式中，VDip，VTrend，VN為完整的彈跳速度矢量，VDip為傾角方向上的速度分量，VTrend為法向方向上的速度分量，RN、RT為恢復(fù)系數(shù)，范圍在[0，1]。

滾石沿坡面運動過程中，根據(jù)牛頓運動定律對加速和速度變換等模擬計算，其受力分析(圖5)和計算公式如下。

m石a=Gsinθ-F阻.

(5)

V末=V初-at.

(6)

3 研究結(jié)果

3.1 崩塌落石特征

(1)物源區(qū)巖體特征

根據(jù)建模生成的報告可知模型平面精度3 cm，高程精度6 cm，滿足災(zāi)害調(diào)查精度要求。從三維模型觀察到研究區(qū)共有3處物源區(qū)，從低到高依次劃分為A區(qū)、B區(qū)和C區(qū)，如圖6(a)所示。根據(jù)模型統(tǒng)計出A區(qū)物源方量為710 m3，B區(qū)為1 530 m3，C區(qū)為450 m3。3處巖體整體受3組主要結(jié)構(gòu)面J1、J2、J3控制，被切割成塊狀(圖6(b))。其中A區(qū)塊度最大，出露巖面也最新，B區(qū)次之，C區(qū)最小。物源B區(qū)巖體表面X形節(jié)理極其發(fā)育(圖6(c))，可見3處裂縫，由三維模型統(tǒng)計了長度，其中受J3節(jié)理控制的兩處裂縫分別長1.5 m和0.75 m，J2控制的一條裂縫長2.5 m(圖6(d))。

圖6 巖體結(jié)構(gòu)面Fig.6 Structural plane of rock mass

除了3組主控結(jié)構(gòu)面，巖體還發(fā)育部分次生結(jié)構(gòu)面。由于結(jié)構(gòu)面強度低于巖塊強度，不同結(jié)構(gòu)面之間的組合形式產(chǎn)生了多種破壞方式，而且結(jié)構(gòu)面容易積水，高寒地帶冰雪凍融交替作用使得結(jié)構(gòu)面間隙擴張，導(dǎo)致進一步加重了崩塌落石的發(fā)生。

軟、硬巖風(fēng)化速度不同而形成空腔導(dǎo)致巖體失穩(wěn)是一種常見的巖體破壞特征。物源區(qū)巖體巖性整體為粉砂質(zhì)泥巖夾粉砂巖，呈“軟-硬-軟”夾層狀。由于粉砂巖的抗風(fēng)化能力比粉砂質(zhì)泥巖要強，在受到暴雨沖刷和高原冰雪凍融侵蝕破壞后，粉砂質(zhì)泥巖風(fēng)化剝蝕速度加快，從而形成空腔(圖7)，粉砂巖失去支撐力，同時又被結(jié)構(gòu)面切割成不穩(wěn)定塊體，在自身重力和外力作用到一定程度時，危巖塊體失穩(wěn)崩落。

(2)巖體結(jié)構(gòu)面特征

巖體結(jié)構(gòu)面是巖體的主要特征，對巖體穩(wěn)定性起控制作用。巖體結(jié)構(gòu)面在三維點云數(shù)據(jù)中一般都表現(xiàn)為規(guī)則的平面，可以通過提取平面上不在同一條直線上的構(gòu)成范圍較大，且內(nèi)角均勻的任意n個特征點的三維坐標(biāo)(n≥3)來計算，傾向和傾角誤差可以控制在5°以內(nèi)[32-33]。

為了驗證由三維模型獲取的巖體產(chǎn)狀的精度，從坡腳選取5處平整的面(圖8)并提取每個面的4組三維坐標(biāo)，基于Matlab語言編寫計算程序算出每個結(jié)構(gòu)面的3組產(chǎn)狀，求出平均值，然后用手工羅盤實測每個面的產(chǎn)狀，將2種方式算得的產(chǎn)狀值進行對比，如表1所示。結(jié)果發(fā)現(xiàn)傾向最大誤差11°，最小誤差7°，平均誤差8.6°，傾角誤差最大為6°，最小0°，平均誤差2.2°，均符合地質(zhì)災(zāi)害調(diào)查要求。

表1 巖體產(chǎn)狀驗證Table 1 Verification of rock mass occurrence

驗證了由三維模型獲取的巖體產(chǎn)狀精度符合調(diào)查要求后，通過提取J1、J2、J3每組結(jié)構(gòu)面的四組三維坐標(biāo)，同樣基于Matlab語言編寫的計算程序算出每個結(jié)構(gòu)面的3組產(chǎn)狀，求出平均值，并統(tǒng)計結(jié)構(gòu)面間距，如表2所示。

表2 物源區(qū)巖體產(chǎn)狀Table 2 Occurrence of rock mass in provenance area

(3)落石堆積特征

調(diào)查發(fā)現(xiàn)堆積區(qū)塊石都是零散崩落堆積在坡腳。統(tǒng)計了粒徑大于1 m的塊石共有15處，最大粒徑2.3 m。堆積最遠(yuǎn)的一處塊石粒徑2.2 m，距離物源A區(qū)斜距51 m，距G349國道3.5 m，坡腳處未做任何防治措施，如圖9所示。從三維模型統(tǒng)計的結(jié)構(gòu)面間距來看，3組結(jié)構(gòu)面切割形成的塊石中，粒徑大于1 m的均在物源A區(qū)。3處物源區(qū)中A區(qū)有明顯崩落的痕跡，參考上述圖6(b)，并且現(xiàn)場可以觀察到A區(qū)出露的巖面較新，無明顯風(fēng)化現(xiàn)象，再結(jié)合塊石堆積的位置，綜合判斷坡腳堆積的大塊石均來自于物源A區(qū)。確定現(xiàn)有的堆積大塊石來自于A區(qū)，若B區(qū)和C區(qū)的塊石失穩(wěn)崩落，能威脅多大范圍，需要進行數(shù)值模擬分析判斷。

圖9 崩塌堆積區(qū)Fig.9 Partial accumulation area of rockfalls

3.2 巖體穩(wěn)定性分析

將J1、J2、J3的產(chǎn)狀和斜坡坡向、坡度(J4)繪制到極射赤平投影圖上，如圖10所示。其中J1和J2組合交線的交點位于J4的外側(cè)，因此交線傾角小于坡角，且交線傾向與坡向相同，所以切割體不穩(wěn)定；J1、J3的組合交線傾向和J2、J3組合交線傾向方向與坡向相反，切割體相對穩(wěn)定。影響巖體穩(wěn)定性的除了主控結(jié)構(gòu)面之外，部分次生結(jié)構(gòu)面也可導(dǎo)致切割體失穩(wěn)?；贖OEK楔形體模型計算出危巖穩(wěn)定系數(shù)，在降雨充水狀況下為0.99，干燥狀況下為2.2，由此說明在干燥狀況下穩(wěn)定，在降雨充水狀況下不穩(wěn)定，穩(wěn)定系數(shù)參考依據(jù)如表3所示。危巖塊石失穩(wěn)后的運動路徑復(fù)雜且難以確定，需要用數(shù)值模擬軟件進行分析預(yù)測。

圖10 赤平投影圖Fig.10 Stereographic projection

表3 穩(wěn)定系數(shù)Table 3 Stability coefficient

3.3 落石數(shù)值模擬

(1)參數(shù)獲取

利用Rockfall Analyst對落石模擬分析時，DEM精度、恢復(fù)系數(shù)(RN、RT)和內(nèi)摩擦角(Deg)的取值以及落石源點的設(shè)置對模擬結(jié)果的準(zhǔn)確性影響較大。高精度DEM由無人機成果獲取，落石源點根據(jù)高精度三維模型和正射影像結(jié)合現(xiàn)場調(diào)查確定。確定研究區(qū)恢復(fù)系數(shù)和內(nèi)摩擦角時根據(jù)歷史崩塌落石的堆積位置，采用參數(shù)反演的方法確定：先根據(jù)參數(shù)取值范圍確定大概數(shù)值，然后進行試算模擬，最后將模擬結(jié)果與現(xiàn)場調(diào)查的落石堆積位置作對比，最終確定該區(qū)域斜坡下墊面的恢復(fù)系數(shù)和內(nèi)摩擦角[19]。

通過調(diào)整參數(shù)，模擬A區(qū)落石生成的運動軌跡如圖11所示，模擬的運動軌跡最遠(yuǎn)到達坡腳，沒有危及到G349國道。對比發(fā)現(xiàn)模擬的軌跡終點與現(xiàn)場調(diào)查的塊石堆積體分布情況基本吻合(圖12)，由此確定生成該軌跡的下墊面參數(shù)作為其它物源區(qū)塊石的模擬參數(shù)，如表4所示。在對A區(qū)落石模擬時，除了對已崩落的塊石進行反演模擬，還模擬了大量潛在崩落的塊石，生成的軌跡沒有對G349國道造成破壞，因此不再對A區(qū)潛在崩落塊石的速度和彈跳高度進行模擬分析。

表4 下墊面參數(shù)Table 4 Parameters of underlying surface

(2)運動特征分析

根據(jù)模擬得出的參數(shù)對物源B區(qū)和C區(qū)塊石進行模擬，最終得到落石的運動軌跡如圖13(a)、(b)和圖14所示，彈跳高度的柵格圖如圖13(c)、(d)所示。

青島工學(xué)院工程管理類本科的教學(xué)模式同樣存在重理論輕實踐、授課方式單一、課程內(nèi)容狹窄等問題，本研究主要從以下三個方面探討工程管理專業(yè)教學(xué)模式創(chuàng)新[5]：

圖13 B、C區(qū)落石和彈跳高度軌跡Fig.13 Trajectory of rockfall and bounce height in areas B and C

從圖13(a)、(b)可以看出，兩處物源區(qū)落石都危及到了G349國道和鄉(xiāng)道。其中B區(qū)西側(cè)落石最遠(yuǎn)堆積到了國道外側(cè)，東側(cè)落石大部分堆積到鄉(xiāng)道，沒有危及到果園；C區(qū)西側(cè)落石基本都堆積到斜坡中部，對國道沒有大的威脅，東側(cè)部分落石不僅危及省道和果園，最遠(yuǎn)堆積到國道外側(cè)公共衛(wèi)生場所附近。彈跳高度柵格圖如圖13(c)、(d)所示，B區(qū)落石運動過程中最大彈跳高度為7 m，C區(qū)落石最大彈跳高度為6 m。結(jié)合落石的運動三維軌跡圖(圖14)可以看出，B區(qū)堆積到國道外側(cè)的落石都是從坡腳處以弧線運動軌跡飛躍過去，C區(qū)部分落石直接從斜坡中部彈起飛躍鄉(xiāng)道進入果園。結(jié)合現(xiàn)場調(diào)查，鄉(xiāng)道一側(cè)因修路切坡形成一高4 m的陡坎，因此模擬的彈跳高度符合實際。

圖14 落石三維軌跡Fig.14 Three dimensional trajectory of rockfall

兩處物源區(qū)落石的運動速度柵格圖如圖15、16所示。B區(qū)落石速度最大為9 m/s，C區(qū)落石速度最大為15 m/s。結(jié)合野外調(diào)查將B區(qū)塊石的平均質(zhì)量設(shè)為800 kg，C區(qū)塊石的平均質(zhì)量設(shè)為200 kg，算得B區(qū)落石的最大沖擊動能為32.4 KJ，C區(qū)落石的最大沖擊動能為22.5KJ，如表5所示。

表5 模擬值Table 5 Analog value

3.4 防治建議

結(jié)合現(xiàn)場調(diào)查和運動學(xué)模擬來看，B、C兩區(qū)落石對公路和果園產(chǎn)生的威脅大，應(yīng)及時做防治措施。在靠近G349國道這一側(cè)，大部分落石到坡腳處速度為0，如果修筑擋墻，受到落石的沖擊力最小，但是無法避免部分飛躍的落石，仍然會危及國道；如果在落石飛躍之前做防治，雖受到落石沖擊力最大，但防護效果最好。在鄉(xiāng)道這一側(cè)，如果在坡腳處修擋墻能防止B區(qū)的落石危及鄉(xiāng)道，但是C區(qū)的落石將直接飛入果園，因此建議修筑被動網(wǎng)才能做到最有效的防治。

4 結(jié)論與討論

根據(jù)現(xiàn)場調(diào)查結(jié)合無人機傾斜攝影測量技術(shù)，實現(xiàn)了對昌都市察雅縣一處崩塌落石的精細(xì)化調(diào)查，并利用Rockfall Analyst軟件完成了危巖塊石威脅范圍的數(shù)值模擬。主要結(jié)論如下。

(1)分析了物源區(qū)結(jié)構(gòu)面的基本信息，得到三組結(jié)構(gòu)面產(chǎn)狀：J1：217°∠71°，J2：82°∠38°，J3：315°∠63°，并統(tǒng)計了結(jié)構(gòu)面間距。物源區(qū)巖體變形破壞特征主要表現(xiàn)在不同結(jié)構(gòu)面組合形式上的破壞和地層巖性上的差異性分化形成的破壞。

(2)統(tǒng)計了堆積區(qū)粒徑大于1 m的塊石，為數(shù)值模擬過程中的參數(shù)反演確定了塊石來源。通過參數(shù)反演確定了研究區(qū)下墊面的參數(shù)，反演結(jié)果和現(xiàn)場進行調(diào)查分析的結(jié)果基本吻合，保證了后期模擬結(jié)果的準(zhǔn)確性。

(3)從三維模擬結(jié)果來看，物源B區(qū)和C區(qū)的落石均能危及到G349國道、鄉(xiāng)道和果園。B區(qū)落石距地面的最大彈跳高度為7 m，C區(qū)為6 m；B區(qū)落石的最大動能為32.4 KJ，C區(qū)為22.5 KJ。通過模擬確定了滾石能夠威脅的范圍、彈跳高度以及運動速度，有效評估了災(zāi)害影響范圍。將運動軌跡以三維的形式展示，不僅消除了二維數(shù)值模擬軟件人為控制因素的影響，而且能直觀反映落石運動過程中彈跳和飛行的空間變化，有利于防護位置的選擇。通過計算落石威脅到公路時的最大沖擊動能，在防護強度的選擇上提供了數(shù)據(jù)依據(jù)。

致謝：感謝中國地質(zhì)科學(xué)院探礦工藝研究所提供的實習(xí)機會，以及各位老師、師兄給與的指導(dǎo)與幫助！感謝蘭恒星老師團隊無償提供的Rockfall Analyst軟件，感謝仉義星師兄在軟件使用過程中的無私指導(dǎo)！