孔德明，曹 帥，沈 閱，周逸人

(1. 燕山大學(xué) 電氣工程學(xué)院， 河北 秦皇島 066004； 2. 河北燕大燕軟信息系統(tǒng)有限公司， 河北 秦皇島 066004)

1 引 言

在新型智能化散雜貨港口中，散雜貨料堆的快速、高效的無(wú)人化機(jī)械堆取是港口未來(lái)的發(fā)展方向，全天候、快速、便捷和穩(wěn)定地得到料堆的離散點(diǎn)云數(shù)據(jù)和插值構(gòu)建數(shù)字高程模型(digital elevation model, DEM)是實(shí)現(xiàn)堆取料機(jī)全自動(dòng)化和體積測(cè)量的關(guān)鍵。

在現(xiàn)有的技術(shù)手段下，由于料堆形狀過(guò)高、紋理特征貧乏和材質(zhì)反射率低等特性，如何快速、高精度獲取被測(cè)物體表面的三維信息是測(cè)量問(wèn)題的重點(diǎn)[1]。離散點(diǎn)云數(shù)據(jù)的獲取主要有激光測(cè)量[2～4]、計(jì)算機(jī)視覺(jué)測(cè)量[5～7]等。楊德山等[8]提出采用手提測(cè)量系統(tǒng)，通過(guò)圍繞堆體行走時(shí)激光掃描儀動(dòng)態(tài)測(cè)量堆體表面的三維數(shù)據(jù)；張望等[9]提出操作者僅需手提掃描設(shè)備繞行被測(cè)貨堆一周便可快速獲得細(xì)節(jié)特征明顯的三維模型的方法；呂小寧等[10]基于斷面的激光精細(xì)測(cè)量系統(tǒng)，對(duì)某地下儲(chǔ)油洞室群進(jìn)行了現(xiàn)場(chǎng)實(shí)驗(yàn)，得到儲(chǔ)油洞室的三維數(shù)據(jù)；趙其杰等[11]提出激光掃描傳感器在傳動(dòng)平臺(tái)驅(qū)動(dòng)下，沿著堆積物料一側(cè)運(yùn)動(dòng)并進(jìn)行大范圍掃描采集被測(cè)點(diǎn)的空間坐標(biāo)信息。但是在露天港口的工作環(huán)境下，基于激光和視覺(jué)的方法不能夠在霧雨雪天氣以及沙塵環(huán)境下正常工作，因此，新的研究方法應(yīng)能克服惡劣環(huán)境的影響，使得自動(dòng)化堆取料機(jī)系統(tǒng)隨時(shí)能夠正常獲取離散點(diǎn)云數(shù)據(jù)并且能夠達(dá)到一定精度。

77 GHz毫米波雷達(dá)具有頻帶寬、波長(zhǎng)短、體積小、功耗低和穿透力強(qiáng)等特點(diǎn)，相比于激光紅外探測(cè)，其穿透性強(qiáng)的特點(diǎn)可以保證雷達(dá)能夠在霧雨雪以及沙塵環(huán)境工作，受天氣影響較小[12]。隨著毫米波雷達(dá)技術(shù)的不斷發(fā)展，在大型散雜貨料堆的測(cè)量應(yīng)用場(chǎng)景中，毫米波雷達(dá)能夠像二維激光雷達(dá)一樣通過(guò)掃描獲取得到目標(biāo)的位置信息。

對(duì)在港口中獲取料堆DEM的研究已有許多報(bào)道，毫米波雷達(dá)點(diǎn)云數(shù)據(jù)的逐點(diǎn)插值是構(gòu)建DEM最有效途徑之一[13]，常見(jiàn)的插值算法有：普通Kriging (ordinary kriging，OK)插值算法[14,15]、反距離加權(quán)(inverse distance weighted，IDW)插值算法[16]、基于三角剖分的線(xiàn)性插值(triangulation with linear interpolation，TLI)算法和自然鄰域(natural neighbor，NN)插值算法[17]等，針對(duì)測(cè)量裝置獲取的毫米波雷達(dá)離散點(diǎn)云數(shù)據(jù)存在空洞、不均勻和散亂度高的問(wèn)題，利用現(xiàn)有的插值方式對(duì)毫米波雷達(dá)點(diǎn)云進(jìn)行插值，得到的可視化模型和DEM精度較低，故應(yīng)改進(jìn)現(xiàn)有逐點(diǎn)插值方法以獲取更高精度的料堆DEM。

