梁 波 梁曉鵬 王 健 曲相屹

(1. 山東科技大學(xué) 測繪與空間信息學(xué)院 山東 青島 266590;2. 山東黃金礦業(yè)(萊州)有限公司焦家金礦, 山東 煙臺 261441)

0 引言

在礦山的整個生命周期中,采剝工程量始終是一個重要的數(shù)據(jù)依據(jù),因此，采剝工程量的計算也是測量驗收工作中的一項重要工作。在傳統(tǒng)的露天礦測量當(dāng)中,葉錦強[1]針對露天礦測繪中控制點經(jīng)常被破壞的問題,采用了全站儀自由設(shè)站法來解決該問題,雖然能夠很好地解決問題,但是該方法仍需要大量人力物力,而且測量作業(yè)時間較長。黃軍等[2]采用免像控?zé)o人機航測技術(shù),通過提高精度冗余和影像重疊度的方法來提高航測成果的精度,平面精度達到5.4 cm,高程精度達到8.7 cm。但是遮擋地區(qū)的影像仍會發(fā)生畸變。蔣鳳保等[3]利用無人機航攝系統(tǒng)對露天礦進行了數(shù)據(jù)采集，并制作了數(shù)字正射影像圖(Digital Orthophoto Map,DOM)、數(shù)字線劃地圖(Digital Line Graphic,DLG)等。

隨著測量手段的不斷進步和發(fā)展,三維激光掃描技術(shù)也開始應(yīng)用于礦山測量中。王森等[4]將尺度不變特征變換(Scale-invariant feature transform,SIFT)算法和經(jīng)典最近點迭代(Iterative Closest Point,ICP)算法結(jié)合，將站式掃描儀所獲得的兩期點云數(shù)據(jù)精確配準(zhǔn),計算得到采剝工程量。提高了露天礦驗收測量的效率和精度。解決了傳統(tǒng)測繪中單點地形圖測繪速度慢、測量密度低、計算方量不夠精確、測繪過程不安全等問題。但是該配準(zhǔn)算法過于復(fù)雜,不適合大規(guī)模應(yīng)用。潘勇等[5]運用三維激光掃描技術(shù)對露天礦山進行采礦權(quán)核查和動態(tài)儲量檢測,成果更為豐富、直觀、數(shù)字化。但是在點云顯示方面還有待提高。并且在點云數(shù)據(jù)預(yù)處理時耗時較長;人工參與較多,導(dǎo)致最終的數(shù)據(jù)處理結(jié)果受人為因素影響較大。

本文探討的是三維激光掃描技術(shù)在露天礦測繪中的應(yīng)用。在獲取的機載點云基礎(chǔ)上,采用布料模擬算法和人工去噪相結(jié)合的方法剔除點云數(shù)據(jù)噪聲點,減少了人工參與的情況。最后根據(jù)三角網(wǎng)體積計算法計算采剝工程量。

1 技術(shù)原理

1.1 機載激光雷達技術(shù)原理

機載激光雷達(Light Detection and Ranging,LiDAR)系統(tǒng)由激光測距儀(Laser Range Finder,LRF)、全球定位系統(tǒng)(Global Positioning System,GPS)、慣性傳感器(Inertial Measurement Unit,IMU)、計算機控制導(dǎo)航系統(tǒng)(Computer Controlled Navigation System,CCNS)、數(shù)據(jù)存儲單元、電荷耦合器件(Charge Coupled Device,CCD)相機組成。圖1為一個機載LiDAR系統(tǒng)的基本組成。其中,LRF用來發(fā)射和接收激光并計算距離,GPS用來確定掃描點的坐標(biāo),IMU用來測量機載LiDAR系統(tǒng)的航向角(H)、翻滾角(R)、俯仰角(P),CCNS用來控制數(shù)據(jù)通信和飛行器的導(dǎo)航,CCD相機用來同步采集地面的真彩色影像數(shù)據(jù)。

圖1 機載LiDAR系統(tǒng)的基本組成

1.2 布料模擬濾波算法

布料模擬濾波(Cloth Simulation Filtering,CSF)算法的計算過程可視為模擬布料的物理過程[6]。

計算機中模擬布料的過程就是海量離散節(jié)點構(gòu)成一個格網(wǎng)的過程,并且該格網(wǎng)中的節(jié)點相互聯(lián)系,節(jié)點之間由連接線進行連接。同時將各個格網(wǎng)點均設(shè)置為具有恒定質(zhì)量且無尺寸的質(zhì)點,因此，這種生成格網(wǎng)的方式也稱作質(zhì)點彈簧模型。如果想了解模擬的布料在重力作用下某一瞬時時刻的形狀,就需要知道各布料點在那一時刻的空間位置。假設(shè)各布料節(jié)點的位移只發(fā)生在垂直方向上,而且只受重力與鄰近節(jié)點的相互作用力。當(dāng)只考慮重力作用時,根據(jù)牛頓第二定律,由式(1)來確定布料點某一時刻的空間位置和那一時刻所受的作用力之間的關(guān)系：

