劉寧,袁維*,高嶺,白雪山,翟星,李建朋,王安禮,徐江

(1.石家莊鐵道大學(xué)土木工程學(xué)院，石家莊 050043；2.河北省交通規(guī)劃設(shè)計(jì)院有限公司，石家莊 050300；3.河北省地質(zhì)環(huán)境監(jiān)測(cè)院，石家莊 050022；4.河北省地質(zhì)資源環(huán)境監(jiān)測(cè)與保護(hù)重點(diǎn)實(shí)驗(yàn)室，石家莊 050022；5.貴州省喀斯特地區(qū)水資源開發(fā)利用工程技術(shù)研究中心，貴陽(yáng) 550002)

在大型工程建設(shè)中，無(wú)論是公路鐵路還是水利水電都難免會(huì)跨越或者下穿巖質(zhì)邊坡，所以需要保證巖質(zhì)邊坡具有良好的穩(wěn)定性。巖質(zhì)邊坡穩(wěn)定性分析一般需要首先進(jìn)行野外勘察收集數(shù)據(jù)，然后進(jìn)行相關(guān)試驗(yàn)分析獲取力學(xué)參數(shù)，最后運(yùn)用數(shù)值模擬等方法對(duì)巖質(zhì)邊坡進(jìn)行計(jì)算評(píng)估。野外勘察收集數(shù)據(jù)通常采用人工測(cè)量和原地取試樣的方法[1]，該方法會(huì)受到復(fù)雜的地形地貌條件制約，使正常的勘察調(diào)查工作無(wú)法開展，并會(huì)威脅測(cè)量人員的安全，而試驗(yàn)環(huán)節(jié)獲取參數(shù)的過程漫長(zhǎng)，且存在人為操作的誤差，影響試驗(yàn)結(jié)果，因此，需要找到一種快捷、高效的獲取地質(zhì)數(shù)據(jù)并得到相關(guān)參數(shù)的方法。

國(guó)內(nèi)外學(xué)者對(duì)上述問題進(jìn)行了廣泛的研究，無(wú)論是地形地貌的測(cè)量方面，還是根據(jù)地質(zhì)數(shù)據(jù)進(jìn)行相關(guān)參數(shù)的提取方面都取得了不少成果。梁玉飛等[2]把激光三維掃描技術(shù)應(yīng)用到地質(zhì)工程中，通過激光三維掃描獲取研究區(qū)域滑坡的點(diǎn)云數(shù)據(jù)，以此為基礎(chǔ)進(jìn)行節(jié)理信息的提取，為地質(zhì)災(zāi)害分析提供了數(shù)據(jù)支撐。嵇美偉等[3]利用近景攝影測(cè)量技術(shù)提取了巖石節(jié)理信息，結(jié)合Python編寫的兩種算法實(shí)現(xiàn)了用非測(cè)量相機(jī)完成高精度信息的提取。激光三維掃描技術(shù)和近景攝影測(cè)量技術(shù)近些年被應(yīng)用在各種工程，同時(shí)其缺點(diǎn)也暴露了出來(lái)，如有障礙物遮擋或者在險(xiǎn)峻的地形下，這兩種方法實(shí)現(xiàn)的過程就變得復(fù)雜。隨著小型無(wú)人機(jī)的發(fā)展，其以快捷、高效地獲取數(shù)據(jù)等優(yōu)點(diǎn)被廣泛應(yīng)用于地質(zhì)調(diào)查中。Goncalves等[4]采用無(wú)人機(jī)攝影測(cè)量方法，通過設(shè)置地面控制點(diǎn)，進(jìn)行繪制和監(jiān)測(cè)了沿海地區(qū)沙丘和海灘的情況，再結(jié)合三維重建軟件對(duì)沙丘和海灘的地形進(jìn)行模型重構(gòu)，結(jié)果表明無(wú)人機(jī)可獲取高精度的圖像，構(gòu)建的三維模型也符合實(shí)際。周成凱等[5]利用小型無(wú)人機(jī)對(duì)高位危巖進(jìn)行了數(shù)據(jù)采集，通過采用赤平投影和圖像識(shí)別相結(jié)合的方法對(duì)危巖穩(wěn)定性進(jìn)行評(píng)價(jià)，為地質(zhì)信息獲取提供了一種新方法。連會(huì)青等[6]利用無(wú)人機(jī)傾斜攝影測(cè)量技術(shù)獲取地質(zhì)數(shù)據(jù)，然后利用Context Capture軟件進(jìn)行三維重構(gòu)，再把三維模型導(dǎo)入軟件Acute 3D viewer中實(shí)現(xiàn)了對(duì)地質(zhì)信息的識(shí)別提取，結(jié)果表明得到的數(shù)據(jù)和野外地質(zhì)調(diào)查的數(shù)據(jù)一致性較好，具有實(shí)用意義。戴均豪等[7]利用小型無(wú)人機(jī)對(duì)甘肅北山地區(qū)進(jìn)行數(shù)據(jù)采集，再利用Photoscan軟件進(jìn)行三維模型的合成，對(duì)三維模型進(jìn)行解譯標(biāo)志和地質(zhì)解譯，得到了精細(xì)的地質(zhì)圖。