目前，關(guān)于基于毫米波雷達(dá)獲取料堆表面幾何信息的研究報(bào)道較少，本文在文獻(xiàn)[8,9]利用激光掃描獲取料堆離散點(diǎn)云數(shù)據(jù)的基礎(chǔ)上，提出了一種利用77 GHz毫米波雷達(dá)、差分北斗和角編碼器等傳感器獲取料堆表面離散點(diǎn)云數(shù)據(jù)的方法，分析離散點(diǎn)云數(shù)據(jù)誤差來(lái)源，通過(guò)改進(jìn)Kriging插值算法提高插值精度，一方面是可以補(bǔ)足缺失的空洞，另一方面是通過(guò)對(duì)插值區(qū)域的控制得到均勻的數(shù)據(jù)，采用交叉驗(yàn)證的方式對(duì)比分析不同離散點(diǎn)云數(shù)據(jù)插值方法的插值精度。實(shí)驗(yàn)結(jié)果證明：本文提出的裝置和方法具有很好的環(huán)境適應(yīng)性、穩(wěn)定性、便捷性和低成本的特點(diǎn)，能夠滿(mǎn)足實(shí)際的項(xiàng)目需求。

2 獲取料堆離散點(diǎn)云數(shù)據(jù)

2.1 數(shù)據(jù)采集

2.1.1 基于77 GHz毫米波雷達(dá)數(shù)據(jù)采集裝置

以堆料機(jī)為例，毫米波雷達(dá)傳感器安裝在出料口前方，獲取料堆表面坐標(biāo)信息；雙天線(xiàn)差分北斗傳感器安裝在回轉(zhuǎn)平臺(tái)上，獲取回轉(zhuǎn)中心在堆場(chǎng)內(nèi)的位置數(shù)據(jù)；角編碼器獲取堆料機(jī)大臂的俯仰角和回轉(zhuǎn)角數(shù)據(jù)。77 GHz毫米波雷達(dá)是能夠在全天候場(chǎng)景下快速感知0～200 m范圍內(nèi)周邊環(huán)境內(nèi)目標(biāo)信息的長(zhǎng)距離雙波束雷達(dá)，遠(yuǎn)距離波束的掃描范圍為水平方向45°，近距離波束的的掃描范圍是60°，并且可以72 ms完成一次近距和遠(yuǎn)距的測(cè)量，毫米波雷達(dá)工作在集群模式，其掃描面垂直向下掃描料堆的橫截面，能夠返回測(cè)量點(diǎn)的橫向距離、縱向距離、雷達(dá)的反射截面積、虛警概率等信息，其掃描示意圖如圖1所示；采用精確到cm級(jí)的差分北斗定位傳感器；采用精確到0.15°的角編碼器。

圖1 毫米波雷達(dá)掃描示意圖Fig.1 Schematic diagram of millimeter wave radar scanning

2.1.2 堆場(chǎng)模型

在長(zhǎng)方形露天的自動(dòng)化散雜貨堆場(chǎng)，建立堆場(chǎng)全局坐標(biāo)系O0-X0Y0Z0，以水平向北為Y0軸，水平向東為X0，Z0軸垂直于長(zhǎng)方形堆場(chǎng)；定義堆料機(jī)在堆場(chǎng)軌道上的位置為堆料機(jī)坐標(biāo)系O1-X1Y1Z1；其中X1O1Z1平面隨著大臂的轉(zhuǎn)動(dòng)進(jìn)行轉(zhuǎn)動(dòng)；定義毫米波雷達(dá)坐標(biāo)系O-XYZ。料堆的堆場(chǎng)模型由圖2所示，θ是大臂的俯仰角，L是毫米波雷達(dá)到回轉(zhuǎn)中心的長(zhǎng)度，H是回轉(zhuǎn)中心到堆料機(jī)坐標(biāo)系原點(diǎn)的高度，φ是大臂的回轉(zhuǎn)角。

