李 翠 李杰林 張良兵 楊承業(yè) 徐繼業(yè) 周科平
(1.攀鋼集團礦業(yè)有限公司,四川 攀枝花 617000;2.中南大學資源與安全工程學院,湖南 長沙 410083;3.金屬礦山安全與健康國家重點實驗室,安徽 馬鞍山 243000)
結構面是巖體中強度較低的不連續(xù)面,是控制巖質邊坡穩(wěn)定性的重要因素,巖質邊坡的變形和破壞通常都是受控于巖體介質特性并沿著巖體結構面發(fā)生的。在礦巖地質條件、開采擾動、外部風化條件等影響下,結構面力學性能弱化,導致巖質邊坡發(fā)生滑坡、巖石塊體滑落等災害[1-3]。因此,精確獲取巖質邊坡的地形地貌及其巖體結構面參數是開展邊坡穩(wěn)定性分析與評價的基礎[4-5]。
目前,針對露天礦山邊坡地形測量及巖體結構面調查手段趨向于遠程化與智能化,傳統(tǒng)的人工接觸式測量工作效率低,且受限于復雜地形限制,已逐漸被三維激光掃描、攝影測量等非接觸式地質調查技術取代[6-8]。針對高陡邊坡地勢復雜險峻的特點,國內外許多學者利用無人機設備搭載激光掃描儀或攝像設備對邊坡進行了測繪與分析,陳昌富等[9]利用無人機貼近攝影獲取了邊坡的巖體結構特征;JIA 等[10]利用無人機對露天高陡邊坡進行了攝影測量,并利用數值分析方法獲取了邊坡巖體的穩(wěn)定性;KONG 等[11]基于無人機攝影測量結合結構運動(SfM)技術,研究了高分辨率邊坡數字模型生成方法;LIU 等[12]提出了一種將3D-DDA 與無人機攝影測量相結合的集成系統(tǒng),用于塊狀巖體邊坡的穩(wěn)定性評價。基于無人機測量技術的優(yōu)勢在于不受巖質邊坡的險要地形限制,可以充分獲取研究區(qū)域的圖像、結構及紋理特征,結合空中三角測量技術[12-13]可實現不同角度、不同方位的攝像信息的逆向建模,從而獲取測量邊坡區(qū)域的精細三維模型。
基于無人機測量所獲取的邊坡三維模型與點云數據,可以實現邊坡巖體結構的三維數字化還原。隨著計算機技術的發(fā)展,基于巖體結構數字模型開展室內巖體工程地質調查已成為主流的研究手段[14-17],相比于現場實地測量,可使得測量人員擁有更優(yōu)良的工作環(huán)境、更清晰的觀測條件,從而獲取數據量更大且更精確的巖體結構面數據。同時現有的計算機算法可輔助測量人員開展工程地質調查,進一步增強了作業(yè)的速度與可靠性。因此,基于無人機測量數據的結構面信息提取是建立從無人機邊坡測量到數字化建模,最后實現巖體結構數字化分析的關鍵。
為此,本研究基于無人機傾斜攝影測量與貼近攝影測量技術,進行露天巖質邊坡測量,獲取邊坡點云模型,并開展基于點云數據的邊坡巖體結構面識別研究,形成了一種基于無人機測量高精度點云數據的巖質邊坡結構面自動識別方法,最后將該方法應用于攀枝花鐵礦尖山礦區(qū)北部邊坡。
基于無人機精確測量是開展露天礦巖質邊坡結構面識別的前提,本研究采用了大疆精靈Phantom 4RTK 無人機開展了高精度的邊坡攝影測量,如圖1所示。