可見，小型無(wú)人機(jī)在地質(zhì)調(diào)查中扮演的角色越來(lái)越重要了，同時(shí)也暴露出現(xiàn)階段的研究存在著不足：①?gòu)暮綔y(cè)照片到三維模型的重建，再到地質(zhì)信息的識(shí)別與提取，其過程要經(jīng)過多個(gè)軟件的導(dǎo)入與導(dǎo)出，十分麻煩；②上述的三維模型重構(gòu)軟件多為商用軟件，精度達(dá)不到巖土工程的要求；③地質(zhì)信息識(shí)別與提取的軟件較少，且不一定是符合工程需求。這些不足之處也是目前重要的研究方向，學(xué)者們做了很多嘗試。楊帥等[8]引入多視圖立體視覺算法，通過對(duì)無(wú)人機(jī)位姿數(shù)據(jù)的研究提出了一種無(wú)人機(jī)影像匹配的新算法，結(jié)果表明新的方法無(wú)論在查準(zhǔn)率還是查全率上都有明顯的替身。王楠[9]引入了區(qū)域生長(zhǎng)法用于特征點(diǎn)的匹配，并對(duì)用于三維模型重構(gòu)的隨機(jī)抽樣一致性算法進(jìn)行了改進(jìn)，減小了三維重構(gòu)過程中的誤差，使得點(diǎn)云集的匹配結(jié)果更加準(zhǔn)確。Riquelme等[10]提出了基于MATLAB的結(jié)構(gòu)面識(shí)別算法，采用了領(lǐng)域?qū)c(diǎn)的方法，搜索相同特性的點(diǎn)分為一組并對(duì)其進(jìn)行去噪處理，可以高效找到巖體結(jié)構(gòu)面并計(jì)算其產(chǎn)狀。余宏明等[11]利用無(wú)人機(jī)攝影測(cè)量的方法，獲取了測(cè)區(qū)高精度的航拍影像，通過空三計(jì)算得到了測(cè)區(qū)地形地貌的三維模型，以此為基礎(chǔ)結(jié)合GeoStudio數(shù)值模擬軟件，對(duì)其滲流場(chǎng)和穩(wěn)定性進(jìn)行分析，為計(jì)算降雨條件下滑坡的穩(wěn)定性提供了一種快捷的方法。宋珺敏[12]以圖像數(shù)據(jù)為基礎(chǔ)，結(jié)合點(diǎn)云空間數(shù)據(jù)，得到以平面形態(tài)出露的結(jié)構(gòu)面，再利用點(diǎn)云聚類的方法得到了巖體結(jié)構(gòu)面節(jié)理信息，該方法的到的節(jié)理信息較符合實(shí)際。Liu等[13]利用無(wú)人機(jī)攝影技術(shù)采集邊坡數(shù)據(jù)，并以此為基礎(chǔ)建立三維模型，然后利用三維非連續(xù)變形分析軟件模擬危巖的移動(dòng)，最終得到邊坡穩(wěn)定性結(jié)果。金愛兵等[14]利用小型無(wú)人機(jī)獲取邊坡地形數(shù)據(jù)并建立三維模型，將三維模型與有限差分法軟件FLAC3D相結(jié)合，得到邊坡穩(wěn)定性結(jié)果。

綜上分析，現(xiàn)今學(xué)者的相關(guān)研究有兩方面不足，首先大多數(shù)的研究主要是利用無(wú)人機(jī)攝影技術(shù)建立三維模上的研究和基于點(diǎn)云數(shù)據(jù)對(duì)巖體進(jìn)行節(jié)理面信息提取算法的研究，而從航拍影像到節(jié)理面信息提取的整個(gè)過程的研究相對(duì)較少，即使有也大多利用商業(yè)軟件輔助完成。其次，當(dāng)測(cè)區(qū)存在大面積的植被時(shí)，三維建模的難度大大提升，有關(guān)采用無(wú)人機(jī)攝影測(cè)量技術(shù)對(duì)野外露頭巖質(zhì)邊坡進(jìn)行三維重構(gòu)的研究相對(duì)較少。現(xiàn)針對(duì)以上問題運(yùn)用無(wú)人機(jī)采用傾斜攝影測(cè)量技術(shù)對(duì)野外露頭巖質(zhì)邊坡進(jìn)行數(shù)據(jù)采集，基于增量式運(yùn)動(dòng)結(jié)構(gòu)恢復(fù)算法和深度圖融合的多視點(diǎn)立體視覺算法(multi-view stereo，MVS)對(duì)巖質(zhì)邊坡進(jìn)行三維模型重構(gòu)，然后再此基礎(chǔ)上進(jìn)行結(jié)構(gòu)面的識(shí)別，提取節(jié)理信息。根據(jù)上述過程，開發(fā)設(shè)計(jì)基于傾斜攝影的裂隙巖體邊坡三維模型重構(gòu)系統(tǒng)，為無(wú)人機(jī)在邊坡工程中的應(yīng)用提供了完整方法。

1 傾斜攝影技術(shù)

1.1 傾斜攝影技術(shù)要點(diǎn)

