利用激光點(diǎn)云數(shù)據(jù)進(jìn)行露天礦采剝工程量的計(jì)算方法

梁 波 梁曉鵬 王 健 曲相屹

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

0 引言

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

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

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

1 技術(shù)原理

1.1 機(jī)載激光雷達(dá)技術(shù)原理

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

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

1.2 布料模擬濾波算法

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

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

(1)

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

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

(2)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(3)

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

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

2 工程實(shí)例

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

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

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

本文采用人機(jī)交互的方式對(duì)遠(yuǎn)離地表的離群點(diǎn)進(jìn)行去噪,之后采用布料模擬算法將地表和近地表的“噪聲點(diǎn)”進(jìn)行分離。

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

(a)濾波前

圖4 分離的非地面點(diǎn)

(a)第一期

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

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

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

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

(a)第一次掃描

2.3 采剝工程量的計(jì)算

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

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

圖8 差量模型

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

(4)

表1 對(duì)比分析

式中,σ表示相對(duì)精度;m表示測(cè)量值;n表示真實(shí)值。

3 結(jié)束語

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