賈象陽，黃先鋒，牛文淵，張 帆，高云龍，楊 沖

武漢大學(xué)測繪遙感信息工程國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079

在測量工程中，回歸算法經(jīng)常被用于解決含有噪聲的數(shù)據(jù)模型參數(shù)估計(jì)問題，其中，經(jīng)典的RANSAC算法[1]憑借其實(shí)現(xiàn)簡單和穩(wěn)健等優(yōu)勢(shì)在圖像匹配[2-3]、點(diǎn)云定向[4-6]、三維模型修復(fù)[7]和機(jī)器人定位抓取[8]等領(lǐng)域有著廣泛應(yīng)用。

根據(jù)經(jīng)典RANSAC算法的假設(shè)-驗(yàn)證策略，所有RANSAC改進(jìn)算法針對(duì)3個(gè)方面進(jìn)行優(yōu)化，包括：①優(yōu)化假設(shè)模型[9-11]，②改進(jìn)驗(yàn)證策略[12]，③降低時(shí)間成本[13]。LO-RANSAC[9]、OR-RANSAC[14]和U-RANSAC[15]算法通過改進(jìn)采樣的方式優(yōu)化初始假設(shè)模型，提高模型參數(shù)估計(jì)精度的同時(shí)，也降低了時(shí)間成本。與之類似的算法還有NAPSAC[11]、PROSAC[10]、GOODSAC[16]及SCRAMSAC[17]等，他們都在內(nèi)群點(diǎn)比噪聲點(diǎn)更加緊湊的前提下對(duì)模型參數(shù)進(jìn)行估計(jì)。在文獻(xiàn)[12]基于內(nèi)群點(diǎn)服從無偏高斯分布以及噪聲服從均勻分布的假設(shè)提出MLESAC算法之后，MAPSAC[18]、RMLESAC[19]、AMLESAC[20]、u-MLESAC[21]和SASAC[22]等算法對(duì)內(nèi)群點(diǎn)比率(γ)和方差(σ2)的精度和穩(wěn)健性進(jìn)行了優(yōu)化，試驗(yàn)效果顯著?；谖墨I(xiàn)[12]的假設(shè)，文獻(xiàn)[23]提出一種基于NFA(number of false alarms)的對(duì)抗策略，保證了內(nèi)群點(diǎn)比例自適應(yīng)的同時(shí)，對(duì)采樣策略和模型評(píng)定準(zhǔn)則也進(jìn)行了改進(jìn),但是，仍然需要人工給定閾值(如：NFA閾值)。另外，有些學(xué)者對(duì)傳統(tǒng)RANSAC算法進(jìn)行優(yōu)化，不考慮數(shù)據(jù)分布情況，比如：文獻(xiàn)[24]通過窮舉一定范圍內(nèi)的誤差容忍閾值，并采用固定迭代次數(shù)的方式估計(jì)模型參數(shù)方差，最小方差對(duì)應(yīng)的誤差容忍閾值與模型參數(shù)作為最優(yōu)估計(jì)。

但是，現(xiàn)有的RANSAC改進(jìn)算法在尋求最優(yōu)回歸模型時(shí)都需要人工給定先驗(yàn)信息，這對(duì)于算法精度、效率、穩(wěn)健性及普適性都會(huì)造成不同程度的影響?；诖耍疚奶岢隽艘环N參數(shù)自適應(yīng)的RANSAC算法，不需要給定任何先驗(yàn)信息，在保證不降低迭代有效性的情況下，模型回歸精度要優(yōu)于現(xiàn)有算法，并分別通過仿真模擬和真實(shí)場景進(jìn)行了對(duì)比試驗(yàn)分析。

1 算法原理

基于文獻(xiàn)[11]的假設(shè)：相較于噪聲，內(nèi)群點(diǎn)在空間中分布更加緊湊(相比于噪聲，①內(nèi)群點(diǎn)平均密度或數(shù)量通常更大；②內(nèi)群點(diǎn)分布與給定模型更加吻合)。本文提出一種空間密度函數(shù)，用于描述點(diǎn)云空間分布與密度之間的關(guān)系，并自適應(yīng)地估計(jì)模型參數(shù),技術(shù)路線見圖1。