攝影測(cè)量是通過采集、處理、提取和表達(dá)調(diào)查區(qū)域的影像信息來(lái)獲取現(xiàn)場(chǎng)信息的方法。在傳統(tǒng)攝影測(cè)量中，一般要通過測(cè)區(qū)的分布總圖布置拍攝點(diǎn)，然后手工拍攝，這樣不僅耗時(shí)還十分煩瑣，并且傳統(tǒng)的攝影測(cè)量主要采用的是正直影像，僅對(duì)地形地貌上部區(qū)域有較好識(shí)別的信息，而無(wú)法獲取沿高度方向的地貌信息。傾斜攝影的出現(xiàn)，彌補(bǔ)了正直影像的不足，通過同一架無(wú)人機(jī)實(shí)現(xiàn)垂直、傾斜等不同角度的同時(shí)拍攝，進(jìn)而獲取更加準(zhǔn)確的測(cè)區(qū)信息，其生成模型的精度是傳統(tǒng)攝影測(cè)量精度的2～3倍[15]。在傾斜攝影時(shí)一般需要對(duì)下視、前視、后視、左視和右視5個(gè)方向進(jìn)行拍攝，現(xiàn)采用單鏡頭小型無(wú)人機(jī)進(jìn)行傾斜拍攝，無(wú)人機(jī)通過規(guī)劃不同的航線進(jìn)行采集數(shù)據(jù)，得到測(cè)區(qū)的完整場(chǎng)景信息。為了確保得到的場(chǎng)景信息能滿足三維建模的要求，需要拍攝的影像之間具有較高的重疊度和一致的地面分辨率。其中重疊度可以通過規(guī)劃無(wú)人機(jī)航線的旁向和航向的交叉大小來(lái)控制，因此地面分辨率是能否準(zhǔn)確反應(yīng)出地表真實(shí)信息的重要參數(shù)之一。正直影像地面分辨率的計(jì)算公式為

(1)

式(1)中：GSD為正直影像時(shí)的地面分辨率；δ為像元尺寸；h為無(wú)人機(jī)拍攝高度；f為相機(jī)焦距。再結(jié)合無(wú)人機(jī)相機(jī)主光軸旋轉(zhuǎn)角度，可得到傾斜攝影時(shí)近點(diǎn)、中心的和遠(yuǎn)點(diǎn)的大地分辨率GSDj、GSDz和GSDy，公式為

(2)

(3)

(4)

式中：a和b分別為傾斜角度和視場(chǎng)角的一半，其中b的計(jì)算公式為

(5)

式(5)中：γ表示相機(jī)內(nèi)部圖像傳感器的大小在一個(gè)方向的計(jì)算。

傾斜攝影的幾何關(guān)系如圖1所示，其中T′為影像中的一個(gè)點(diǎn)，T為T′在的實(shí)際場(chǎng)景中地面的對(duì)應(yīng)點(diǎn)，PP′為影像中的主光軸與像平面的交點(diǎn)，β為T和傾斜攝影中心O的連線和ON豎線的夾角。從圖1中可以看出，影響大地分辨率的因素主要包括無(wú)人機(jī)拍攝的高度、焦距、像素值、傾角等，其中焦距和像素值是由無(wú)人機(jī)攜帶的相機(jī)確定的；飛行高度越低，大地分辨率的值越小，則精度越高；結(jié)合式(4)和圖1可知，當(dāng)傾角越大時(shí)，遠(yuǎn)點(diǎn)的分辨率的數(shù)值越小且遠(yuǎn)點(diǎn)和近點(diǎn)分辨率差值越大。

1.2 數(shù)據(jù)采集

航測(cè)區(qū)域位于石家莊市井陘縣南障城鎮(zhèn)，周圍群山環(huán)繞，地形復(fù)雜。測(cè)區(qū)長(zhǎng)度為437.0 m，寬度為232.9 m，邊坡高140.0 m。收集數(shù)據(jù)選用大疆御2專業(yè)版小型無(wú)人機(jī)，無(wú)人機(jī)及地面遙控器如圖2所示，無(wú)人機(jī)飛行器配備先進(jìn)的全方位視覺系統(tǒng)及紅外傳感系統(tǒng)，對(duì)周圍障礙物有感應(yīng)識(shí)別的功能，這樣就使得無(wú)人機(jī)可以安全完成飛行任務(wù)。無(wú)人機(jī)上搭載高精度相機(jī)，并且配合高穩(wěn)定性云臺(tái)，足以滿足工程的拍攝要求。地面遙控器采用高清圖傳技術(shù)，可以控制無(wú)人機(jī)在最大8 km通訊距離內(nèi)完成飛行任務(wù)，無(wú)人機(jī)技術(shù)參數(shù)見表1。結(jié)合無(wú)人機(jī)操作端自帶的 DJI Pilot軟件，可以實(shí)現(xiàn)對(duì)飛行航線的規(guī)劃，得到垂直影像和傾斜影像。本文中在測(cè)區(qū)總共設(shè)置五條航線，第一條航線鏡頭豎直向下，其他航線鏡頭傾角60°，飛行航線如圖3所示，每種顏色代表一條完整的航線，傾斜攝影具體飛行軌跡是無(wú)人機(jī)從航線1起點(diǎn)出發(fā)，飛到指定高度完成航線1的任務(wù)，然后從航線1的最后位置飛向航線2的起點(diǎn)繼續(xù)拍攝任務(wù)，直到所有航線任務(wù)完成回到航線1起點(diǎn)降落，由于邊坡最高140 m，故飛行高度設(shè)置為200 m。在傾斜攝影測(cè)量過程中，為了能得到詳盡的數(shù)據(jù)，航向重疊度設(shè)置為80%，旁向重疊度設(shè)置為70%，等距離拍照，具體航線規(guī)劃設(shè)置和航線參數(shù)信息如表2所示。