(1)

式中,m為布料點的質(zhì)量;g為重力加速度常數(shù);X(t)為某一時刻格網(wǎng)節(jié)點的位置;Δt為時長。

當(dāng)知道時長和格網(wǎng)節(jié)點的初始位置時,通過式(1)就可以得到格網(wǎng)節(jié)點的當(dāng)前位置。為了減小粒子在倒置的表面反轉(zhuǎn)的影響,還需考慮由粒子內(nèi)部因素引起的位移。在格網(wǎng)節(jié)點中隨機選取2個相鄰的粒子,若兩者均可移動,那么就將兩者沿著相反方向移動相同的距離;若兩者均具有相同的高度,那么就不進行移動;若其中一個粒子不可移動,那么就移動另一個。位移量可由式(2)進行計算：

(2)

布料模擬濾波算法的詳細步驟如下：

第一步：在第三方軟件中采用人機交互的方式檢查原始點云數(shù)據(jù)中是否存在明顯的噪點。

第二步：將點云數(shù)據(jù)進行翻轉(zhuǎn)。

第三步：設(shè)置初始的布料格網(wǎng),并設(shè)定合適的格網(wǎng)分辨率和節(jié)點數(shù)量。將最初的“布料”位置設(shè)置在最高點上方。

第四步：將所有點云數(shù)據(jù)和格網(wǎng)粒子投影到同一個水平面上,遍歷每個節(jié)點粒子對應(yīng)的最近鄰點,并將其投影前的高程記錄下來。

第五步：計算每個格網(wǎng)粒子在重力影響下產(chǎn)生的位移,并將其與粒子對應(yīng)最鄰近點投影前的高程進行比較,若粒子的高程低于或者等于其投影前的高程,則把粒子的高度設(shè)置為其投影前的高程并將該粒子設(shè)置為固定點。

第六步：計算每個格網(wǎng)“粒子”在內(nèi)部驅(qū)動因素影響下所產(chǎn)生的位移。

第七步：重復(fù)上述第五步、第六步,直到所有粒子的最大高程變化達到預(yù)設(shè)值或者迭代次數(shù)超出閾值,則停止模擬過程。

第八步：計算點云數(shù)據(jù)與格網(wǎng)粒子之間的高度差異。

第九步：分離地面點與非地面點,若點云數(shù)據(jù)中的點與格網(wǎng)粒子之間的距離小于之前設(shè)置的閾值,則認(rèn)為其是地面點,反之則認(rèn)為其是非地面點。

1.3 Delaunay三角網(wǎng)生長算法原理

三角網(wǎng)生長算法的主要思想是先在點云數(shù)據(jù)中找出距離最短的兩點,使其構(gòu)成一條三角形初始邊,然后根據(jù)Delaunay三角網(wǎng)準(zhǔn)則確定第三個點生成第一個三角形,再對另外兩條邊繼續(xù)尋找第三點進行生長,向周圍點云區(qū)域擴展,重復(fù)此過程,最終構(gòu)成Delaunay三角網(wǎng)[7]。

步驟1：隨機選一點作為起始點,尋找與起始點最近的數(shù)據(jù)點,將兩點構(gòu)成初始邊,以這條邊作為基線;

步驟2：找出與基線構(gòu)成Delaunay三角形的最優(yōu)點形成最優(yōu)三角形;

步驟3：重復(fù)步驟1至步驟2,直到所有點都加入三角網(wǎng)格中。

1.4 三角網(wǎng)法計算原理

三角網(wǎng)法就是根據(jù)變化前后的兩期坐標(biāo)數(shù)據(jù),構(gòu)成相鄰的三角形,將這些三角形組成不規(guī)則三角網(wǎng)(Triangulated Irregular Network，TIN)。對所需區(qū)域采用三棱柱體積計算的方法來計算體積,最后累計求和得到最終所需總體積[8-9]，如式(3)所示：

(3)

式中,Vi為第i個三棱柱的體積;H1，H2，H3分別為三角形三個頂點的高程值;Si為三棱柱投影后的底面積。

根據(jù)計算原理可知,該方法具有很強的魯棒性,能夠適應(yīng)多樣的地形地貌。并且能夠最大程度還原原始的地形地貌。但該方法所需的數(shù)據(jù)量較大,運算時對計算機的配置要求較高[10-11]。

2 工程實例

