密集匹配點云下的礦石塊度提取

崔欣， 王井利*， 吳冬， 趙鑫

(1. 沈陽建筑大學交通與測繪工程學院， 沈陽 110168； 2. 中鐵十九局集團礦業投資有限公司， 北京 100161)

目前，對礦石塊度測定方法的研究大多分直接測量法與間接測量法。直接測量方法包括篩分法、二次爆破統計法、爆堆直接測量法等。間接測量方法包括相關數據測量法、視覺分析法、攝影測量法[1]。文獻[2-3]使用三維高分辨率激光掃描技術實現了爆破產生的礦石的全尺度量化。程雄等[4]采用數字圖像處理技術實現對礦石塊度參數的在線測量，并確定礦石塊度參數的分布。張建立等[5]基于工業攝像機獲取到的傳送帶上的礦石顆粒圖像，利用礦石顆粒像素標定以及像素面積檢測的方法研究了礦石的粒度分布。Ma等[6]通過對礦石破碎過程的分析，利用計算機視覺和圖像處理技術，實現了礦石粒度的在線檢測和粒度分布的計算。同樣的，文獻[7-9]也采用了圖像處理、圖像識別等圖像分析技術對礦石的塊度進行測定與統計。此類方法對獲取的礦石影像要求較高，礦石塊度提取精度受圖像分辨率的限制。Maron[10]提出了一種通過安裝新的粒度跟蹤(particle size tracking，PST)測量技術來估算礦石粒徑的方法，實現了在線粒度測量，由于數據量較大，算法處理實時性較差。吳開興等[11]使用水洗篩分法研究了單顆粒粒度分布特征，該算法不適用于大量礦石參數測定。Krauze等[12]基于震動信號處理方法開發一種基于振動的松散固體粒度和流動在線測量方法，實現了顆粒粒徑和流速的估算。在利用圖像識別塊度領域，Noy等[13]引入礦石堆的立體相機成像為基于圖像的測量算法提供了改進的自動粒子分割；孫厚廣等[14]對爆破后的礦山形態進行表面全景垂直拍照，通過軟件將照片劃分規則格網，通過對規則格網的分析，來評價爆破的塊度問題。此類基于圖像處理技術的塊度測定方法對相機分辨率要求較高。喻戰江等[15]利用三維激光掃描技術獲取爆破碎塊的幾何尺寸，實現了塊度分級，但并未真正意義上達到塊度識別的目的，方法精度有待提高。謝博等[16]利用3D激光掃描儀獲得巖塊表面模型，在內業處理中，通過不同的超體素聚類的閾值設置和權重因子分配，得到最佳的爆堆點云超體素聚類效果，之后利用區域生長法的思想進行分割達到塊度識別的目的。上述方法多采用圖像識別或者三維激光掃描技術進行礦石塊度測定，其主要問題在于礦石塊度精度低、提取效率差、僅適用于少量礦石塊度測定。為此，現采用無人機攝影測量密集匹配點云進行礦石塊度估算，旨在提高礦石塊度估算的精度與效率。

1 數據預處理

采用大疆Phantom 4 RTK無人機獲取了新疆磁海鐵礦某礦石堆的影像數據，利用ContextCapture Master對影像進行空中三角測量以及密集匹配點云模型的構建，構建出的三維模型如圖1所示。在完成三維點云模型構建之后針對礦石點云與地面點云相互粘黏問題，采用布料模擬濾波算法對點云進行濾波處理，實現地面點與非地面點的分離。經實驗發現在該算法中當布料分辨率為0.5，分類閾值為0.2，最大迭代次數取500時濾波效果最為理想。為了便于后期能夠快速地提取出每個礦石堆的礦石塊度，在分離出的礦石堆點云中采用基于歐氏距離的點云分割算法實現單個礦石堆點云的有效分割，最終經過預處理后的礦石堆分割結果如圖2所示。

圖1 研究區三維點云模型Fig.1 Three-dimensional point cloud model of the study area

圖2 礦石堆分割結果Fig.2 Results of ore pile segmentation

2 礦石堆坡度估算

坡度(slope)是地表單元陡緩的程度，通常把坡面的垂直高度h和水平方向的距離d的比叫作坡度，θ為該坡度下對應的坡度角，如圖3所示。其數學表達式為

(1)

圖3 地形坡度Fig.3 Topographic slope

依據坡度的定義式可知，尋找出礦石堆某一坡面處礦石點的最高點以及最低點即可完成對該坡面的坡度估計。對于同一礦石堆而言，最高點是唯一的，即礦石堆中高程最大值所對應的點。由于三維點云中點的三維坐標是通過(x，y，z)進行表示的，在點云處理過程中默認將三維坐標中的z分量作為高程值。所以本文研究在進行最高點確定時，首先遍歷所有坐標點的z值；其次比較了所有坐標點對應的z值，將z值最大的坐標點(xi，yi，zmax)作為礦石堆的最高點P，其中xi、yi表示點云中z值最大點對應的x與y軸分量。