注: 灰色矩形框表示本文的改進(jìn)和優(yōu)化之處，整體為原始RANSAC的流程。

(1) 遍歷原始數(shù)據(jù)集，統(tǒng)計(jì)并計(jì)算每個(gè)點(diǎn)的鄰域半徑(詳見1.1節(jié)解釋)以及搜索半徑，選取落在搜索半徑內(nèi)的點(diǎn)個(gè)數(shù)作為采樣權(quán)重。

(2) 改變傳統(tǒng)RANSAC算法僅依賴內(nèi)群點(diǎn)數(shù)量判斷模型好壞的方式，根據(jù)內(nèi)群點(diǎn)的空間密度信息表征，聯(lián)合誤差與內(nèi)群點(diǎn)分布判斷模型的好壞。

(3) 利用空間密度函數(shù)計(jì)算其終止條件，其算法見算法1。

算法1 基于空間密度函數(shù)的自適應(yīng)RANSAC算法

變量：模型:M,最優(yōu)模型:best_M,內(nèi)群點(diǎn):inliers,最優(yōu)內(nèi)群點(diǎn):best_inliers;最優(yōu)誤差容忍閾值:best_α;本次誤差容忍閾值:α;上一次誤差容忍閾值:α′;散亂點(diǎn)云:spc;容忍失敗率:ε;迭代次數(shù):t;本次空間密度:ρ;上一次空間密度:ρ′;點(diǎn)到當(dāng)前M的符號(hào)距離:d;最大和最小符號(hào)距離:dmax,dmin;內(nèi)群點(diǎn)個(gè)數(shù):Nin;點(diǎn)云個(gè)數(shù)：N;點(diǎn)到當(dāng)前M的跨度:S;模型評(píng)價(jià)參數(shù):ω

初始化：ε=0.01，t=1，α′=+，ρ′=0.0

輸入：spc

輸出：best_M,inliers

1. 計(jì)算spc中每個(gè)點(diǎn)的權(quán)重因子(見1.1小節(jié))；

2. 估計(jì)最優(yōu)模型和內(nèi)群點(diǎn)集：

whilet>0

利用加權(quán)隨機(jī)采樣并估計(jì)初始假設(shè)模型M;

根據(jù)α′確定內(nèi)群點(diǎn)inliers;

ifNin/N<0.1

t=t-1;

continue;%重新初始化假設(shè)模型M

end if

利用最小二乘方法對(duì)inliers重新估計(jì)M;

計(jì)算dmax,dmin,S,Nin,N,α,ρ,其中：α=(max(|dmax|,|dmin|)+mean(|d|))/2.0;

ifρ≥ρ′

ift==1

t=N(N-1)(N-2);

else

end if

best_M=M;

best_inliers=inliers;

best_α=α′;

α′=α;

ρ′=ρ;

else

t=t-1;

end if

end while

1.1 加權(quán)采樣

傳統(tǒng)RANSAC算法中每次采樣都具有獨(dú)立隨機(jī)性，且所有樣本都被等概率采樣。然而，在噪聲多于內(nèi)群點(diǎn)時(shí)，內(nèi)群點(diǎn)被采樣的概率明顯會(huì)低于噪聲，這就會(huì)無法估計(jì)出良好的初始模型。內(nèi)群點(diǎn)平均點(diǎn)密度或內(nèi)群點(diǎn)數(shù)不會(huì)在很大程度上小于噪聲，反之，噪聲就會(huì)被認(rèn)定為內(nèi)群點(diǎn)，原始內(nèi)群點(diǎn)被看作是噪聲。本文首先通過統(tǒng)計(jì)每個(gè)目標(biāo)點(diǎn)pobj的k個(gè)最近鄰點(diǎn)(k=1用于所有試驗(yàn))，記為：ψ={p1,p2,…,pk}，并計(jì)算k個(gè)最近鄰點(diǎn)到pobj之間的歐幾里得距離，記為：η={η1,η2,…,ηk}，其中ηi=‖pi-pobj‖,i∈1,2,…,k，然后計(jì)算鄰域半徑r=median(η)，最后，遍歷整個(gè)數(shù)據(jù)集，統(tǒng)計(jì)Ω={r1,r2,…,rN}，N為點(diǎn)云總個(gè)數(shù)，計(jì)算搜索半徑R=mean(Ω)，并統(tǒng)計(jì)所有目標(biāo)點(diǎn)在搜索半徑R范圍內(nèi)的鄰域點(diǎn)個(gè)數(shù)，將其作為目標(biāo)點(diǎn)權(quán)重因子w={w1,w2,…,wN}，保證在每次迭代中估計(jì)初始模型參數(shù)時(shí)，能夠最大可能地在內(nèi)群點(diǎn)中選取最小樣本，降低時(shí)間代價(jià)。另外， 為了提高權(quán)重因子統(tǒng)計(jì)的效率，k-d樹被用于最近鄰搜索和范圍鄰域點(diǎn)個(gè)數(shù)的統(tǒng)計(jì)過程。