該工程實例采用的是某一露天礦點云數(shù)據(jù)。該試驗區(qū)東西長約1.3 km,南北長約1 km,該地區(qū)礦坑呈現(xiàn)階梯狀,從北向南依次遞減,有多個臺階,每個臺階大概有20 m的高度,等高線地形如圖2所示。大型礦車具有很大的視線盲區(qū),使用傳統(tǒng)的測繪方法對工作人員具有一定的安全威脅，且傳統(tǒng)測量方式周期較長不能滿足整個礦區(qū)的測量驗收需求。因此，本項目采用激光機載雷達進行施測。兩期數(shù)據(jù)間隔半個月。每期數(shù)據(jù)測量時間均在半天左右,相比于傳統(tǒng)測量,外業(yè)測量時間大大減少。

圖2 測區(qū)等高線地形圖

2.1 點云數(shù)據(jù)的去噪與濾波

本文采用人機交互的方式對遠離地表的離群點進行去噪,之后采用布料模擬算法將地表和近地表的“噪聲點”進行分離。

選取一小部分地勢較為平坦的點云進行濾波處理,設(shè)置布料分辨率為1,處理結(jié)果如圖3所示。圖3標(biāo)記部分為近離群“噪聲點”,通過布料模擬算法進行濾波后能夠很好地分離地面點和非地面點。將該算法用于露天礦點云中,露天礦呈階梯狀地形起伏比較大,會產(chǎn)生過度濾波的結(jié)果,如圖4所示,該算法將坡度較大的區(qū)域當(dāng)成“噪聲點”。主要分為坡度點云和遠離群點點云,故將這部分點云結(jié)果進行簡單的手動去噪。最后將兩次處理結(jié)果合并得到完整的地表。最終兩期點云數(shù)據(jù)的去噪效果如圖5所示。

(a)濾波前

圖4 分離的非地面點

(a)第一期

2.2 基于三角網(wǎng)生長法的曲面建模

本文采取隨機采樣法提取平均點云密度為2 m的點云生成三角網(wǎng),局部的三角網(wǎng)細節(jié)如圖6所示。

圖6 三角網(wǎng)模型局部細節(jié)

將三角網(wǎng)生成結(jié)果與實測結(jié)果對比,發(fā)現(xiàn)誤差主要為坡度較大地區(qū)。故將隨機采樣提取的平均點密度設(shè)置為1 m,結(jié)果發(fā)現(xiàn)坡度較大地區(qū)與實測地區(qū)最大誤差為0.3 m?？蓾M足精度要求,最終得到兩期數(shù)據(jù)的地形模型。如圖7所示。

(a)第一次掃描

2.3 采剝工程量的計算

在地形起伏比較大的地區(qū),采用三角網(wǎng)法計算采剝工程量時,需將測點間距設(shè)置為3～10 m,這樣才能保證計算結(jié)構(gòu)偏差為0～0.44%[12-13]。而機載點云數(shù)據(jù)則是能夠?qū)y點間距保持在分米級。但是點的數(shù)量越多，生成的三角網(wǎng)越密集,所耗費的時間越長。

本文綜合考慮計算偏差和耗費時間,將測點密度設(shè)置為1 m左右,建立三角網(wǎng)模型。將兩期的三角網(wǎng)模型進行對比，以重疊部分的高程值設(shè)置為0，最終得到差量模型,如圖8所示。

圖8 差量模型

根據(jù)三角網(wǎng)體積計算法,對差量模型進行計算,得到采剝工程量。將未進行布料模擬算法處理的兩期原始點云同樣進行三角網(wǎng)模型建立和三角網(wǎng)體積計算,得到采剝工程量。與露天礦實際采剝工程量進行對比分析,最終結(jié)果如表1所示。由結(jié)果可知,采剝工程量與真實值相比有所減少,主要原因是采用布料模擬算法和人工去噪將近地表的部分點云當(dāng)成噪聲點去除了。運用公式(4)進行精度評定可知,相對精度由原來的0.44%提高到了0.19%,使得結(jié)果更接近真實值。

(4)

表1 對比分析

式中,σ表示相對精度;m表示測量值;n表示真實值。

3 結(jié)束語

本文在點云數(shù)據(jù)計算采剝工程量過程中,運用布料模擬算法和Delaunay三角網(wǎng)生長算法,通過設(shè)定一定的參數(shù)對三角網(wǎng)進行構(gòu)網(wǎng)，使其能夠滿足露天礦測繪的要求,通過工程實例證明該方法所得結(jié)果更接近于真實值。同時使用三維激光掃描技術(shù)，提高了數(shù)據(jù)采集的效率,更直觀地表達了被測物體的三維信息。但是在點云數(shù)據(jù)預(yù)處理過程中,用布料模擬算法進行去噪,會將坡度較大的地區(qū)誤認(rèn)為“噪聲點”,導(dǎo)致過度去噪,從而增加后期人工去噪的工作量和難度。后續(xù)對去噪方法仍需要繼續(xù)研究。