礦石堆不同坡面的最低點并不是唯一確定的，最低點的選取與坡面有關。在本節中最低點是指在計算礦石堆某一坡向處礦石堆坡度時，該坡面的坡腳所對應的高程最低點。

為確定礦石堆不同坡面所對應的最低點，筆者以礦石堆點云的最高點為起始點，基于k近鄰搜索以及坐標方位角進行最低點搜尋，其整個實現步驟如下。

步驟1利用k近鄰搜索法搜尋距離最高點P的距離在[r，r+Δr]的所有離散點Ω。

步驟2將獲取到的離散點Ω與最高點P均投影至XOY平面上得到點P′和Ω′，在平面中計算環形點云Ω′中每一點與P′點連線的坐標方位角。

步驟3選取步驟2中計算的方位角與0°0′0″最為接近的點，以該點的投影前的坐標點作為下一搜索的中心點P1，重復步驟1～步驟3繼續進行最低點搜索，當同一方位角搜索出來的前后兩坐標點Pn、Pn+1的高程值(即z值)之差小于設定的閾值ε，即|Pn+1z-Pnz|<ε時，停止搜索，點Pn即為礦石堆在0°0′0″方向處的最低點。

步驟4同理選取步驟2中計算的方位角與kΔα(k=1，2，…且kΔα≤360°，Δα取值越小其后期分割結果越理想，本文中Δα取6°)最為接近的點，類似于尋找在0°0′0″方向處的最低點尋找出各個方位處礦石堆的最低點。

整個最低點的搜索過程剖面示意圖與俯視圖分別如圖4和圖5所示。

圖4 最低點搜尋剖面示意圖Fig.4 Schematic diagram of the lowest point search profile

圖5 最低點搜尋俯視圖Fig.5 The top view of the lowest point search

在完成礦石堆點云最高點以及最低點提取后，根據兩最值點三維坐標計算出礦石堆坡度i。計算公式為

(2)

式(2)中：xi、yi、zmax表示點云中最高點對應的x、y、z軸分量；xj、yj、zmin表示搜尋到的最低點對應的x、y、z軸分量。研究中將坡度計算結果與人工測量的坡度結果進行對比，分析坡度計算的誤差，其誤差分布直方圖如圖6所示。從圖6中可以看出計算得到大多數坡度的絕對誤差在0.1以內，且測定誤差呈現正態分布。

圖6 礦石堆坡度計算誤差直方圖Fig.6 Histogram of calculation error of ore pile slope

3 基于坡度信息的礦石點云分割

通過對礦石堆坡面進行逐層切割尋找出相鄰礦石間的交界點，以此獲得礦石點云邊界線，最后依據得到的邊界線完成對礦石點云的分割。

3.1 基于坡度信息的礦石堆點云切割

礦石堆點云的切割平面主要由礦石堆的坡度i以及最高點與最低點的三維坐標所決定。其基本原理如圖7所示。

利用求得礦石堆點云最高點P、最低點Q以及坡度i確定一個過最高點P與最低點Q同時與地面夾角為坡度角α的平面β0。將平面β0作為初始礦石堆點云裁剪的平面。其中平面β0的表達式為

(3)

式(3)中：xP、yP、zP分別為最高點P對應的三維坐標x、y、z分量；xQ、yQ、zQ分別為最高點Q對應的三維坐標x、y、z分量。

圖7 裁剪平面確定示意圖Fig.7 Schematic diagram of determining the clipping plane

但是由于礦石堆的坡面上會存在著多個凸起或凹陷的礦石，利用單個裁剪平面β0進行裁剪無法保證坡面上所有的礦石點云都被裁剪下來。另一方面，單一平面無法保證礦石點云與該平面相交所截出的平面為該礦石最大的截面，從而影響礦石塊度計算的精度。針對上述問題，提出了多層切片的裁剪方式。通過沿著平面β0的法向量方向上下平移平面β0得到更多的裁剪平面。

平面β0上下移動的距離以及最終平移的最大距離并非無限制的。而平面β0最終平移的最大距離由坡面表面礦石隆起的最大高度以及礦石凹陷的最大深度所決定的。本文以該坡面上所有離散點到平面β0的距離作為平面β0最終平移的最大距離判定依據。其實現步驟如下。

步驟1計算坡面內所有離散點(xk，yk，zk)，k=1，2，…到平面β0的距離dis：

(4)

式(4)中：xP、yP、zP為最高點P對應的三維坐標x、y、z分量；xQ、yQ、zQ為最高點Q對應的三維坐標x、y、z分量；xk、yk、zk(k=1，2，…)為第k個離散點對應的三維坐標x、y、z分量。