搜索半徑R通過影響權(quán)重因子w，使得采樣樣本有更大的概率來自于內(nèi)群點(diǎn)。但是，如果R過大，每個(gè)點(diǎn)的權(quán)重差異性就會(huì)相對(duì)減弱，如果R過小，對(duì)于誤差較小的噪聲作用不明顯。而對(duì)于不含噪聲的數(shù)據(jù)而言，最近鄰點(diǎn)數(shù)k對(duì)搜索半徑R的影響較小，即：R對(duì)k不敏感，k可以在較大的閾值范圍內(nèi)取值。為了減弱噪聲對(duì)搜索半徑R的影響，本文引入互為最近鄰方法抑制噪聲，其原理是最近鄰點(diǎn)落在目標(biāo)點(diǎn)搜索的范圍內(nèi)，同時(shí)目標(biāo)點(diǎn)也落在k個(gè)最近鄰點(diǎn)的鄰域范圍內(nèi)，如圖2所示，藍(lán)色實(shí)心圓圈不在目標(biāo)點(diǎn)的最近鄰范圍內(nèi)，而綠色實(shí)心圓圈落在了藍(lán)色實(shí)心圓圈的最近鄰范圍內(nèi)。

注：綠色實(shí)心圓圈為目標(biāo)點(diǎn)，藍(lán)色實(shí)心圓圈為噪聲點(diǎn)，實(shí)線指向目標(biāo)點(diǎn)鄰域點(diǎn)，虛線指向噪聲點(diǎn)鄰域點(diǎn)。

1.2 空間密度函數(shù)

當(dāng)場景尺度未知時(shí)，誤差容忍閾值無法通過人工給定的方式保證高精度估計(jì)模型參數(shù)，但是，一般情形下，內(nèi)群點(diǎn)比率與誤差容忍閾值往往存在相關(guān)性，其適用于不同場景模型。

以平面為例(下同)，假設(shè)數(shù)據(jù)集在容忍誤差內(nèi)沒有噪聲，其表現(xiàn)為更加緊湊，跨度S將會(huì)較小(圖3(a))，因此，在最優(yōu)平面(圖3的紅色平面)投影面積一定的情況下，每個(gè)點(diǎn)的占有平均體積較小。而在圖3(b)中，由于存在噪聲且噪聲往往分布較為分散，跨度S的增加使得每個(gè)點(diǎn)占有的平均體積明顯增加。另外，噪聲還有一個(gè)明顯的特征是相較于內(nèi)群點(diǎn)，其偏離最優(yōu)平面的距離較大。

注:紅色平面表示最優(yōu)平面，綠色立方體盒子表示數(shù)據(jù)集包圍盒，S表示點(diǎn)到最優(yōu)平面的距離跨度。

令平面P的方程如式(1)

Ax+By+Cz+D=0

(1)

(2)

式中，S=max(dmax,dmax-dmin)表示每一次迭代的距離跨度；max(*,*)表示兩者中最大者；S′表示原始數(shù)據(jù)的距離跨度。由式(2)可知：ρ與S成反比，與Nin成正比。如果減少或增加某個(gè)點(diǎn)時(shí)，ρ出現(xiàn)明顯的變化，則說明點(diǎn)是噪聲的可能性較大，反之，ρ基本無變化，則說明點(diǎn)在容忍誤差閾值范圍內(nèi)。