圖1 無人機飛控系統(tǒng)Fig.1 UAV flight control system
精靈Phantom 4RTK 是一款小型多旋翼高精度航測無人機,具備RKT 厘米級導航定位系統(tǒng)和高性能成像系統(tǒng),通過控制器可預設、調節(jié)機身的姿態(tài)、飛行軌跡、飛行速度以及攝像頭的角度、拍攝頻率,從而獲取不同角度、不同分辨率的攝影圖像。通過顯示器可實時獲取攝像頭的回傳影像以及飛機的飛行參數,從而控制無人機完成飛行測量作業(yè)。該無人機系統(tǒng)具有精度高、靈活、輕便以及測量速度快等特點,適用于小范圍的巖質邊坡精確測量。
巖質邊坡無人機測量及建模包括粗略建模與精細建模兩個部分,粗略建模即利用無人機對研究區(qū)域的總體形貌特征進行還原,獲取研究區(qū)的總體模型架構。精細建模則是在總體模型架構的基礎上,通過特征匹配技術對模型細節(jié)進行補充建模,從而獲取測量對象更精確的紋理結構特征。本研究通過DJI Terra軟件實現基于無人機巖質邊坡的粗、細建模,該款軟件基于空中三角測量技術對粗、細建模的攝影圖像數據進行耦合,實現更精細化的模型紋理及結構展現[13]。基于上述方法,本研究對無人機精確測量及建模的具體流程進行如下分析。
1.2.1 粗略建模
利用精靈Phantom 4RTK 無人機的傾斜攝影測量數據開展巖質邊坡的粗略建模。首先,根據測量區(qū)域的工程資料,通過現場勘察確定區(qū)內地形及周邊構筑物特征,開展無人機傾斜攝影的航線規(guī)劃,設置無人機的飛行測量參數。飛行測量參數主要包括無人機的飛行高度、相機傾斜角度、圖像重疊率等,精靈Phantom 4RTK 無人機即可根據航線及設置參數自動在邊坡上方保持一定的高度開展傾斜攝影測量,如圖2 所示。粗建模的主要目的是獲取測量邊坡區(qū)域的地形形貌特征,得到區(qū)域整體模型,進而得到邊坡區(qū)域的模型架構[18]。

圖2 無人機傾斜攝影測量示意Fig.2 Schematic of UAV oblique photogrammetry
1.2.2 精細建模
根據圖2 可以看出,對于有高低起伏的邊坡,無人機在傾斜攝影的固定航線上與邊坡表面的相對高度會不斷變化,所獲取的影像分辨率同樣存在差異,因此無人機傾斜攝影只能獲取邊坡的整體形貌,無法獲取邊坡復雜地形的高分辨率、精度均勻形態(tài),對于結構面參數識別研究,需要高精度的邊坡模型數據作為支撐。
利用無人機貼近攝影測量技術[19]進行無人機精細建模,如圖3 所示。無人機貼近攝影測量的具體操作步驟為:①利用無人機對邊坡表面進行近距離拍攝(5~50 m),同時保證無人機在豎直方向上的移動角度與邊坡角β一致,并且無人機攝影云臺的云臺角β1與邊坡角為余角關系(圖3(a)),使得無人機在沿航線移動時機身與坡面距離保持不變,并且云臺的攝影角度始終與坡面保持垂直;② 無人機沿著航線飛行攝影時的旁向重疊度不小于70%,航向重疊度不小于80%(圖3(b))。通過無人機貼近攝影測量技術可使無人機貼合邊坡表面的坡度獲取高精度、高分辨率的數據,能識別出更精細的邊坡坡面紋理與結構面信息,相對于邊坡粗略模型可展現出更精確、全面的巖體結構。

圖3 無人機貼近攝影測量Fig.3 UAV Nap-of-the-object Photography
基于本研究無人機粗建模以及細建模流程,可獲取研究區(qū)邊坡的精密點云數據,該類點云數據具有精準坐標并能精確地反映邊坡巖體的幾何特征,從而實現了利用點云數據通過計算機程序進行邊坡巖體結構面的幾何識別。
基于無人機精確測量獲得的點云數據可精確反映露天邊坡巖體結構面的幾何形貌特征,通過點云數據的平面分割算法即可獲取邊坡巖體結構面的信息。本研究主要采用區(qū)域生長法實現點云數據的平面分割,從而快速識別出巖質邊坡點云數據中的結構面信息。
2.1.1 點云數據的體素濾波
通過無人機高分辨率攝像獲取的點云數據通常較為密集且不均勻,使得點云數據平面分割計算的魯棒性變差,點云數據的體素濾波即是對點云進行體素化,通過體素采樣方式使得邊坡的點云數據規(guī)模降低且平整化,而整體邊坡的幾何及拓撲特性基本不發(fā)生變化[20],如圖4 所示。