圖2 堆場(chǎng)料堆模型示意圖Fig.2 Schematic diagram of stockyard model

將毫米波雷達(dá)獲取到的數(shù)據(jù)P(x,y,z)、北斗定位數(shù)據(jù)和角度傳感器數(shù)據(jù)通過(guò)公式(1)解算到堆場(chǎng)坐標(biāo)系O0-X0Y0Z0，得到料堆表面的三維數(shù)據(jù)P0(x0,y0,z0)。

(1)

2.1.3 測(cè)量過(guò)程

測(cè)量過(guò)程：通過(guò)PLC控制堆料機(jī)的大臂俯仰到一定角度，以不同的回轉(zhuǎn)角往復(fù)低速走行采集料堆表面的幾何信息，同時(shí)存儲(chǔ)北斗傳感器和角編碼器的數(shù)據(jù)，將多傳感器數(shù)據(jù)發(fā)送到數(shù)據(jù)處理服務(wù)器，通過(guò)對(duì)多傳感器的數(shù)據(jù)解算和坐標(biāo)系轉(zhuǎn)換得到料堆多組三維數(shù)據(jù)，對(duì)多組點(diǎn)云數(shù)據(jù)進(jìn)行濾波和精簡(jiǎn)，之后進(jìn)行配準(zhǔn)得到料堆點(diǎn)云數(shù)據(jù)，采用點(diǎn)云數(shù)據(jù)投影區(qū)域的逐點(diǎn)插值算法獲取均勻的點(diǎn)云數(shù)據(jù)，料堆DEM獲取過(guò)程如圖3所示。

圖3 料堆DEM模型獲取過(guò)程Fig.3 Stock model acquisition process diagram

3 Kriging插值算法及其改進(jìn)

3.1 普通Kriging插值算法

將采集到的離散點(diǎn)云數(shù)據(jù)投影到二維平面并求取凸包，根據(jù)情況對(duì)凸包內(nèi)區(qū)域劃分出網(wǎng)格點(diǎn)，生成查詢(xún)點(diǎn)，查找到查詢(xún)點(diǎn)周?chē)阎叱讨档狞c(diǎn)，通過(guò)普通Kriging插值算法得到查詢(xún)點(diǎn)的高程值。目前已有學(xué)者研究基于普通Kriging算法的點(diǎn)云數(shù)據(jù)的插值方法[18]，普通Kriging插值算法原理如下：

對(duì)于查詢(xún)點(diǎn)附近的n個(gè)臨近點(diǎn)t1,t2,t3,…,tn，其觀(guān)測(cè)值為w(t1),w(t2),w(t3),…,w(tn)，對(duì)于查詢(xún)點(diǎn)t0的估計(jì)值可以通過(guò)一個(gè)線(xiàn)性關(guān)系來(lái)估值：

(2)

式中λi為權(quán)值。

由于Kriging插值是一種無(wú)偏最優(yōu)化的插值，無(wú)偏性和最小化方差是確定λi的兩個(gè)條件如下：

(3)

當(dāng)區(qū)域的變量屬于二階平穩(wěn)假設(shè)和本證假設(shè)，由公式(3)的兩個(gè)條件可推導(dǎo)Kriging方程組為

(4)

半方差函數(shù)模型γ(h)通過(guò)式(5)進(jìn)行計(jì)算，

(5)

式中：N(h)是臨近點(diǎn)的點(diǎn)對(duì)數(shù)；h是滯后距；w(ti)、w(ti+h) 分別是ti和ti+h位置的觀(guān)測(cè)值。

