基于密集匹配點(diǎn)云的露天礦臺(tái)階線提取方法研究

王繼野

(中鐵十九局集團(tuán)礦業(yè)投資有限公司,北京 100161)

露天礦臺(tái)階線是礦山開采過程中產(chǎn)生的重要邊坡跡線,是露天礦測(cè)量的主要對(duì)象之一。傳統(tǒng)的露天礦臺(tái)階線提取主要是利用全站儀[1]、GPS[2]等測(cè)量方法。隨著測(cè)繪科學(xué)技術(shù)的發(fā)展,出現(xiàn)了利用數(shù)字高程模型(Digital Elevation Model,DEM)[3]、三維激光點(diǎn)云數(shù)據(jù)、無人機(jī)傾斜攝影測(cè)量[4]等技術(shù)來提取露天礦臺(tái)階線的現(xiàn)代測(cè)量方法。其中,基于激光掃描與攝影測(cè)量技術(shù)的臺(tái)階線提取方法主要是利用現(xiàn)代化測(cè)繪手段構(gòu)建三維點(diǎn)云模型,并在點(diǎn)云模型中提取臺(tái)階特征線。

目前主要的點(diǎn)云特征線提取方法可以歸納為3種:① 基于點(diǎn)云曲率或者法向量[5]。CHEN等[6]根據(jù)向量分布和聚類的性質(zhì)提取點(diǎn)云特征點(diǎn),采用三次B樣條曲線擬合算法得到點(diǎn)云特征線。該方法適用于規(guī)則物體特征線的提取,對(duì)于一些地形起伏的場(chǎng)景點(diǎn)云,該方法難以提取出相應(yīng)的特征點(diǎn)。② 基于鄰近點(diǎn)云特征相似性。FU等[7]提出了一種基于點(diǎn)空間幾何結(jié)構(gòu)的特征線提取方法,提高了點(diǎn)云特征線提取的準(zhǔn)確性和效率;WANG等[8]提出了一種利用激光線結(jié)構(gòu)光獲取三維數(shù)據(jù)的特征提取方法,并被應(yīng)用于鐵路特征線提取?？傮w來看,基于鄰近點(diǎn)云特征相似性提取點(diǎn)云時(shí),對(duì)點(diǎn)云顏色等信息要求較高,點(diǎn)云間相似性度量標(biāo)準(zhǔn)的選擇難以把握,故而較難廣泛使用。③基于點(diǎn)云分割的特征線提取。NI等[9]基于區(qū)域增長(zhǎng)和模型擬合的混合方法實(shí)現(xiàn)點(diǎn)云特征線跟蹤;GILANI等[10]通過點(diǎn)云分割的方式完成了建筑物特征線提取,解決了在遮擋和低采樣點(diǎn)情況下的特征線提取難題。該類基于點(diǎn)云分割的特征線提取操作需要在分割的基礎(chǔ)上實(shí)現(xiàn),對(duì)于一些表面不規(guī)則的物體分割較為困難,難以保證特征線的提取精度,因而該方法應(yīng)用具有一定的局限性。

上述特征線提取方法雖然可以較好地完成點(diǎn)云特征線提取,但是大部分主要是針對(duì)工業(yè)零部件、建筑物點(diǎn)云,該類點(diǎn)云模型形狀較為規(guī)整,特征線以直線為主。目前研究中較少針對(duì)露天礦區(qū)不規(guī)則臺(tái)階線提取方法的研究,東北大學(xué)王植等[4]曾提出一種顧及鄰域幾何屬性的三維邊緣檢測(cè)與曲率指數(shù)加權(quán)方法,實(shí)現(xiàn)了臺(tái)階線的自動(dòng)化提取,但提取出的臺(tái)階線出現(xiàn)了少量不連續(xù)現(xiàn)象。由于露天礦臺(tái)階線具有不規(guī)則性且礦區(qū)地形復(fù)雜,因而利用點(diǎn)云數(shù)據(jù)進(jìn)行露天礦臺(tái)階線提取具有一定的難度。鑒于此,本研究首先基于自適應(yīng)鄰域內(nèi)高程變化特征提取臺(tái)階線特征點(diǎn),然后采用三次B樣條曲線擬合出露天礦臺(tái)階線,取得了較好的提取效果。

