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

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

(1.西南石油大學 土木工程與測繪學院，四川 成都 610500; 2.中國地質調查局地質災害防治技術中心，四川 成都 611734)

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

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

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

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

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

1 研究區概況

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

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

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

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

2 研究方法

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

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

圖3 調查流程Fig.3 Investigation process

2.2 Rockfall Analyst運動學模擬

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

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

(1)

(2)

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

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

(3)

(4)

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

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

m石a=Gsinθ-F阻.

(5)

V末=V初-at.

(6)

3 研究結果

3.1 崩塌落石特征

(1)物源區巖體特征

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

圖6 巖體結構面Fig.6 Structural plane of rock mass

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

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

(2)巖體結構面特征

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

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

表1 巖體產狀驗證Table 1 Verification of rock mass occurrence

驗證了由三維模型獲取的巖體產狀精度符合調查要求后，通過提取J1、J2、J3每組結構面的四組三維坐標，同樣基于Matlab語言編寫的計算程序算出每個結構面的3組產狀，求出平均值，并統計結構面間距，如表2所示。

表2 物源區巖體產狀Table 2 Occurrence of rock mass in provenance area

(3)落石堆積特征

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

圖9 崩塌堆積區Fig.9 Partial accumulation area of rockfalls

3.2 巖體穩定性分析

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

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

表3 穩定系數Table 3 Stability coefficient

3.3 落石數值模擬

(1)參數獲取

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

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

表4 下墊面參數Table 4 Parameters of underlying surface

(2)運動特征分析

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

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

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

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

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

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

表5 模擬值Table 5 Analog value

3.4 防治建議

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

4 結論與討論

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

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

(2)統計了堆積區粒徑大于1 m的塊石，為數值模擬過程中的參數反演確定了塊石來源。通過參數反演確定了研究區下墊面的參數，反演結果和現場進行調查分析的結果基本吻合，保證了后期模擬結果的準確性。

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

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