常見(jiàn)的擬合半方差函數(shù)模型有球狀模型、高斯模型、冪函數(shù)模型、空穴函數(shù)模型等，可作為半方差函數(shù)的球狀函數(shù)模型：

(6)

式中：c0是塊金值；a是變程；c是拱高。

3.2 基于量子化鴿群優(yōu)化的Kriging插值算法

由于毫米波雷達(dá)點(diǎn)云數(shù)據(jù)散亂程度較大，造成Kriging插值算法擬合半方差函數(shù)時(shí)精度不高，所以采量子化鴿群算法[19,20]優(yōu)化擬合半方差函數(shù)?；绝澴尤簝?yōu)化算法包括地圖指南針?biāo)阕雍偷貥?biāo)算子兩部分。地圖和指南針?biāo)阕邮悄７绿?yáng)和地球磁場(chǎng)這兩種導(dǎo)航工具對(duì)鴿子的影響。鴿群中第j只鴿子在第s代的速度信息和位置信息的更新策略如下：

Vj(s+1)=e-RtVj(s)+ru[(Ggb(s)-G(s)]

(7)

Gj(s+1)=Gj(s)+Vj(s+1)

(8)

式中：R是羅盤(pán)算子，取值范圍是0到1；Vj(s)是鴿子的速度信息；Gj(s)是鴿子的位置信息；ru表示隨機(jī)數(shù)；Ggb(s)鴿群中所有鴿子的全局最優(yōu)的位置。

當(dāng)?shù)螖?shù)達(dá)到一定值時(shí)停止地圖指南針?biāo)阕拥墓ぷ?，進(jìn)入地標(biāo)算子的迭代過(guò)程中繼續(xù)工作，每一次迭代鴿子數(shù)量減半，位置更新策略為

Np(s)=Np(s-1)/2

(9)

(10)

Gi(s)=Gi(s-1)+ru(GC(s)-Gi(s-1))

(11)

式中：Np是當(dāng)前鴿子數(shù)量；F(Gi(s))是自適應(yīng)度函數(shù)；GC(s)是鴿子群的中心且被當(dāng)做地標(biāo)；ru是0到1之間的隨機(jī)數(shù)。

量子化鴿群優(yōu)化算法引入實(shí)數(shù)的量子表示和量子旋轉(zhuǎn)門(mén)。一個(gè)量子可以通過(guò)“0”和“1”態(tài)表示，量子位狀態(tài)被表示為

|ψi〉=αi|0〉+βi|1〉

(12)

(13)

(14)

式中：Δθ是旋轉(zhuǎn)角度；ε是約束參數(shù)，防止量子困在狀態(tài)“0”和“1”。

4 離散點(diǎn)云數(shù)據(jù)的誤差來(lái)源以及插值算法評(píng)價(jià)指標(biāo)

4.1 誤差來(lái)源

影響料堆表面離散點(diǎn)云數(shù)據(jù)準(zhǔn)確性的原因主要有4類(lèi)：第1類(lèi)是堆取料機(jī)在采集數(shù)據(jù)的過(guò)程中會(huì)發(fā)生輕度的震動(dòng)，造成掃描點(diǎn)位置偏差；第2類(lèi)是傳感器的測(cè)量誤差，主要有毫米波雷達(dá)的測(cè)距和測(cè)角誤差，角編碼器的測(cè)角誤差，差分北斗傳感器的定位誤差以及不同測(cè)量數(shù)據(jù)時(shí)間同步誤差等；第3類(lèi)是測(cè)量設(shè)備安裝的誤差，主要有毫米波雷達(dá)的安裝誤差角和位置偏移量等；第4類(lèi)是料堆的形狀不同可能會(huì)造成料堆表面數(shù)據(jù)缺失，造成插值誤差。其中，第3類(lèi)測(cè)量設(shè)備的安裝誤差可通過(guò)采集數(shù)據(jù)的直線(xiàn)特征提取并結(jié)合結(jié)構(gòu)參數(shù)列出約束方程求出初始狀態(tài)標(biāo)定值，再利用初始狀態(tài)值和測(cè)量過(guò)程中獲取的系統(tǒng)位置姿態(tài)變化量實(shí)現(xiàn)在線(xiàn)標(biāo)定和測(cè)量。