1 露天礦臺(tái)階線提取原理

1.1 點(diǎn)云局部特征判定

對(duì)于點(diǎn)云樣本數(shù)據(jù)集而言,空間維度為3,此時(shí),點(diǎn)云樣本數(shù)據(jù)集S= { s1,s2,…,sj,…,sn}?R3,j=1,2,…,n,n為樣本點(diǎn)數(shù)量。取數(shù)據(jù)集S中任意樣本點(diǎn)s,其K近鄰域點(diǎn)集Ns= { s1,s2,…,sb,…,sk},b=1,2,…,k,點(diǎn)s的協(xié)方差矩陣可表示為[4]

由于協(xié)方差矩陣M為對(duì)稱正定陣,因此其特征值λ0、λ1、λ2具有非負(fù)性,同時(shí)可用協(xié)方差矩陣M對(duì)應(yīng)的特征值對(duì)點(diǎn)周邊鄰域的維度特性進(jìn)行描述,即:① 當(dāng) λ0＞＞λ1、λ2＞0 時(shí),局部鄰域呈線性特征分布;② 當(dāng) λ0、λ1＞＞λ2＞0 時(shí),局部鄰域呈平面特征分布;③ 當(dāng) λ0≈λ1≈λ2時(shí),局部鄰域呈三維曲面特征分布。

1.2 自適應(yīng)鄰域內(nèi)高程變化的特征點(diǎn)提取

露天礦臺(tái)階線包括臺(tái)階坡頂線和臺(tái)階坡底線[11]。如圖1所示,臺(tái)階上部平盤與下部平盤均是一個(gè)平面,平盤內(nèi)各處的高程值大致相等,臺(tái)階坡面高程呈現(xiàn)漸變趨勢(shì)。依據(jù)該特征可以進(jìn)行臺(tái)階特征線提取。

圖1 露天礦臺(tái)階組成示意Fig.1 Schematic of the composition of the open-pit steps

如圖2所示,假設(shè)點(diǎn)P為臺(tái)階線中某一點(diǎn),P1、P2,…,PN是以點(diǎn)P為中心、r為半徑的鄰近點(diǎn)。記鄰近點(diǎn)數(shù)量為N。理論上,鄰近點(diǎn)P1、P2,…,PN中高程(z分量)與中心點(diǎn)P高程相等的點(diǎn)的數(shù)量約占鄰近點(diǎn)數(shù)量N的50%。但是平盤地形起伏以及臺(tái)階坡面隆起凹陷將會(huì)造成以下影響:① 位于平盤上的鄰近點(diǎn)高程值并非處處相等。位于平盤上的鄰近點(diǎn)P1、P2,…,Pn(圖2中藍(lán)色點(diǎn))與中心點(diǎn)P的高程(z分量)之差小于某一高程差閾值Th,即 zPi-zP≤Th(zPi、zP分別表示鄰近點(diǎn)Pi與中心點(diǎn)P的高程,i=1,2,…,n,n＜N),同時(shí)記滿足該條件的點(diǎn)(平盤上的鄰近點(diǎn))的數(shù)量為n(n＜N)。② 平盤上鄰近點(diǎn)的數(shù)量n與總鄰近點(diǎn)數(shù)量N的比值在50%左右波動(dòng),即n/N∈ ( 0.5-l,0.5+l ),其中l(wèi)為波動(dòng)幅度。

圖2 臺(tái)階線特征點(diǎn)示意Fig.2 Schematic of the characteristic points of step line

確定出高程差閾值Th以及波動(dòng)幅度l后,即可對(duì)臺(tái)階線進(jìn)行提取。如圖3所示,紅色點(diǎn)為中心點(diǎn),r為搜索半徑,中心點(diǎn)上下兩條點(diǎn)劃線表示到中心點(diǎn)高程差在Th范圍的分界線,黃色點(diǎn)表示近鄰點(diǎn)云中到中心點(diǎn)高程差小于Th的離散點(diǎn)。由圖3可知:對(duì)于平坦區(qū)域(位置一)和陡峭區(qū)域(位置三)選用同樣的閾值,會(huì)將位置一處的中心點(diǎn)誤提取為特征點(diǎn)。