步驟2判斷離散點與平面的位置關系，并對距離dis進行方向標定，如式(5)所示。其中“+”表示離散點位于平面之上，“-”表示離散點位于平面之下。

(5)

式(5)中：zk(k=1，2，…)為第k個離散點對應的三維坐標z分量；z′k(k=1，2，…)為第k個離散點的x、y分量對應于平面該點處的高程值。

步驟3分別求取平面上下采樣點到平面的最大距離max_dis+、max_dis-，將兩個最大距離作為平面β0沿著法向量方向上下平移的極限值max_l+、max_l-。

在已知平面β0沿著法向量方向上下平移的極限值max_l+、max_l-后，可計算得到最高點及最低點上下移動的極值max_h+、max_h-分別為

(6)

根據最高點及最低點上下移動的極值max_h+、max_h-可知，最高點高程值的變化范圍為(zP- max_h-，zP+max_h+)。

由此可以得出第k個裁剪平面經過點(xP，yP，zP- max_h-+kΔh)。故第k個裁剪平面的平面方程為

[z-(zP-max_h-+kΔh)]=0

(7)

式(7)中：xP、yP、zP為最高點P對應的三維坐標x、y、z分量；xQ、yQ、zQ為最高點Q對應的三維坐標x、y、z分量；k為第k個裁剪平面；Δh為初始平面每次平移量；max_h+為平面上移的最大距離。

點云裁剪平面確定完成后即可根據點云中離散點的空間坐標以及裁剪平面的空間位置進行點云裁剪。

3.2 礦石點云交界線粗提取

以估算得到的坡度i為基礎對礦石堆點云進行多層切割處理。對前后兩次不同裁剪平面裁剪得到的點云進行分割處理，同時計算分割得到的每個礦石點云的中心點，其結果如圖8所示。依據前后兩次計算得到的中心點的數量以及位置關系進行礦石邊界線粗提取。

在圖8中，若兩次裁剪后的礦石點云中中心點數量以及位置關系相等或相近，則判斷前后兩次分割點云均來自同一個礦石，如圖8(a)～圖8(c)所示3個不同裁剪平面對同一坡面點云進行裁剪，其裁剪后的礦石點云中心點均為4個且4個中心點的位置一一對應[圖8(e)～圖8(h)]，則認為這3次裁剪出來的點云中不含有邊界點云；若兩個裁剪出的礦石點云中心點數量不一致，同時中心點位置也無法一一對應，則認為兩次切割點云中含有交界點，如圖8(c)～圖8(d)所示，兩個裁剪平面裁剪后的礦石點云分別形成了4個和3個點云，且點云中心點存在不對應現象，如圖8(g)和圖8(h)，表明圖8(d)裁剪出的點云中包含兩個礦石點云交界點。

圖8 不同裁剪平面切割出的點云Fig.8 Point clouds cut by different cutting planes

為了確定出礦石的交界點，將8(g)～圖8(h)結果中中心點不對應的點云進行作差處理，其結果如圖9所示，圖9中紅色點云即為作差得到的結果Ω。

兩點云作差結果中不僅包括了兩礦石的交界點，還包括礦石非交界點。為了剔除非交界點，計算圖8(h)結果中的兩中心點P1、P2到點云Ω中每一點的距離，并按從小到大進行排序，分別選取距離P1、P2點距離在某一閾值(實驗中選取為前60%)的點Ω1、Ω2，如圖10所示，圖10中紅色點云表示距離P1點距離在前60%的點Ω1，藍色點云表示距離P2點距離在前60%的點Ω2。最后對Ω1、Ω2進行求交運算，求交結果即為粗提取得到的兩礦石交界點，圖10中橘黃色點即為交界點。

圖9 點云作差結果Fig.9 Point cloud difference result

圖10 粗提取得到的交界點Fig.10 Junction points obtained by rough extraction

3.3 礦石點云交界線精提取

在3.2節中，受閾值選取的影響，提取的交界點中仍存在著非交界點，為提高礦石的分割精度，對交界點進行了進一步的精提取。

3.3.1 灰度信息提取

受光照條件的影響，礦石點云中礦石交界線處點云照比其他點云較黑，故在研究中將點云的紅(R)、綠(G)、藍(B)三個顏色通道信息轉換成灰度信息，其灰度值Gray的計算公式為

Gray=0.299R+0.587G+0.114B

(8)

式(8)中：Gray為點云灰度值；R、G、B為點云的紅、綠、藍3個顏色通道分量。

將計算得到的點云灰度值進行可視化，圖11(a)為礦石原始點云，圖11(b)為點云灰度可視化結果。從圖11(b)可以看出兩礦石交界處的點云呈現藍色或暗藍色，故點云的灰度信息可以用于礦石分割的依據。