4.2 評(píng)價(jià)指標(biāo)

采用交叉驗(yàn)證的方式評(píng)價(jià)不同的逐點(diǎn)插值算法的精度以及料堆表面的還原度，對(duì)處理之后的離散點(diǎn)云數(shù)據(jù)采樣80%的點(diǎn)云數(shù)據(jù)作為內(nèi)插樣本數(shù)據(jù)，20%的點(diǎn)云數(shù)據(jù)作為實(shí)測(cè)數(shù)據(jù)來(lái)驗(yàn)證插值精度。誤差的大小通過(guò)高程值的差異進(jìn)行反應(yīng)，通過(guò)均方誤差(MSE)、均方根誤差(RMSE)和預(yù)測(cè)吻合度R2評(píng)價(jià)插值的精度，計(jì)算公式如下：

(15)

(16)

(17)

5 算法驗(yàn)證

5.1 點(diǎn)云數(shù)據(jù)統(tǒng)計(jì)信息

利用上述算法，在某煤炭港務(wù)有限公司的自動(dòng)化改造項(xiàng)目中進(jìn)行測(cè)試應(yīng)用，用于料堆表面點(diǎn)云數(shù)據(jù)獲取并通過(guò)插值算法得到料堆的DEM。對(duì)于80 m寬的堆場(chǎng)，堆料機(jī)的大臂俯仰角達(dá)到12.50°，回轉(zhuǎn)角以65°、55°、40°和30°全自動(dòng)走行，同時(shí)采集并存儲(chǔ)多個(gè)傳感器的數(shù)據(jù)，融合解算后得到離散點(diǎn)云數(shù)據(jù)。對(duì)數(shù)據(jù)處理之后共得到6組料堆離散點(diǎn)云數(shù)據(jù)，其每平方米點(diǎn)的個(gè)數(shù)在25個(gè)左右，由峰度可知料堆1和料堆2的高程值更接近正態(tài)分布，料堆3和料堆4稍微陡峭一些，由偏度可知6個(gè)料堆的高程值分布都右偏，料堆都有局部的狹長(zhǎng)空洞，且形狀和面積都不相同，統(tǒng)計(jì)信息包括離散點(diǎn)云數(shù)據(jù)中點(diǎn)的個(gè)數(shù)，離散點(diǎn)云數(shù)據(jù)的投影凸包的面積stt、最大高程值hmax、平均高程值hv、高程值的標(biāo)準(zhǔn)差hσ、高程值的峰度、高程值的偏度。料堆的離散點(diǎn)云數(shù)據(jù)詳細(xì)的統(tǒng)計(jì)信息見(jiàn)表1。

5.2 不同插值方法得到的DEM

將基于毫米波雷達(dá)多傳感器融合得到料堆1的離散點(diǎn)云數(shù)據(jù)經(jīng)過(guò)去噪、去除地面點(diǎn)、精簡(jiǎn)和配準(zhǔn)之后,得到具有狹長(zhǎng)空洞的離散點(diǎn)云數(shù)據(jù)見(jiàn)圖4中的原始點(diǎn)云數(shù)據(jù)。選擇0.1 m分辨率的格網(wǎng)尺寸，采用的改進(jìn)的Kriging插值算法(improve ordinary kriging，IOK)、OK、IDW、TLI和NN插值算法獲取DEM如圖4所示，圖4中的插值算法都能夠填補(bǔ)原始點(diǎn)云數(shù)據(jù)的漏洞得到完整的DEM，它們是獲取地形地勢(shì)、建筑物DEM的主要方式。

表1 料堆離散點(diǎn)云數(shù)據(jù)的描述參數(shù)統(tǒng)計(jì)Tab.1 Description parameter statistics of the discrete point cloud data of the stockpile