圖3 臺(tái)階線特征點(diǎn)提取示意Fig.3 Schematic of the feature point extraction of step line

為此提出了一種基于自適應(yīng)鄰域內(nèi)高程變化的特征點(diǎn)提取方法,依據(jù)前文1.1節(jié)中判定出的局部點(diǎn)云的特征維度選用不同的閾值。基本步驟如下:

(1)選定某一點(diǎn)P作為中心點(diǎn),采用K近鄰搜索(K Nearest Neighbor Search)法搜尋到中心點(diǎn)距離在半徑r(本研究中半徑r取2 m)以內(nèi)的鄰近點(diǎn)Pi(i=1,2,…,N,N為搜索到的鄰近點(diǎn)數(shù)量)。

(2)依據(jù)式(1)判斷中心點(diǎn)的鄰近點(diǎn)云的局部鄰域維度特性,按照中心點(diǎn)局部鄰域特性選擇不同的高程差閾值Th和波動(dòng)幅度l。

(3)遍歷點(diǎn)P的鄰近點(diǎn),比較鄰近點(diǎn)云中各離散點(diǎn)Pi的z分量與中心點(diǎn)P的z分量之差,記滿足zPi-zP≤Th的點(diǎn)的數(shù)量為n;計(jì)算n與N的比值,判定該中心點(diǎn)是否為臺(tái)階線特征點(diǎn)。

(4)以點(diǎn)云中各個(gè)離散點(diǎn)為中心點(diǎn),循環(huán)重復(fù)步驟(1)、(2)、(3),直至執(zhí)行完點(diǎn)云中所有離散點(diǎn)。

1.3 臺(tái)階線平滑

由于臺(tái)階平盤局部區(qū)域存在凹陷和突起現(xiàn)象,使得提取到的臺(tái)階線特征點(diǎn)中存在誤提取現(xiàn)象。在特征點(diǎn)平滑連接之前需要剔除誤提取的特征點(diǎn)。首先采用歐氏距離分割算法[12]對(duì)點(diǎn)云進(jìn)行聚類,其次在完成聚類分割后計(jì)算各個(gè)點(diǎn)云簇中離散點(diǎn)的數(shù)量Ni,同時(shí)設(shè)定一個(gè)數(shù)量閾值Nthresh。若Ni＜Nthresh,則對(duì)該類進(jìn)行剔除;反之,予以保留。本研究中數(shù)量閾值Nthresh經(jīng)試驗(yàn)驗(yàn)證最終設(shè)為200。在試驗(yàn)中發(fā)現(xiàn)經(jīng)剔除后的特征點(diǎn)中不僅包含有臺(tái)階線特征點(diǎn),還含有礦區(qū)道路邊線點(diǎn)云。

由于礦區(qū)車輛運(yùn)輸?shù)葘?shí)際需要,進(jìn)出采坑的道路均具有一定的坡度;而露天礦平盤是一較為規(guī)則的水平面,地形變化較小。采用點(diǎn)云的高程值作為區(qū)分道路邊線與臺(tái)階線特征點(diǎn)的依據(jù)。對(duì)于經(jīng)歐氏距離分割得到的某一簇道路邊線特征點(diǎn)Ω1而言,其高程值隨坡度不斷變化,該簇點(diǎn)云的高程值離散程度偏大,極差偏大;對(duì)于經(jīng)歐氏距離分割得到的某一簇臺(tái)階線特征點(diǎn)Ω2而言,該簇點(diǎn)云的高程值離散程度較小,高程極差較小。研究中手動(dòng)提取出了多簇道路邊線特征點(diǎn)與臺(tái)階線特征點(diǎn),經(jīng)計(jì)算得到如下規(guī)律:當(dāng)點(diǎn)云簇中各點(diǎn)高程值的方差小于0.75且極差值小于2 m時(shí),判定該點(diǎn)云簇為臺(tái)階線特征點(diǎn);反之,為道路邊線特征點(diǎn)。

利用上述規(guī)律對(duì)提取到的特征點(diǎn)進(jìn)行進(jìn)一步篩選,將篩選后的特征點(diǎn)作為最終提取的露天礦臺(tái)階線特征點(diǎn)。