圖4 點云數據體素濾波(單位:m)Fig.4 Voxel filtering of point cloud data
2.1.2 點云數據特征值提取
開展點云數據平面分割的前提是獲取點云數據的特征值[21],包括點云數據的法向量n及曲率特征φ。點云數據的法向量n是進行點云平面擬合的關鍵,而曲率特征φ可以確定平面的起伏狀態(tài),從點云數據中識別具有一定粗糙度的結構面時,選取合理的曲率特征φ值十分關鍵。對于點云數據中單個點云數據的法向量,本研究通過近鄰搜索算法將單個點云數據及其相鄰的k個點云數據作為一個合集,設合集中包含i個點云數據,則點云合集在最小二乘意義上所表示的平面P可表示為
式中,p為平面P的法向量;d為平面P到坐標原點的距離;ei為第i個點云數據的坐標。
根據上述分析,可認為由點云數據集合構成的平面P的法向量p,是由單個點云數據的法向量n通過主成分分析法(PCA)[22]求解得到,具體思路是其通過分解得到點云合集方差矩陣的特征值來計算法向量p。點云合集方差矩陣M及其特征值具有如下關系:

式中,λi與vj(j=0,1,2)分別為特征值及特征值向量。
若特征值λ0<λ1<λ2,則特征值中的最小值λ0即為單個點云數據的法向量n,同樣單個點云數據的曲率特征φ可由特征值計算得到:

2.1.3 區(qū)域生長法的點云數據平面切割與擬合
獲取點云數據的單個點云數據法向量n及曲率特征φ后,即可通過區(qū)域生長法[23]進行點云數據的平面切割,從而實現不同產狀結構面的分割識別,如圖5 所示。區(qū)域生長法是通過設置初始種子點,利用點云數據的法向量n及曲率特征φ作為判斷生長條件對種子點進行生長。若某一點云數據與種子點不屬于同一平面,則其法向量與種子點法向量的夾角以及該點的曲率必然大于某一閥值,通過設置這一閥值即可限定種子點的生長邊界,使得邊界內的點云屬于同一平面。

圖5 區(qū)域生長法原理Fig.5 Principle of region growing method
區(qū)域生長法的具體步驟為:
(1)根據點云數據的曲率值對點云進行排序,將曲率最小的點叫作初始種子點。
(2)設置一空的聚類序列C和空的種子點序列Q,將選好的初始種子點加入種子點序列,并搜索該種子點附近的k個鄰近點云數據。
(3)設置兩個點云數據之間的法向量夾角閾值a以及曲率閾值b,計算k個鄰近點云數據中每個點云法向量nk與種子點法向量n之間的夾角。當小于設定的法向量夾角閾值a時,首先將該點加入聚類序列C中,同時判斷該鄰域點的曲率值是否小于曲率閾值b,將小于曲率閥值的點云數據加入種子點序列Q中;然后在Q中重新選擇新的種子點,重復上述步驟,直到序列Q為空。
(4)按照曲率從小到大排序,取不存在于聚類序列C中的點云數據重新作為初始種子點,重復步驟(1)至(3),可實現基于區(qū)域生長法的點云數據平面分割。
2.1.4 結構面產狀計算
基于區(qū)域生長法的點云數據平面分割,可使得點云數據按照不同產狀的平面實現分離,利用主成分分析法即可得到各結構面的法向量D(A,B,C),通過平面法向量與結構面產狀的關系[24],即可計算出結構面傾角α與傾向β,公式為
一切仿佛是昨天,記憶是流淌的河,深入藏地,遇見他們讓我感到無比幸福,那些鮮活的面孔伴隨著樸素的名字,宛如河底多彩的石,閃動著美妙的色彩,縈繞在溫暖的思緒里。至今我常去甘南草原看看,想念他們成了慣性,每畫,總沉迷。
式中,T為角度常量,取值為
基于無人機傾斜攝影測量與貼近攝影測量獲取的巖質邊坡精細三維點云數據,結合上述點云數據的處理算法從而形成了巖質邊坡結構面自動識別方法,如圖6 所示。