圖3 航線規(guī)劃圖Fig.3 Route plan

表1 無(wú)人機(jī)技術(shù)參數(shù)Table 1 Technical parameters of UAV

表2 航線參數(shù)Table 2 Route parameters

2 三維模型重建及節(jié)理信息識(shí)別

2.1 基本流程

由于小型無(wú)人機(jī)搭載相機(jī)相幅較小，且對(duì)旁向重疊度和航向重疊度要求很高，得到的航片照片數(shù)量較多，為了滿足能同時(shí)識(shí)別較多的照片，本文中選擇基于增量式運(yùn)動(dòng)結(jié)構(gòu)恢復(fù)算法(structure-from-motion，SFM)[16]和深度圖融合的MVS算法對(duì)巖質(zhì)邊坡進(jìn)行三維模型重構(gòu)，在此基礎(chǔ)上進(jìn)行節(jié)理信息數(shù)字識(shí)別。具體技術(shù)流程見圖4。

圖4 技術(shù)流程圖Fig.4 Technology flow chart

在開源軟件COLMAP的框架的基礎(chǔ)上，對(duì)其進(jìn)行二次開發(fā)，利用了COLMAP軟件開源庫(kù)里的增量式SFM算法、深度圖融合的MVS算法和PCL軟件開源庫(kù)里的隨機(jī)抽樣一致算法(random sample consensus，RANSAC)、區(qū)域生長(zhǎng)法，將其重新編譯為既能進(jìn)行三維模型重構(gòu)又能進(jìn)行巖體結(jié)構(gòu)面識(shí)別提取的基于傾斜攝影的裂隙巖體邊坡三維模型重構(gòu)系統(tǒng)，實(shí)現(xiàn)了由傾斜攝影測(cè)量技術(shù)得到的數(shù)據(jù)生成三維模型和結(jié)構(gòu)面自動(dòng)識(shí)別的功能，軟件界面如圖5所示。

2.2 稀疏重構(gòu)

采用增量運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)(增量式的SFM算法)進(jìn)行稀疏重構(gòu)[17]，該算法通過相機(jī)的運(yùn)動(dòng)同時(shí)恢復(fù)相機(jī)的參數(shù)和場(chǎng)景結(jié)構(gòu)，其主要流程包括：特征點(diǎn)的提取、特征點(diǎn)的匹配和捆綁調(diào)整優(yōu)化相機(jī)姿態(tài)和場(chǎng)景結(jié)構(gòu)。

(1)特征點(diǎn)提取：小型無(wú)人機(jī)進(jìn)行拍攝任務(wù)時(shí)會(huì)因不可控因素使得拍攝的照片產(chǎn)生畸變，畸變的影響會(huì)影響到三維重構(gòu)的精度。選用SIFT(scale invariant feature transform)算法對(duì)畸變的照片進(jìn)行處理，能較好地完成特征點(diǎn)的提取任務(wù)[18]。SIFT算法是通過在尺度空間中計(jì)算極值點(diǎn)來(lái)檢測(cè)穩(wěn)定的點(diǎn)的特征。其中尺度空間定義為

L(x,y,σ)=G(x,y,σ)I(x,y)

(6)

(7)

式中:(x,y)為點(diǎn)的坐標(biāo);σ為尺度大小；G(x,y,σ)為尺度σ的可變高斯函數(shù)；I(x,y)為原圖像大小。SIFT點(diǎn)特征檢測(cè)在高斯差分尺度空間中進(jìn)行，目標(biāo)函數(shù)為

D(x,y,σ)=[G(x,y,kσ)-G(x,y,σ)]I(x,y)

(8)

式(8)中：k為比例系數(shù)。

利用高斯差分對(duì)數(shù)據(jù)進(jìn)行計(jì)算，得到取得極大值與極小值時(shí)的三維點(diǎn)，這些三維點(diǎn)為離散的潛在特征點(diǎn)，然后將這些潛在特征點(diǎn)進(jìn)行曲面擬合，擬合曲面的峰谷或者峰底即為最終特征點(diǎn)。通過本文研發(fā)的軟件檢測(cè)的特征點(diǎn)結(jié)果如圖6所示(見圖中紅色點(diǎn)為特征點(diǎn))。

圖6 SIFT特征檢測(cè)結(jié)果圖Fig.6 SIFT feature detection result graph

(2)特征點(diǎn)匹配：由于采用無(wú)人機(jī)進(jìn)行采集數(shù)據(jù)，數(shù)據(jù)中包含定位系統(tǒng)坐標(biāo)(global positioning system，GPS)和飛行姿態(tài)，故選擇空間匹配的方法?？臻g匹配的原理是首先給定一個(gè)查詢點(diǎn)，然后利用K最近鄰算法計(jì)算距離查詢點(diǎn)距離最近的點(diǎn)記為候選點(diǎn)，最后判斷候選點(diǎn)到查詢點(diǎn)距離是否超過了最大距離，沒超過最大距離的候選點(diǎn)加入特征點(diǎn)集，超過的最大距離的候選點(diǎn)剔除，完成特征匹配。