臺(tái)階線是由提取到的無數(shù)個(gè)特征點(diǎn)按序連接而成的。但是上述方法提取到的特征點(diǎn)是無序的,因此需要對(duì)特征點(diǎn)進(jìn)行排序[13]。本研究主要依據(jù)離散點(diǎn)間距與偏轉(zhuǎn)角進(jìn)行點(diǎn)云排序。

如圖4所示,以P1點(diǎn)作為初始中心點(diǎn),依據(jù)最近點(diǎn)優(yōu)先原則,P2點(diǎn)為P1點(diǎn)的連接點(diǎn);以P2點(diǎn)作為中心點(diǎn),P3點(diǎn)為P2點(diǎn)的連接點(diǎn);以P3點(diǎn)作為中心點(diǎn)時(shí),由于P3到P4的距離與P3到P5的距離相等,即P3的最近點(diǎn)為P4和P5點(diǎn),此時(shí)最近點(diǎn)優(yōu)先原則失效,改用偏轉(zhuǎn)角來確定P3的連接點(diǎn)。從圖4中可以看出,P1P2和 P2P3之間的夾角 θ1為順時(shí)針轉(zhuǎn)角;P2P3和P3P4之間的偏轉(zhuǎn)角為θ2,且θ2為順時(shí)針轉(zhuǎn)角;P2P3和P3P5之間的偏轉(zhuǎn)角為θ3,且θ3為逆時(shí)針轉(zhuǎn)角。由圖4可進(jìn)一步看出,θ1與θ2轉(zhuǎn)向一致,大小相近;而θ1與θ3轉(zhuǎn)向相反,大小相差較大,故將P4點(diǎn)作為P3點(diǎn)的連接點(diǎn)。

圖4 離散點(diǎn)排序規(guī)則[13]Fig.4 Sorting rules for discrete points

經(jīng)排序得到有序點(diǎn)云后,采用三次B樣條函數(shù)對(duì)特征點(diǎn)進(jìn)行平滑連接處理[14]。三次B樣條曲線的方程為

式中,Fi,k(t)為基函數(shù),即:

因此,每4個(gè)離散點(diǎn)就可擬合一段曲線,例如P0、P1、P2、P3點(diǎn)可以擬合一段光滑曲線,P1、P2、P3、P4點(diǎn)可以擬合下一段,相鄰兩段曲線是平滑過渡的,以此類推,N個(gè)點(diǎn)可以擬合出(N-3)段平滑相接的曲線。

2 試驗(yàn)結(jié)果與分析

2.1 數(shù)據(jù)獲取及預(yù)處理

本研究使用的數(shù)據(jù)是由大疆Phantom 4 RTK無人機(jī)采集的鞍山齊大山鐵礦以及鞍千礦的影像數(shù)據(jù)。在Smart 3D軟件中對(duì)獲取到的影像進(jìn)行空中三角測(cè)量、多視影像密集匹配、紋理映射等操作得到研究區(qū)的密集匹配點(diǎn)云模型[15],如圖5所示。受到無人機(jī)拍照角度的影響,所構(gòu)建的密集匹配點(diǎn)云中各處的點(diǎn)云密度存在差異性。利用CloudCompare軟件計(jì)算出構(gòu)建的點(diǎn)云平均點(diǎn)密度約為48個(gè)/m2,相鄰點(diǎn)間最大高程差為0.048 9 m。

圖5 研究區(qū)三維點(diǎn)云模型Fig.5 3D point cloud model of the study area

本研究使用C++語言與PCL點(diǎn)云庫編寫算法程序,計(jì)算機(jī)配置參數(shù)見表1。

表1 試驗(yàn)所用的計(jì)算機(jī)配置Table 1 Computer configuration for the test

2.2 參數(shù)確定

由前文1.2節(jié)可知,本研究方法中影響露天礦臺(tái)階線特征點(diǎn)提取精度的參數(shù)有高程差閾值Th和波動(dòng)幅度l,且在不同地形條件下這兩個(gè)參數(shù)的取值應(yīng)是不同的。為了確定合理的參數(shù)取值,本研究基于鞍山齊大山鐵礦數(shù)據(jù)集對(duì)參數(shù)合理取值進(jìn)行試驗(yàn)分析。

2.2.1 高程差閾值Th選取