圖6 巖質邊坡結構面自動識別方法Fig.6 Automatic identification method of rock slope discontinuities
攀枝花鐵礦位于攀枝花市北側,年產1 300 萬t礦石量,是攀鋼集團鐵礦石原料生產的主要基地之一。該礦包括朱家包包礦區(qū)與尖山礦區(qū)兩個開采區(qū)域,兩個礦區(qū)之間沿東西走向相連,前期均采用露天開采。目前尖山礦區(qū)已轉為地下開采,朱家包包礦區(qū)也轉入深凹露天開采,兩個礦區(qū)的最終露天邊坡高達600 m,為典型的高陡邊坡。由于節(jié)理裂隙的分布,在開采擾動、裂隙滲流、巖石風化作用等綜合因素影響下,露天采場邊坡時有發(fā)生滑坡及臺階垮塌等災害,對露天礦的安全生產造成了重要影響,如圖7 所示。

圖7 攀枝花鐵礦邊坡滑坡災害情況Fig.7 Landslide disaster of Panzhihua Iron Mine
本研究以尖山礦區(qū)北部邊坡(圖8)為例,開展基于無人機的精確測量作業(yè)。該區(qū)域主要巖性為輝長巖,邊坡角約為45°,由于風化、地質作用已產生了嚴重的臺階垮塌,若采用傳統(tǒng)的接觸式結構面測量方式,不僅存在較大的安全風險,并且測量人員難以到達上部區(qū)域,采用三維激光掃描儀器則難以找到合適的架站點進行測量,因此該區(qū)域更適用于利用無人機進行測量。
測量時,首先預設4 個控制點,利用精靈Phantom 4RTK 無人機的五向飛行進行邊坡整體區(qū)域的傾斜攝影,設置航向重疊度為80%,旁向重疊度為70%,無人機與邊坡的相對飛行高度設置為200 m。根據測量結果,利用DJI Terra 軟件得到研究區(qū)邊坡的整體模型。根據粗略模型及現場勘查所反映的研究區(qū)形貌,考慮到區(qū)域覆蓋的用電線路復雜,若采用線路規(guī)劃的自動飛行進行貼近測量難度較大,因此本研究通過手動平飛方式進行邊坡的貼近測量。無人機與邊坡的垂直距離保持在30 m 左右,無人機的云臺攝影角度設置為45°,可以看出,相對于邊坡的粗略模型,通過貼近測量得到的邊坡精細模型可以更細致地展現出邊坡的紋理與巖體結構等信息(圖8(a))。最后,基于無人機貼近測量得到的邊坡精細模型,并通過DJI Terra 軟件可輸出研究區(qū)邊坡的精密點云數據(圖8(b)),坐標N 代表正北方向,研究區(qū)垂直坡高為93 m,水平長度為268 m。

圖8 研究區(qū)模型構建Fig.8 Model establishment of the study area
3.3.1 邊坡整體模型精度驗證
基于4 個控制點對整體邊坡模型的精確度進行驗算,參考國家大地坐標系數學精度檢測規(guī)范[25],本研究通過誤差絕對值取平均來進行模型的精度檢驗,公式為
通過式(5)可計算出整體邊坡模型的4 個控制點在X、Y、Z方向上的平均誤差為ωx=0.93 cm、ωy=0.88 cm、ωz=1.03 cm,同時計算4 個控制點之間的平均距離誤差ωd=1.13 cm。可以看出,通過無人機測量得到的整體邊坡模型與實際模型坐標誤差約1 cm,可以滿足工程應用需求。
3.3.2 結構面識別算法精度驗證
由于測量人員無法靠近研究區(qū)域,為驗證結構面識別的精確度,選取邊坡整體模型中已噴漿加固且較為平整的臺階坡面對結構面識別的準確度進行分析,結果如圖9 所示。