1.3 終止條件

對(duì)于迭代性回歸問題，終止條件對(duì)于算法精度和運(yùn)行效率都非常重要。冗余迭代將會(huì)帶來較高的時(shí)間代價(jià)，而迭代不足將可能會(huì)導(dǎo)致精度較低。因而，針對(duì)不同尺度的場景，自適應(yīng)的終止條件能夠保證在時(shí)間成本低的情況下，有效估計(jì)出模型參數(shù)。

傳統(tǒng)RANSAC算法提出一種概率估計(jì)的方法計(jì)算迭代次數(shù)，其假設(shè)在真實(shí)誤差容忍閾值αt范圍內(nèi)，迭代次數(shù)為t，則期望E(t)=γ-m，其中，m為最小采樣樣本。同時(shí)，令容忍失敗率為ε時(shí)，通過式(3)對(duì)迭代次數(shù)t進(jìn)行了保守估計(jì)

(3)

但是，γ是作為后驗(yàn)知識(shí)計(jì)算得到的，即：γ與α具有正相關(guān)性。已知α為人工給定的誤差容忍閾值，β表示最優(yōu)估計(jì)內(nèi)群點(diǎn)比率，則在足夠多的迭代次數(shù)時(shí)，滿足式(4)

(4)

由此可以看出，在傳統(tǒng)RANSAC算法中，模型參數(shù)估計(jì)精度依賴于α的人工給定精度，并且無法通過增加迭代次數(shù)t的方式減少α的誤差。

為了權(quán)衡效率與精度，本文利用空間密度估計(jì)迭代次數(shù)t，式(5)

(5)

2 試驗(yàn)分析與應(yīng)用

2.1 試驗(yàn)數(shù)據(jù)

為了驗(yàn)證方法的穩(wěn)健性，本文模擬仿真了不同內(nèi)群點(diǎn)比例的點(diǎn)云數(shù)據(jù)，用于平面擬合。試驗(yàn)數(shù)據(jù)采用正態(tài)分布與均勻分布偽隨機(jī)函數(shù)分別生成100個(gè)內(nèi)群點(diǎn)和不同比例的噪聲。其中，內(nèi)群點(diǎn)服從

(6)

式中，σ取為1.0。噪聲通過隨機(jī)數(shù)方式確定分布在平面兩側(cè)的點(diǎn)個(gè)數(shù)，其坐標(biāo)服從

(7)

式中，γ分別取為0.1～0.9。另外，本文對(duì)每個(gè)算法運(yùn)行1000次，取其平均值用于對(duì)比與分析。

在點(diǎn)云拼接或定向中，人工標(biāo)靶(如平面標(biāo)靶、球形標(biāo)靶)常被用作拼接不同掃描站或提供控制點(diǎn)，標(biāo)靶幾何中心精度將會(huì)直接影響點(diǎn)云拼接或定向的精度[25-27]。為了驗(yàn)證方法在真實(shí)場景中的可行性，本文選取球形標(biāo)靶作為研究對(duì)象，對(duì)其表面點(diǎn)云進(jìn)行球面擬合，并通過點(diǎn)云定向精度對(duì)本文算法進(jìn)行試驗(yàn)對(duì)比分析。試驗(yàn)數(shù)據(jù)使用Rigel VZ-400掃描儀進(jìn)行采集，一共有11個(gè)自制球形標(biāo)靶，其半徑有0.161 m和0.106 m兩種。球形標(biāo)靶中心到掃描儀中心最遠(yuǎn)距離為78.688 m，最近距離為4.051 m。球形標(biāo)靶中心的工程指定坐標(biāo)系采用CGCS2000坐標(biāo)系，高程采用85高程基準(zhǔn)，其坐標(biāo)見表1的第2—4列。

表1 球形標(biāo)靶中心在CGCS2000坐標(biāo)系和掃描儀坐標(biāo)系的坐標(biāo)

2.2 試驗(yàn)對(duì)比

本文分別在查全率、查準(zhǔn)率、模型參數(shù)誤差(見式(8))和迭代次數(shù)4個(gè)方面與現(xiàn)有算法(MLESAC[12]、RMLESAC[19]、AMLESAC[20]、u-MLESAC[21]和StaRSac[23])進(jìn)行對(duì)比分析