在保持區(qū)間波動(dòng)幅度l恒定(波動(dòng)幅度l=0.05)的情況下,確定高程差閾值Th的取值。由于所構(gòu)建的密集匹配點(diǎn)云中相鄰點(diǎn)間最大高程差為0.048 9 m,對(duì)該值進(jìn)行取整,以0.05 m為間隔選取Th,分別取Th=0.05、0.10、0.15、0.20、0.25、0.30 m,不同高程差閾值下的露天礦臺(tái)階線特征點(diǎn)提取結(jié)果見圖6。由圖6可知:當(dāng)Th取0.05 m時(shí),研究區(qū)內(nèi)地形平坦區(qū)域存在著大量誤提取的點(diǎn)云;隨著Th的增大,平坦區(qū)域內(nèi)誤提取的點(diǎn)云逐漸減少。當(dāng)Th取0.10 m時(shí),地形平坦區(qū)域中誤提取的點(diǎn)較圖5(a)中明顯減少;當(dāng)Th取0.30 m時(shí),由圖6(e)可見平坦區(qū)域內(nèi)幾乎沒有誤提取點(diǎn)。但是在Th增大的過程中,臺(tái)階上(下)部平盤中一些突起區(qū)域中逐漸出現(xiàn)了被提取出的點(diǎn)云,然而這些點(diǎn)云并非臺(tái)階線特征點(diǎn),如圖6(d)、圖6(e)、圖6(f)中橢圓標(biāo)注的位置。從圖6(c)中可以看出,臺(tái)階坡底線特征點(diǎn)提取結(jié)果優(yōu)于Th=0.10 m和0.20 m的提取結(jié)果。

為了保證臺(tái)階線特征點(diǎn)的提取精度,依據(jù)前文1.1節(jié)中鄰近點(diǎn)協(xié)方差矩陣的特征值進(jìn)行了高程差閾值Th選取,地勢(shì)平坦區(qū)域Th取0.30 m,地形復(fù)雜區(qū)域Th取0.15 m。

2.2.2 波動(dòng)幅度l選取

利用圖像構(gòu)建出的密集匹配點(diǎn)云中相鄰點(diǎn)間最大高程差為0.048 9m。隨機(jī)選取了數(shù)據(jù)集中5個(gè)點(diǎn)作為中心點(diǎn),同時(shí)分別搜索到與這5個(gè)中心點(diǎn)的距離為2m的鄰近點(diǎn),記錄鄰近點(diǎn)與中心點(diǎn)高程值之差小于(0.048 9t)m的離散點(diǎn)數(shù)量T,其中t=1,2,…,k且滿足 0.048 9k≤2。然后,計(jì)算高程差每增加0.048 9T m的增長(zhǎng)比例。

以選取的5個(gè)中心點(diǎn)為研究對(duì)象,試驗(yàn)發(fā)現(xiàn)高程差每增加0.048 9 m,滿足條件的離散點(diǎn)數(shù)量T平均增長(zhǎng)比例分別為 2.06%、3.64%、3.32%、2.51%、2.93%。為了保證提取的精度以及計(jì)算效率,本研究中將2%作為波動(dòng)幅度l的調(diào)整幅度,即l=0.02k,且滿足0.02k≤0.5。圖7為部分波動(dòng)幅度l取值的特征點(diǎn)提取結(jié)果。

圖7 不同波動(dòng)幅度l取值下的特征點(diǎn)提取結(jié)果Fig.7 Feature point extraction results under different fluctuation amplitude l values

由圖7可知:隨著波動(dòng)幅度l的增大,提取到的特征點(diǎn)不斷增多,特征點(diǎn)逐步從臺(tái)階坡線位置向臺(tái)階坡面上移動(dòng),開始遠(yuǎn)離臺(tái)階線所在的實(shí)際位置。當(dāng)波動(dòng)幅度l=0.12時(shí),提取到的非臺(tái)階線特征點(diǎn)增多,如圖7(d)所示;當(dāng)波動(dòng)幅度l=0.06時(shí),提取到的特征點(diǎn)存在斷續(xù)現(xiàn)象,漏提取的特征點(diǎn)較多。當(dāng)波動(dòng)幅度l=0.10時(shí)提取效果較為理想。