3.3.2 鄰域特征分析

通過對兩礦石交界處點云進行分析發現，交界點處的點云高程值小于其周邊的大多數點云高程。

圖11 點云灰度可視化前后對比Fig.11 Comparison before and after point cloud gray scale visualization

依據這一特征可以對提取出的礦石交界處點云進行驗證，具體方法如下。

首先分別以每個粗提取出的礦石點云邊界點為中心進行K鄰近搜索尋找鄰近點云；其次比較近鄰點云高程與中心點高程的大小，記錄高程大于中心的鄰近點數量；若大于中心的鄰近點數量占鄰近點總數量的比值超過某一閾值，則判定該點為礦石交界點，反之則為非交界點。

在完成礦石交界點提取后，根據礦石交界點的連線完成對單個礦石的分割。圖12(a)為最終提取出的礦石點云的交界點，圖12(b)為依據提取出的礦石交界點分割得到的單個礦石結果。從圖12(b)中可以看出分割結果良好，但也存在著少量欠分割和過分割現象。

4 塊度測定及分布規律研究

礦石塊度(lumpiness)是指采用爆破等手段所開采的礦石碎塊的幾何尺寸的大小，常用礦石碎塊兩端最大距離表示。在圖像處理領域，學者們多采用等效匹配法對礦石參數進行估算，該方法中最為常用的算法包括最小外接矩形以及最小擬合橢圓等。

將分割得到的礦石點云進行投影變換得到近似于礦石最大截面的點云，采用了基于形狀特征的礦石塊度參數測定方法。即根據礦石截面的形狀特征選用不同的圖形對礦石截面邊界點云進行等效匹配，使得計算出的截面周長與截面面積誤差最小。本文采用了圓度的概念對礦石截面的形狀進行描述。圓度(circ_rity)是指礦石最大截面接近理論圓的程度。當圓度circ_rity越接近1時，表示截面形狀越近似于圓，此時采用最小擬合橢圓對礦石的最大截面進行等效匹配；否則采用最小外接矩形進行等效匹配。

圖12 基于礦石交界線的點云分割Fig.12 Point cloud segmentation based on ore boundary

圓度的定義式為

(9)

式(9)中：S為礦石截面面積；C為礦石截面周長。

實驗中選取了一個礦石堆中的32個礦石進行分析，計算出了32個礦石利用最小擬合橢圓和最小外接矩形進行等效匹配時面積、周長誤差與礦石最大截面圓度的關系。其結果如圖13(a)和圖13(b)所示。

從圖13(a)可知，當圓度circ_rity≥0.55時，利用橢圓擬合計算得到的周長誤差普遍小于利用最小外接矩形計算得到的結果，且大部分誤差在10%以內。從圖13(b)中可以看出，當circ_rity≥0.55時，利用橢圓擬合計算得到的面積誤差普遍小于利用最小外接矩形計算得到的結果。綜上所述，在計算礦石參數前，首先進行圓度計算，當圓度circ_rity≥0.55時，利用橢圓擬合的方式進行等效匹配；否則

圖13 周長、面積誤差隨圓度變化Fig.13 Circumference and area error change with roundness

利用最小外接矩形的方式進行等效匹配。并將擬合橢圓的長半軸或者最小外接矩形的長度作為礦石塊度l：

(10)

式(10)中：C為礦石最大截面的周長；S為礦石最大截面的面積；circ_rity為礦石最大截面圓度。

將提取出的礦石最大截面的長度作為最終提取出的礦石塊度，并利用統計學方式對礦石塊度進行統計分析，其結果如表1所示。同時繪制了礦石長度分布直方圖，如圖14所示。

從表1以及圖14中可以看出，礦石塊度介于0.5～1.0 m的比例最大，約為41.63%。塊度在4 m以上的礦石數量較少，最大礦石的塊度為7.674 m。

表1 礦石塊度統計表Table 1 Statistical table of ore lumpiness

圖14 礦石塊度分布直方圖Fig.14 Histogram of ore lumpiness distribution

5 結論

針對礦石塊度測定困難、精度低等問題，提出了基于密集匹配點云的礦石塊度計算方法。通過對本文方法進行試驗驗證，得到以下結論。

(1)利用搜尋礦石堆最高點、最低點的方式計算出的礦石堆坡度具有一定的可靠性，其計算誤差在0.1以內，且滿足正態分布。

(2)利用礦石堆的坡度信息、點云灰度信息以及點云鄰域特征可以較好地實現單個礦石點云的有效分割，且分割結果滿足后期礦石塊度計算的要求。

(3)研究區內的礦石塊度較為均勻，礦石塊度大多介于0.5～1.0 m。