圖4 不同插值方法得到的DEMFig.4 Elevation model obtained by different interpolation methods

5.3 不同插值方法精度分析

對(duì)5種的離散點(diǎn)云數(shù)據(jù)插值算法進(jìn)行參數(shù)優(yōu)選，IOK初始以地圖和指南針?biāo)阕訛閷?dǎo)向的迭代次數(shù)設(shè)置為70，地標(biāo)算子為導(dǎo)向的迭代次數(shù)為30次，種群數(shù)量為80；IDW插值算法的冪次設(shè)置為2次，搜索25個(gè)臨近點(diǎn)；TLI插值算法以Delaunay三角化準(zhǔn)則構(gòu)造三角形為基礎(chǔ)進(jìn)行插值；NN插值算法以Delaunay三角化構(gòu)造Voronoi多邊形為基礎(chǔ)進(jìn)行插值。為了評(píng)價(jià)不同插值算法的插值精度，采樣料堆離散點(diǎn)云數(shù)據(jù)的20%當(dāng)作真值，得到不同料堆的不同插值算法的有效性交叉驗(yàn)證結(jié)果，6組不同料堆的均方誤差MSE和均方根誤差RMSE交叉驗(yàn)證結(jié)果如圖5、圖6所示。

圖5 均方誤差的結(jié)果Fig.5 Results of mean square error

圖6 均方根誤差的結(jié)果Fig.6 Results of root mean square error

由圖5和圖6所示的6組數(shù)據(jù)可知IOK算法RMSE低于0.37 m，MSE低于0.14 m, 預(yù)測(cè)吻合度R2達(dá)到0.984以上；OK插值算法獲取的6組數(shù)據(jù)的RMSE大于0.5 m，MSE大于0.32 m，R2低于0.977 7，其中4組數(shù)據(jù)R2低于0.96；IDW、NN、TLI的插值算法獲取的交叉驗(yàn)證結(jié)果較差，RMSE大于0.5 m，R2低于0.958 5。

由圖5和圖6所示的MSE和RMSE說(shuō)明，基于毫米波雷達(dá)獲取的料堆離散點(diǎn)云數(shù)據(jù)采用不同的逐點(diǎn)插值算法獲取料堆DEM的有效性依次是IOK插值算法、OK插值算法、IDW插值算法、NN插值算法、TLI插值算法，其中IDW、TLI、NN算法誤差浮動(dòng)較大，可知處理不同料堆的離散點(diǎn)云數(shù)據(jù)時(shí)IOK插值算法都能取得較好的效果，具有較好的穩(wěn)定性，IOK插值算法的RMSE穩(wěn)定在0.2 m到0.4 m之間，計(jì)算每一組數(shù)據(jù)的IOK插值算法和OK插值算法RMSE的相對(duì)偏差，可知IOK的RMSE相比OK算法降低了39.9%，改進(jìn)后的插值算法更加適合毫米波雷達(dá)點(diǎn)云數(shù)據(jù)獲取DEM。

6 結(jié) 論

在某煤炭港口的自動(dòng)化堆取料系統(tǒng)實(shí)際應(yīng)用毫米波雷達(dá)、差分北斗傳感器和角編碼器傳感器集成的技術(shù)，全自動(dòng)全天候獲取料堆的DEM，基于毫米波雷達(dá)取點(diǎn)云數(shù)據(jù)的方式可以彌補(bǔ)雨雪以及大量灰塵工況下激光和視覺(jué)方式測(cè)量的短板；插值算法的有效性依次是IOK插值算法、OK插值算法、IDW插值算法、NN插值算法、TLI插值算法，改進(jìn)后的Kriging插值算法提高了毫米波雷達(dá)點(diǎn)云數(shù)據(jù)的插值精度，其RMSE低于0.37 m，MSE低于0.14 m，RMSE相比OK插值算法降低了39.9%，R2達(dá)到0.984以上；得到的DEM能夠滿(mǎn)足堆取料機(jī)自動(dòng)化項(xiàng)目對(duì)精度的需求。