(3)增量的運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)：首先選取一對(duì)圖像進(jìn)行重建，再利用三角量測(cè)的方法對(duì)相機(jī)進(jìn)行旋轉(zhuǎn)、平移和三角化計(jì)算得到一個(gè)空間三維點(diǎn)，再把這個(gè)三維點(diǎn)重新投影到兩個(gè)視角中計(jì)算重投影誤差，誤差較小則重建成功。其中重投影誤差函數(shù)為

(9)

式(9)中：Pj為相機(jī)參數(shù)；Mi為三維點(diǎn)坐標(biāo)；mij為第j幅影像上的第i個(gè)觀測(cè)點(diǎn)；f(Pj,Mi)為觀測(cè)點(diǎn)對(duì)應(yīng)的投影點(diǎn)。

在此基礎(chǔ)上增加新的圖像重復(fù)上述過程，直至計(jì)算完所有照片，得到由空間三維點(diǎn)組成的稀疏場(chǎng)景模型，增量式SFM過程見圖7。

圖7 增量式SFM過程圖Fig.7 Incremental SFM process diagram

完成上述所有計(jì)算即可得到稀疏三維點(diǎn)云模型，通過本文研發(fā)的軟件得到的稀疏點(diǎn)云圖如圖8所示，圖8中紅色區(qū)域描述了無(wú)人機(jī)采樣時(shí)的位置及相機(jī)角度。

2.3 稠密三維重建

由圖8可以看出稀疏三維點(diǎn)云模型整體比較粗糙，不能直觀地反映出巖質(zhì)邊坡的狀態(tài)。故以稀疏三維點(diǎn)云模型為基礎(chǔ)，采用基于深度圖融合的方法，對(duì)巖質(zhì)邊坡進(jìn)行稠密三維重建。

基于深度圖融合的方法是通過對(duì)每幅影像重建深度圖進(jìn)而建立點(diǎn)云模型，由于不同影像拍攝的位置不同，看到的場(chǎng)景就會(huì)有所偏差，這些偏差即為視差，視差可以反映出物體在三維空間的深度。深度圖融合的MVS算法的基本流程：首先，在航測(cè)照片集中選取航測(cè)的照片，對(duì)每一張航測(cè)的照片選取若干張照片作為鄰近照片，用稀疏三維重建階段的結(jié)果作為判斷是否具有較高的相似性的標(biāo)準(zhǔn)，確保稠密三維重建的準(zhǔn)確性。然后對(duì)初始的稀疏特征點(diǎn)開始深度估計(jì)并且把這些稀疏特征點(diǎn)作為種子點(diǎn)，用PatchMatch算法進(jìn)行深度圖估計(jì)，該算法是通過像素點(diǎn)q的像素坐標(biāo)(qx,qy)和一個(gè)參數(shù)fq(aq,bq,cq)描述視差qd，公式為

qd=aqqx+bqqy+cp

(10)

PatchMatch算法中，每個(gè)像素初始化深度值是指用隨機(jī)數(shù)為fq的三個(gè)參數(shù)設(shè)置初始值，然后把每個(gè)像素的深度信息傳遞到鄰近像素，根據(jù)匹配代價(jià)判斷是否更新鄰近像素深度信息，最后不斷循環(huán)迭代得到場(chǎng)景的深度信息。該過程中深度信息優(yōu)化的數(shù)學(xué)模型為

(11)

式(11)中：IR為參考影像深度值；Ik為第k個(gè)視角的深度值；ck為第k個(gè)視角的深度尺度因子；(i,j)為位置坐標(biāo)。

通過上述步驟，計(jì)算出每幅照片的深度圖，如圖9所示(不同顏色代表巖質(zhì)邊坡不同位置距離無(wú)人機(jī)的不同距離)。當(dāng)計(jì)算出數(shù)據(jù)集中所有的圖的深度圖，把他們相互融合即可得到高程圖，如圖10所示，從高程圖中可以看出，最高處的高程為140 m與實(shí)際測(cè)量值吻合。

圖9 深度圖及其原圖Fig.9 Depth map and its original

圖10 高程圖Fig.10 Digital elevation map

最后，將深度圖進(jìn)行融合就能得到稠密三維模型及高程圖。通過本文研發(fā)的軟件得到稠密點(diǎn)云三維模型如圖11所示，可以看出稠密三維重建的模型能較好地還原巖質(zhì)邊坡。

圖11 稠密點(diǎn)云圖和航拍全景圖Fig.11 Dense point cloud and aerial panorama

2.4 節(jié)理信息識(shí)別

經(jīng)上述計(jì)算，得到了巖質(zhì)邊坡三維模型，模型中包含著空間幾何特征信息，因此可以對(duì)三維點(diǎn)云模型進(jìn)行結(jié)構(gòu)面的識(shí)別與提取相關(guān)信息。