圖9 結構面識別算法精度驗證Fig.9 Accuracy verification of structural plane recognition algorithm
該區(qū)域位于尖山礦區(qū)南部,由3 個坡面、1 388 m臺階和1 412 m 臺階組成,由于區(qū)域較為平整,且各個臺階與坡面之間的產狀相同(圖9(a)),可將坡面與臺階面視為平面。通過Python 編程進行結構面識別方法的功能實現,對驗證區(qū)域的平面進行識別與分割,獲得了研究區(qū)域的3 個坡面與2 個臺階平面的識別結果(圖9(b))。將識別擬合的平面產狀與現場羅盤測量結果進行比較(表2),可以看出采用區(qū)域生長法識別的平面產狀與羅盤測量結果誤差均在5°以內,可認為本研究采用的結構面識別方法精確性較好。

表2 結構面識別精度驗證結果Table 2 Accuracy verification results of structural plane recognition
值得注意的是,在驗證分析結果中,坡面與臺階之間轉角處的點云數據未被區(qū)域生長法所采用,主要原因是該區(qū)域點云的法向量與曲率變化過于明顯,從而導致不同平面之間的邊界線存在取樣誤差。為此,在后續(xù)結構面識別時應盡可能選取同一坡面內的點云數據進行分析。
通過對模型與方法的精度驗證,反映出本研究得到的邊坡模型與結構面識別方法可滿足研究區(qū)結構面調查需求,考慮到該區(qū)邊坡風化較為嚴重,選取上部臺階巖體較完整區(qū)域進行了結構面識別分析,結果如圖10 所示。采用上述結構面識別方法,在研究區(qū)共識別出1 536 條結構面(圖10(a))。通過對優(yōu)勢結構組進行編錄,結果表明,該區(qū)域輝長巖主要發(fā)育3 組優(yōu)勢結構面組:第1 組平均傾向為123.16°±22.64°,平均傾角為45.89°±18.25°;第2 組平均傾向為191.86°±25.24°,平均傾角為75.38°±23.19°;第3 組平均傾向為243.30°±18.29°,平均傾角為66.68°±23.76°(圖10(b))。進一步分析可知:研究區(qū)巖體結構較為發(fā)育,在巖石風化、地質作用以及開采擾動的多重因素影響下極易發(fā)生臺階的垮塌破壞,在工程支護中需要重點進行加固和防范。

圖10 研究區(qū)結構面識別與分組Fig.10 Identification and grouping of structural planes in the study area
(1)利用無人機傾斜攝影測量以及無人機貼近攝影測量技術,實現了基于無人機的露天礦巖質邊坡由粗到細、由整體到局部的精確測量,獲取了精細的邊坡巖體結構面點云數據模型。通過對點云數據進行體素濾波、特征提取等處理,利用區(qū)域生長法進行基于點云數據特征值的多平面分割,形成了一種基于點云數據的巖質邊坡結構面自動識別與提取方法,通過該方法識別的平面產狀結果驗證誤差均小于5°。
(2)根據所提出的巖質邊坡點云數據結構面快速識別方法,以攀枝花尖山礦區(qū)北部邊坡為例,利用無人機測量獲取了該區(qū)域的精細三維模型與點云數據。基于邊坡點云數據完成了邊坡巖體結構面識別,獲取了優(yōu)勢結構面的特征參數,為礦山邊坡穩(wěn)定性分析提供了基礎數據。
(3)提出了一種通過無人機測量獲取邊坡點云數據、采用計算機算法在巖質邊坡點云數據中自動識別結構面特征的研究思路,后期研究可以從提取結構面算法的魯棒性與準確性方面開展深入分析,通過機器學習對結構面識別樣本進行訓練,解決算法對點云數據中的平面邊界比較敏感的難題,從而獲取更符合工程實際的結構面數據。