綜上所述,高程差閾值Th及波動(dòng)幅度l的取值如下:當(dāng)離散點(diǎn)對(duì)應(yīng)的鄰域點(diǎn)云協(xié)方差矩陣M的特征值 λ0、λ1、λ2存在 λ0＞＞λ1、λ2＞0 或 λ0、λ1＞＞λ2＞0關(guān)系,即局部鄰域呈線性或平面特征分布時(shí),高程差閾值Th取0.30m,波動(dòng)幅度 l=0.10;當(dāng) λ0、λ1、λ2存在λ0≈λ1≈λ2關(guān)系,即局部鄰域呈三維曲面特征分布時(shí),高程差閾值Th取0.15 m,波動(dòng)幅度l=0.10。

2.3 臺(tái)階線提取結(jié)果

在確定出本研究算法中的兩個(gè)重要參數(shù)后,對(duì)研究區(qū)進(jìn)行了臺(tái)階線特征點(diǎn)提取,結(jié)果如圖8所示。圖8(a)是利用統(tǒng)一參數(shù)標(biāo)準(zhǔn)提取的特征點(diǎn),該方法中未考慮地形起伏的影響,整個(gè)研究區(qū)內(nèi)的參數(shù)保持一致,從提取結(jié)果中可以看出存在著較多誤提取的離散點(diǎn)。圖8(b)為自適應(yīng)選用參數(shù)提取到的特征點(diǎn),從中可以看出平坦位置或多或少還存在著誤提取的特征點(diǎn),但是較圖8(a)相比誤提取的點(diǎn)數(shù)量明顯減少,提取效果得到了明顯改善。

圖8 特征點(diǎn)提取結(jié)果對(duì)比Fig.8 Comparison of the feature point extraction results

完成特征點(diǎn)提取后根據(jù)道路邊線特征對(duì)特征點(diǎn)進(jìn)行篩選,剔除道路邊線點(diǎn)。最后采用三次B樣條曲線對(duì)篩選的特征點(diǎn)進(jìn)行曲線擬合得到臺(tái)階線。圖9為最終提取到的兩組數(shù)據(jù)集露天礦臺(tái)階線,從圖中可以看出提取結(jié)果較好,可有效區(qū)分礦區(qū)的臺(tái)階線與道路邊線。

圖9 露天礦臺(tái)階線提取結(jié)果Fig.9 Extraction results of step line in open-pit mines

在提取效率方面,將該方法與傳統(tǒng)的GPS點(diǎn)測(cè)量方法進(jìn)行對(duì)比。對(duì)于1 km2作業(yè)區(qū)來說,使用一臺(tái)GNSS(Global Navigation Satellite System)接收機(jī)每5 m采集一個(gè)碎部點(diǎn)的方式進(jìn)行GPS點(diǎn)測(cè)量,其外業(yè)采集時(shí)間約15 h,內(nèi)業(yè)數(shù)據(jù)處理耗時(shí)約2.5 h;利用本研究方法進(jìn)行臺(tái)階線提取時(shí),飛馬D2000無人機(jī)航飛時(shí)間約22 min,內(nèi)業(yè)處理時(shí)間約4.5 h。從時(shí)間上來看,本研究方法減輕了外業(yè)工作量,同時(shí)極大地提高了工作效率。

3 結(jié) 論

(1)基于無人機(jī)密集匹配點(diǎn)云數(shù)據(jù),提出了自適應(yīng)鄰域內(nèi)高程閾值的露天礦不規(guī)則臺(tái)階線提取方法。試驗(yàn)結(jié)果表明:相比采用統(tǒng)一參數(shù)提取到的特征點(diǎn),該方法能有效減少誤提取的特征點(diǎn)數(shù)量,并可區(qū)分出臺(tái)階線特征點(diǎn)與道路邊線特征點(diǎn),較好地提取出露天礦中不規(guī)則的臺(tái)階線。與傳統(tǒng)臺(tái)階線提取方法相比,在減輕外業(yè)工作量的同時(shí),顯著提升了臺(tái)階線提取效率。

(2)目前該方法在部分區(qū)域由于閾值選取的偏差導(dǎo)致提取到的臺(tái)階線出現(xiàn)了錯(cuò)位現(xiàn)象,后期需要進(jìn)一步進(jìn)行方法優(yōu)化。