結(jié)構(gòu)面識(shí)別過程包括兩個(gè)方面，分別為結(jié)構(gòu)面的區(qū)域搜索和結(jié)構(gòu)面的平面擬合。采用區(qū)域生長(zhǎng)法[19]進(jìn)行結(jié)構(gòu)面區(qū)域搜索，由于輸入的數(shù)據(jù)為稠密點(diǎn)云模型，其三維特征點(diǎn)的數(shù)目龐大，使得計(jì)算時(shí)間過長(zhǎng)且效率低下，所以首先減少三維點(diǎn)云的數(shù)據(jù)且保持三維點(diǎn)云的形狀特征不變，該過程就是對(duì)稠密點(diǎn)云模型進(jìn)行下采樣，下采樣后把候選點(diǎn)中曲率最小的點(diǎn)作為初始生長(zhǎng)點(diǎn)，搜索鄰域，比較鄰域中點(diǎn)的法向量和初始生長(zhǎng)點(diǎn)法線之間的夾角與閾值的大小，小于閾值的點(diǎn)保留到聚類數(shù)據(jù)，再比較聚類數(shù)據(jù)中的點(diǎn)的曲率與曲率閾值的大小，小于閾值的點(diǎn)并入到初始生長(zhǎng)點(diǎn)的區(qū)域，并將該點(diǎn)作為下一個(gè)生長(zhǎng)點(diǎn)向外生長(zhǎng)，直至到?jīng)]有可以歸并的點(diǎn)為止。

經(jīng)過結(jié)構(gòu)面的區(qū)域搜索后可以得到各個(gè)結(jié)構(gòu)面的點(diǎn)集合，然后利用RANSAC(random sample consensus)進(jìn)行巖體結(jié)構(gòu)面擬合和轉(zhuǎn)換[20]。在用RANSAC算法時(shí)，首先要選取有效的初始數(shù)據(jù)，在此基礎(chǔ)上不斷迭代增加數(shù)據(jù)集，最后選取數(shù)據(jù)最多的模型為最終模型。RANSAC平面擬合的流程如下。

(1)在三維點(diǎn)云數(shù)據(jù)集中選取初始點(diǎn)云數(shù)據(jù)，計(jì)算其對(duì)應(yīng)平面方程為

Ax+By+Cz+D=0

(12)

式(12)中：A、B、C、D為平面方程參數(shù)。

(2)計(jì)算所有點(diǎn)至該平面的代數(shù)距離di，選取閾值d，若di

(3)重復(fù)以上步驟，根據(jù)數(shù)據(jù)集計(jì)算出迭代結(jié)束參數(shù)，根據(jù)該迭代結(jié)束參數(shù)判斷迭代是否停止。

在開源軟件COLMAP中加入了PCL編程庫(kù)中Sample-consensus模塊,實(shí)現(xiàn)了RANSAC平面擬合。Sample-consensus模塊中模型定義為平面模型，模型設(shè)置了4個(gè)參數(shù)normal_x、normal_y、normal_z和d分別對(duì)應(yīng)著平面方程參數(shù)A、B、C和D。本文中利用基于傾斜攝影的裂隙巖體邊坡三維模型重構(gòu)系統(tǒng)，得到的結(jié)構(gòu)面識(shí)別圖如圖12所示。

圖12 結(jié)構(gòu)面識(shí)別圖Fig.12 Structural plane identification map

由于三維模型的拓?fù)浣Y(jié)構(gòu)是基于無(wú)人機(jī)傾斜攝影拍攝的照片構(gòu)建的，其坐標(biāo)系和客觀世界坐標(biāo)系相同，故可以根據(jù)空間解析幾何原理可以推導(dǎo)出結(jié)構(gòu)面產(chǎn)狀的關(guān)系[21]，推導(dǎo)關(guān)系如圖13所示，表達(dá)式為

圖13 結(jié)構(gòu)面產(chǎn)狀幾何關(guān)系示意圖Fig.13 Schematic diagram of the geometric relationship of structural planes

(13)

(14)

(15)

式中：θ和φ分別代表結(jié)構(gòu)面的傾角和傾向;ω為同組結(jié)構(gòu)面之間的間距。

從圖13可以看出，該巖質(zhì)邊坡有多組結(jié)構(gòu)面，利用研發(fā)的軟件得到節(jié)理信息數(shù)據(jù)，對(duì)巖質(zhì)邊坡識(shí)別出125個(gè)結(jié)構(gòu)面數(shù)據(jù)，給出部分展示數(shù)據(jù)如圖5(b)所示。整理所有數(shù)據(jù)剔除錯(cuò)誤結(jié)構(gòu)面，得到7組顯著結(jié)構(gòu)面和若干個(gè)單獨(dú)結(jié)構(gòu)面，由于內(nèi)容有限，僅提供顯著結(jié)構(gòu)面數(shù)據(jù)，平面方程參數(shù)見表3。為證明獲得的巖體結(jié)構(gòu)面傾向、傾角的可靠性，用人工現(xiàn)場(chǎng)測(cè)量的方法進(jìn)行比較，如表4所示，可以看出對(duì)于7組結(jié)構(gòu)面的產(chǎn)狀，用兩種方法計(jì)算的最大誤差為3.19，具有一致性。根據(jù)表4中的結(jié)構(gòu)面產(chǎn)狀數(shù)據(jù)，將點(diǎn)云劃分成有限區(qū)域，得到結(jié)構(gòu)面劃分圖，展示了巖質(zhì)邊坡結(jié)構(gòu)面的具體分布狀況，如圖14所示。

圖14 結(jié)構(gòu)面劃分圖Fig.14 Structural plane division diagram