(8)

式中，nin和nout分別為召回的內(nèi)群點(diǎn)和噪聲點(diǎn)個(gè)數(shù)；Nin為真實(shí)內(nèi)群點(diǎn)個(gè)數(shù)；p為內(nèi)群點(diǎn)齊次坐標(biāo)；M為模型參數(shù)；|*,*|表示到M的歐幾里得距離。

如圖4(a)所示，本文方法的查全率隨著γ的增加呈下降趨勢(shì)，其中，在γ=0.7時(shí)，查全率最低，為0.59，并且在γ>0.5時(shí)，所有方法中表現(xiàn)最差。但是，本文方法的查準(zhǔn)率(圖4(b))對(duì)于所有γ優(yōu)于其他方法，尤其是在γ>0.5時(shí)，查準(zhǔn)率接近于1。u-MLESAC與AMLESAC算法在區(qū)分內(nèi)群點(diǎn)和噪聲方面表現(xiàn)最差。在γ<0.7時(shí)，查準(zhǔn)率較低，尤其γ<0.4時(shí)，查準(zhǔn)率與γ相同，即：算法無法剔除掉噪聲點(diǎn)。

如圖4(c)所示，對(duì)于所有γ，本文方法的模型參數(shù)精度優(yōu)于AMLESAC、MLESAC和RMLESAC算法，在γ≥0.3時(shí)，整體上要優(yōu)于u-MLESAC和StaRSac算法。需要說明的是，真實(shí)值是首先通過100個(gè)內(nèi)群點(diǎn)利用最小二乘估計(jì)方法擬合出最優(yōu)模型參數(shù)，然后計(jì)算模型參數(shù)誤差所得。

在時(shí)間成本方面，由于AMLESAC、MLESAC和RMLESAC算法的迭代次數(shù)是采用默認(rèn)固定值，其迭代次數(shù)不變(如圖4(d)所示)，且AMLESAC與MLESAC算法迭代次數(shù)相同。對(duì)于所有γ，本文方法迭代次數(shù)與StaRSac算法接近。為權(quán)衡算法精度與效率，u-MLESAC算法的誤差容忍閾值選取為5，迭代次數(shù)隨著γ的減少而逐漸增加，尤其γ<0.7時(shí)，迭代次數(shù)明顯增加。

圖4 不同內(nèi)群點(diǎn)比率(γ)的4項(xiàng)指標(biāo)對(duì)比Fig.4 The comparison between four indicators for the different γ

2.3 在點(diǎn)云定向中的應(yīng)用

球形標(biāo)靶憑借其簡單、易用、精度高等優(yōu)勢(shì)，廣泛應(yīng)用于大場景三維激光掃描點(diǎn)云定向。為了減弱球形標(biāo)靶表面噪聲對(duì)點(diǎn)云定向精度的影響，本節(jié)首先利用本文方法、u-MLESAC、AMLESAC、MLESAC和RMLESAC算法分別對(duì)不同噪聲比例的球形標(biāo)靶表面點(diǎn)云進(jìn)行定性分析，然后，根據(jù)內(nèi)群點(diǎn)對(duì)球形標(biāo)靶進(jìn)行球面擬合(本文方法擬合得到的掃描儀坐標(biāo)系中球形標(biāo)靶中心坐標(biāo)列于表1的第5—7列)，通過評(píng)定點(diǎn)云定向精度的方式進(jìn)行定量分析。如圖5所示，試驗(yàn)中球形標(biāo)靶點(diǎn)云密度分布服從文中標(biāo)靶表面點(diǎn)云密度大于其他位置。

圖5 點(diǎn)云密度分布熱圖Fig.5 The heat map of point clouds density:the darker the color,the larger the density

2.3.1 定性分析