表3 巖體結(jié)構(gòu)面數(shù)據(jù)Table 3 Rock mass discontinuity data

表4 巖體結(jié)構(gòu)面產(chǎn)狀結(jié)果比對(duì)Table 4 Comparison of rock mass discontinuity occurrence results

3 邊坡數(shù)值模擬計(jì)算

小型無(wú)人機(jī)對(duì)巖質(zhì)邊坡采集好數(shù)據(jù)，經(jīng)過基于傾斜攝影的裂隙巖體邊坡三維模型重構(gòu)軟件計(jì)算，得到地形網(wǎng)格模型及節(jié)理信息數(shù)據(jù)。在此基礎(chǔ)上，利用FLAC3D有限差分法軟件進(jìn)行數(shù)值模擬。數(shù)值模擬中采用Hoek-Brown準(zhǔn)則[22]，其表達(dá)式為

(16)

式(16)中：σ1、σ3分別為巖體最大主應(yīng)力、最小主應(yīng)力；σc為巖體抗壓強(qiáng)度；mb、s、α為巖體的Hoek-Brown常量。

3.1 三維模型建立

基于傾斜攝影的裂隙巖體邊坡三維模型重構(gòu)軟件可以提取點(diǎn)云的X、Y、Z各坐標(biāo)數(shù)據(jù)，利用勞德內(nèi)算法還原地形得到邊坡地形圖，以“ply”格式進(jìn)行保存，如圖15所示。利用Rhino軟件對(duì)地形網(wǎng)格模型進(jìn)行優(yōu)化，生成“.f3grid”格式文件，然后把三維模型導(dǎo)入FLAC3D，三維有限元模型如圖16所示，三維模型長(zhǎng)400 m、寬400 m、高140 m。

圖15 地形網(wǎng)格圖Fig.15 Topographic grid

圖16 三維有限元模型Fig.16 Three dimensional finite element model

3.2 三維模型參數(shù)

在進(jìn)行邊坡穩(wěn)定性分析時(shí)，巖土體均采用Hoek-Brown本構(gòu)模型，計(jì)算時(shí)除了模型上表面以外，對(duì)各邊界均處采用固定端約束。式(16)中mb、s、α三個(gè)常量通常利用地質(zhì)強(qiáng)度指標(biāo)(geological strength index，GSI)對(duì)巖體參數(shù)進(jìn)行確定的方法，其表達(dá)式為

(17)

(18)

(19)

式中：GSI為地質(zhì)強(qiáng)度指標(biāo)；D為巖體擾動(dòng)因子，其取值范圍為0.0～1.0，現(xiàn)場(chǎng)無(wú)擾動(dòng)巖體為0，而非常擾動(dòng)巖體為1.0；mi為巖體霍克布朗常數(shù)，可以通過查表5獲得。研究區(qū)域?yàn)橹械壬皫r，現(xiàn)場(chǎng)無(wú)擾動(dòng)巖體，查表5可知式(19)中mi=19，D=0。

表5 由巖石類型所決定的常數(shù)miTable 5 Constant mi determined by rock type

地質(zhì)強(qiáng)度指標(biāo)GSI 與巖體破裂程度和結(jié)構(gòu)面的狀態(tài)有關(guān)。從上述三維重構(gòu)及節(jié)理信息數(shù)字識(shí)別得出的結(jié)果出發(fā)，選用巖體體積節(jié)理數(shù)JV和結(jié)構(gòu)面條件因子JC對(duì)地質(zhì)強(qiáng)度指標(biāo)GSI進(jìn)行計(jì)算。根據(jù)已有的節(jié)理數(shù)據(jù)推導(dǎo)出適用于本研發(fā)軟件的巖體體積節(jié)理數(shù)的公式。

對(duì)于巖體體積節(jié)理數(shù)JV的統(tǒng)計(jì)，通常計(jì)算公式[22]為

(20)

式(20)要求每組結(jié)構(gòu)面的數(shù)量?jī)蓚€(gè)及兩個(gè)以上，才能算出結(jié)構(gòu)面的間距。在實(shí)際測(cè)量中存在單個(gè)結(jié)構(gòu)面，忽略它們將導(dǎo)致測(cè)量結(jié)果不準(zhǔn)確，假設(shè)這些單個(gè)結(jié)構(gòu)面的間距取各組結(jié)構(gòu)面平均間距累加的平均值，在此基礎(chǔ)上提出了計(jì)算巖體體積節(jié)理數(shù)JV的公式，即

(21)

式(21)中：N為單獨(dú)的結(jié)構(gòu)面?zhèn)€數(shù)；A為測(cè)區(qū)面積。