如圖6所示，本文方法能夠剔除γ≥0.3點(diǎn)云中的噪聲(三腳架和底座)，當(dāng)γ≥0.5時(shí)，AMLESAC和MLESAC算法表現(xiàn)與本文方法基本一致，但是，γ=0.3時(shí)，算法失敗。u-MLESAC算法能夠剔除γ≥0.7時(shí)的噪聲，而RMLESAC算法無法剔除底座處的點(diǎn)云。需要說明的是，StaRSac算法穩(wěn)健性較差，并且，對(duì)于數(shù)據(jù)量較大的點(diǎn)云，效率太低，沒有用于對(duì)比試驗(yàn)。本文方法對(duì)于內(nèi)群點(diǎn)(球形標(biāo)靶)密度大于噪聲(三腳架和底座)密度時(shí)，對(duì)γ有較好的穩(wěn)健性，而其他方法對(duì)于γ更為敏感，即：γ<0.5，算法往往會(huì)失敗。需要說明的是，圖6中的γ取值是一個(gè)粗略估計(jì)值，只用于定性分析。

注：綠色表示篩選出的內(nèi)群點(diǎn)，紅色表示剔除的噪聲，黑色表示整個(gè)三角架的點(diǎn)云數(shù)據(jù)，用于可視化。理想值是球形標(biāo)靶表面點(diǎn)云。

2.3.2 定量分析

本節(jié)首先利用球形標(biāo)靶中心分別在工程指定坐標(biāo)系與掃描坐標(biāo)系的坐標(biāo)進(jìn)行坐標(biāo)轉(zhuǎn)換，統(tǒng)計(jì)空間位置殘差如式(9)所示

(9)

式中，Err3d表示空間位置誤差；Errx、Erry和Errz分別為3個(gè)坐標(biāo)軸誤差，然后與現(xiàn)有方法進(jìn)行對(duì)比分析。其中，“人工去噪”是人工交互刪除噪聲后，使用最小二乘估計(jì)法擬合球形標(biāo)靶中心坐標(biāo)，然后統(tǒng)計(jì)空間位置殘差所得。如圖7所示，本文方法的空間位置誤差整體上明顯低于現(xiàn)有算法及人工去噪后擬合計(jì)算的中心位置誤差，尤其是第5、9、10和11個(gè)標(biāo)靶都在1.0 cm以內(nèi)，u-MLESAC算法精度最低，誤差都在2.0 cm以上。

圖7 空間位置誤差分析Fig.7 The result analysis of spatial position error

綜上所述，在定量分析中，所有方法都可以對(duì)噪聲進(jìn)行不同程度的剔除，尤其是本文方法、AMLESAC和MLESAC算法能夠剔除γ≥0.5時(shí)的噪聲。但是，球形標(biāo)靶表面附近的噪聲點(diǎn)對(duì)球面擬合精度也會(huì)有較大影響。從定性分析結(jié)果可以看出：本文方法擬合得到的標(biāo)靶中心能夠更好地用于點(diǎn)云定向，其精度優(yōu)于現(xiàn)有方法。

3 總 結(jié)

本文利用內(nèi)群點(diǎn)與噪聲之間密度差異性提出了一種參數(shù)自適應(yīng)的RANSAC算法，主要解決了傳統(tǒng)RANSAC算法及其改進(jìn)版本需要先驗(yàn)知識(shí)的不足，能夠適用于不同尺度的模型回歸問題。本文首先利用加權(quán)采樣的方式優(yōu)化初始模型參數(shù)，然后借助空間密度函數(shù)來評(píng)定最優(yōu)估計(jì)模型，并更新迭代次數(shù)和誤差容忍閾值，無須任何先驗(yàn)知識(shí)。

本文針對(duì)不同內(nèi)群點(diǎn)比率的數(shù)據(jù)與已有算法分別在查全率、查準(zhǔn)率、模型參數(shù)誤差及迭代次數(shù)4個(gè)方面進(jìn)行了詳細(xì)的對(duì)比和分析。在查準(zhǔn)率和模型估計(jì)精度方面，本文方法明顯優(yōu)于其他算法，適用于圖像拼接、幾何元素參數(shù)估計(jì)以及機(jī)器人精準(zhǔn)定位等對(duì)精度要求較高的領(lǐng)域；但是，由于內(nèi)群點(diǎn)查全率較低，本文方法不適用于幾何細(xì)節(jié)修復(fù)等領(lǐng)域。