根據(jù)現(xiàn)場(chǎng)勘查可以知道本測(cè)區(qū)面積A=8 000 m2，結(jié)合表3和表4中的數(shù)據(jù)，利用式(21)，求得JV=25.62。結(jié)合Sonmez等[23]得出的地質(zhì)強(qiáng)度指標(biāo)表，如圖17所示，地質(zhì)強(qiáng)度指標(biāo)表中可以看出結(jié)構(gòu)面條件因子JC的大小根據(jù)結(jié)構(gòu)面表面的狀況分為1～18，JC越大表明結(jié)構(gòu)面表面條件越好。根據(jù)如圖11中的航拍實(shí)景圖和稠密點(diǎn)云三維模型，結(jié)合采集數(shù)據(jù)時(shí)的現(xiàn)場(chǎng)勘查結(jié)果可知，該巖質(zhì)邊坡出露面風(fēng)化嚴(yán)重且有大小不一的坑洞，表面條件較差，由此估計(jì)出JC=6。根據(jù)JV和JC的值查圖17，得到GSI的值為52。為了驗(yàn)證結(jié)果的可靠性，根據(jù)Palmstrom的[24]公式進(jìn)行計(jì)算，求得JV=25.60，可見JV值的大小僅差0.02，兩者結(jié)果非常相近，證明了本公式的可靠性。

圖17 GSI量化表Fig.17 GSI quantitative table

利用式(17)～式(19)可計(jì)算出Hoek-Brown準(zhǔn)則中的參數(shù)mb、s和α，結(jié)果見表6。

表6 Hoek-Brown準(zhǔn)則參數(shù)Table 6 Hoek-brown criterion parameters

3.3 數(shù)值計(jì)算結(jié)果

基于強(qiáng)度折減法對(duì)巖質(zhì)邊坡的安全系數(shù)進(jìn)行計(jì)算，將巖質(zhì)邊坡土體的抗剪切強(qiáng)度參數(shù)逐漸減小，直到體系不平衡力與典型內(nèi)力比率小于10-5，認(rèn)為體系達(dá)到了力的平衡，即完成計(jì)算，利用model factor-of-safety命令即可實(shí)現(xiàn)該過程。在求解過程中FLAC3D軟件能自動(dòng)搜索出巖質(zhì)邊坡最危險(xiǎn)的位置，分析其在臨界狀態(tài)下的穩(wěn)定性。計(jì)算結(jié)果表明：圖18為臨界狀態(tài)下巖質(zhì)邊坡總體位移云圖，可以看出，巖質(zhì)邊坡最危險(xiǎn)的位置在靠近山頂?shù)亩钙绿?，最大位移?1 cm。巖質(zhì)邊坡失穩(wěn)通常是沿著剪應(yīng)變?cè)隽孔畲蟮奈恢冒l(fā)生，因此需要找尋巖質(zhì)邊坡剪應(yīng)變的集中帶，圖19為臨界狀態(tài)下巖質(zhì)邊坡最大剪應(yīng)變?cè)隽吭茍D，可以看出，只有靠近坡腳處小部分的剪應(yīng)變?cè)隽枯^大，其他處剪應(yīng)變?cè)隽慷驾^小且變化均勻，說(shuō)明該邊坡整體穩(wěn)定性較好，但是在剪應(yīng)變?cè)隽枯^大的區(qū)域容易發(fā)生剪切破壞，與現(xiàn)場(chǎng)勘察結(jié)果一致。

圖18 臨界狀態(tài)下巖質(zhì)邊坡位移云圖Fig.18 Displacement cloud map of rock slope in critical state

圖19 臨界狀態(tài)下巖質(zhì)邊坡最大剪應(yīng)變?cè)隽吭茍DFig.19 Maximum shear strain increment cloud map of rock slope in critical state

4 結(jié)論

以南障城鎮(zhèn)巖質(zhì)邊坡為工程背景，利用小型無(wú)人機(jī)對(duì)邊坡進(jìn)行拍照，再通過開發(fā)的三維重構(gòu)軟件，得到點(diǎn)云圖、地形網(wǎng)格圖和結(jié)構(gòu)面產(chǎn)狀等。利用得到的數(shù)據(jù)，結(jié)合Hoek-Brown準(zhǔn)則，運(yùn)用數(shù)值模擬的手段分析邊坡的穩(wěn)定性。得到如下結(jié)論。

(1)使用小型無(wú)人機(jī)結(jié)合傾斜攝影測(cè)量技術(shù)，收集到完整有效的研究區(qū)域的數(shù)據(jù)，充分體現(xiàn)了小型無(wú)人機(jī)在地質(zhì)調(diào)查中的高效、便捷和準(zhǔn)確等特點(diǎn)。

(2)基于增量式運(yùn)動(dòng)結(jié)構(gòu)恢復(fù)算法和深度圖融合的MVS算法對(duì)巖質(zhì)邊坡進(jìn)行三維模型重構(gòu)，得到了測(cè)區(qū)的三維模型，該模型能很好的描述出實(shí)際巖質(zhì)邊坡。

(3)基于RANSAC算法和區(qū)域生長(zhǎng)法對(duì)模型的巖體結(jié)構(gòu)面進(jìn)行識(shí)別提取，得到了巖體結(jié)構(gòu)面的識(shí)別圖，提取到了結(jié)構(gòu)面的擬合方程、結(jié)構(gòu)面產(chǎn)狀和間距等參數(shù)，精度基本滿意。

(4)根據(jù)軟件輸出的結(jié)果，提出了一種新的巖體體積節(jié)理數(shù)計(jì)算公式，用該公式計(jì)算地質(zhì)強(qiáng)度指標(biāo)，得到結(jié)果誤差較小。結(jié)合軟件輸出的地形圖可快速建立數(shù)值分析模型，邊坡穩(wěn)定性計(jì)算結(jié)果與現(xiàn)場(chǎng)情況吻合。