林 森，張 強

（1.沈陽理工大學 自動化與電氣工程學院，遼寧 沈陽 110159；2.遼寧工程技術大學 電子與信息工程學院，遼寧 葫蘆島 125105）

1 引 言

隨著三維（3D）激光掃描設備的飛速發(fā)展，點云已經成為表征三維世界的主要數據格式。在實際測量過程中，由于被測物體的幾何形狀和測量方法的限制，測量設備需要從不同角度對被測物體進行多次定位測量，每次測量只能捕獲3D對象的部分點云，為了獲得完整的3D對象，必須使用點云配準技術。點云配準就是把在不同視角下獲取的多個3D點云數據，旋轉平移到同一坐標系下。目前，3D點云配準技術是逆向工程、數字醫(yī)學圖像、表面質量檢測等領域的重要基礎［1］。建立產品三維數據模型屬于逆向工程的一個應用領域，點云配準結果的好壞直接影響建模的正確性。

由于激光設備獲取的3D物體點云部分區(qū)域重疊，且?guī)в性肼暫彤惓｜c，針對該類點云配準問題，一般分為粗配準和精確配準兩個過程。粗配準就是估計兩個點云之間的初始轉換，分為深度學習配準方法和優(yōu)化的配準方法，對于未知產品，若產品與訓練數據存在較大的分布差異，深度學習配準方法性能會急劇下降，而優(yōu)化的配準方法不需要訓練數據，能更好的應用于產品建模［1］。基于優(yōu)化的粗配準又分為三類：手動配準、用投票法則的配準和用局部特征的配準，其中局部特征的配準方法因效率高且能自動配準，所以應用更廣泛。楊穩(wěn)等人［2］提出層次優(yōu)化的顱骨點云配準，該方法根據特征點的特征序列確定初始對應點對，能快速配準完整的顱骨點云，但針對更復雜的部分重疊物體點云，計算特征序列難以精準描述點的特征；Feng［3］等人提出一種基于局部特征描述符（Local Point Feature Histogram，LPFH）的粗配準，相較于常用描述符，其配準精度更高，但是粗配準框架并不剔除異常點，對于數據量大的復雜點云，整體效率低，并且難以完成配準；李新春等人［4］提出一種基于鄰域特征點提取和匹配（Neighborhood Characteristic Point Extraction and Matching，NCCP）的點云配準，提高了精度，并解決了噪聲干擾問題，但該方法在提取和匹配特征點時，兩次使用最小二乘法計算對應點的曲率，耗費了大量時間。

粗配準結果為精確配準提供初始數據。對于精確配準階段，最經典、應用最廣泛的方法是Besl等人［5］在1992年提出的ICP算法，基本原理是，最小化一個點云與另一個變換后點云之間的誤差函數，迭代更新剛體變換，直到獲得收斂為止，最后得到旋轉平移矩陣。但是，ICP算法很難找到需要大變換點云的旋轉平移矩陣，從而陷入局部最小值，計算效率低，并且很難克服噪聲。通過粗配準，可以解決ICP算法陷入局部最小值和魯棒性差的問題，并且對于初始位置好的點云，還可以減少ICP的迭代次數來快速地配準兩片點云，從而提高配準效率。針對ICP算法本身，也有各種改進。Choi等人［6］提出用k維樹改進 ICP（K-dimensional tree Iterative Closest Point，KICP），通過加快搜索最近點的速度，提高了效率；Li等人［7］給ICP迭代過程加入動態(tài)因子，在保證配準精度情況下，提高了ICP的收斂速度；王賓等人［8］在精確配準階段提出基于雙向距離比例的ICP算法，提高了配準精度；Zhang等人［9］提出快速和魯棒的ICP算法，基于Welsch函數改進誤差度量，并使用具有Anderson加速功能的Majorization-Minimization算法將其最小化，提高了配準精度和效率，并且使配準框架具有更強的魯棒性。

由于3D物體點云需將粗配準的匹配點對作為精確配準的輸入數據，才能完成有效的點云配準［10］，因此粗配準結果的好壞，會直接影響精確配準的效率和精度，上述研究直接應用于部分重疊點云，配準效果不好，因此，針對部分重疊點云配準問題，提出應用鄰域點信息描述與匹配的點云配準方法，即基于NDM的KICP算法（NDMKICP）。首先，將曲率變化、測量角度和特征值性質指標組合，提取特征點；其次，改進法向量夾角，點密度和曲率，提出多尺度矩陣描述符；再次，提出用k維樹初步建立匹配關系，并針對點的幾何信息和歐氏距離，兩次剔除錯誤匹配點對，完成粗配準；最后，用k維樹改進ICP算法完成精確配準。該配準方法能完成3D對象表面點云配準，顯著提高了配準的效率和精度，并在斯坦福點云數據庫下，驗證了算法的優(yōu)勢和魯棒性。

2 點云的粗配準

以往的粗配準研究都從單個鄰域點信息出發(fā)，進行高維度的描述，但是配準效率低。本文提出從點云數據的協(xié)方差矩陣出發(fā)，結合點的多個鄰域點幾何信息，提取特征點后，進行低維度的描述。

2.1 特征點提取

提取特征點的目的是在點云數據中獲取具有某些約束的點，作為后續(xù)粗配準工作的輸入點集。提取特征點的關鍵是具有旋轉平移不變性，使用單一的特征往往會導致配準計算量冗余、特征信息不完整［11］。因為曲率變化、測量角度和特征值性質的計算效率高、含有特征信息更多［12-13］，所以將這三個特征組合并改進，提出三尺度特征點（Curvature change，Measuring angle and Eigenvalue property，CME）。

點的曲率反映了曲面偏離平面的程度。對于給定的點pi，鄰域的深度變化越劇烈，以點pi為中心的區(qū)域特征就越明顯。如圖1（a）所示，點A、B和C分別位于平滑、略有變化和急劇變化的區(qū)域，（b）～（d）表示不同半徑曲面的切平面。對點A，在半徑較小的情況下，局部相鄰表面近似為平面，曲率近似為零。而點B，C的同心圓半徑增加，鄰域表面會迅速變化，且曲率增加。因此，可以根據不同半徑的曲率變化選擇表示某些區(qū)域特征的點。

圖1 表面不同點的示例Fig.1 Example of different points in surface

設置查詢點pi，在其半徑r鄰域內構造加權協(xié)方差矩陣：

式中，pj為鄰域點，為權重。

求解公式（1）的特征值λj1≥λj2≥λj3，曲率ci計算如下：

基于三個不同半徑rj(j=1，2，3)計算點pi的三個不同曲率。對于彎曲表面，c1，c2，c3大于零，而對于平坦表面等于零。根據公式（3），圖1（a）中的A點被刪除，而B和C點則被保留。

式中εc=0.1為設定閾值［12］。

當物體的表面急劇變化時，由于測量角度的問題，相鄰點的密度可能會減小，并且點坐標誤差可能會增加［14］，這些點不能視為關鍵點，例如C點，用公式（4）剔除曲率和大于φ=1.6［14］的不良點：

通過求解協(xié)方差矩陣獲得的特征值具有一定的幾何意義［13］，特征值的大小可以看作是橢球軸的長度，特征矢量e1，e2，e3可以看作是局部坐標系的軸，如圖2所示。橢球的形狀是對鄰域內點分布的抽象描述。如果鄰點沿某個方向密集分布，則該方向為橢球第一主方向。點分布稀疏的方向是第二主方向，點分布極為稀疏是第三主方向。為了避免檢測到沿主方向有相似伸展的特征點，選擇特征值變化量大于ρ=1.25［13］的點：

圖2 特征值和特征向量的幾何意義Fig.2 Geometric significance of eigenvalue and eigenvectors

式中，λ11，λ12是第一個半徑r1對應協(xié)方差矩陣的特征值。

對于待配準的輸入點云，給定三個半徑r1、r2、r3，取r2=1.1r1，r3=1.2r1，滿足公式（3）～式（5）的點為CME特征點。

2.2 多尺度矩陣描述符

由于高維度的描述計算復雜度高，而單一的低維度描述特征識別度低，特征信息量少。因此，從點云表面的基本幾何特征出發(fā)［3-4，15］，選取計算復雜度最低的法向量、點密度和表面曲率進行改進，提出多尺度低維度的描述符。

設P是待配準的點云，它包含N個點{p1，p2，...，pN}。給定一個查詢點pi（i∈[1，N]），其鄰域是一個以pi為中心，半徑為r的球體，用pij={pik|k=1，2，3，...，j}表示球體內pi的相鄰點，其中k是pi相鄰點個數。設是pik質心，在球體內建立協(xié)方差矩陣Cov2(pi)：

pi的法向量ni可以用以下等式表示：

式中，λi0，λi1，λi2是特征值，V0是λi0對應的特征向量?；诖?，提出三個尺度的幾何描述：

（1）對于pi和pik，用兩者法向量的cos值作為第一個尺度，用來描述點云表面區(qū)域的起伏變化，如圖3所示，其中θ是向量夾角。計算方法見公式（8）。

圖3 法向量夾角Fig.3 Normal vector included angle

式中，|·|表示向量的模，ni是查詢點pi的法向量，nik是相鄰點pik的法向量，F1(pik)是ni和nik的cos值。F1(pik)的范圍［-1，1］，將此范圍劃分為25個直方圖，并通過計算落入對應直方圖中的數目來形成第一個尺度矩陣。

（2）文獻［16］證明3D投影對噪聲具有魯棒性，因此，對于點pi，其法向量為ni，將半徑為r的球體內點集pik投影到垂直于ni的切面。投影點與pi的歐氏距離作為第二個尺度，用來描述點密度，如圖4所示，圖中L為垂直于ni的切面，計算方法見公式（9）：

圖4 點密度特征圖Fig.4 Illustration of point density feature

式中，||·||表示歐氏距離，pi為查詢點，pik為鄰域點，ni是pi的法向量。最后，將F2(pik)歸一化并劃分為15個直方圖，作為第二個尺度矩陣。

（3）對于點集P中的每個點pi，用表面曲率作為第三個尺度，描述曲面在查詢點處的彎曲程度，計算方法見公式（10）：

式中，λi0≤λi1≤λi2是特征值。將F3(pik)歸一化并劃分為20個直方圖，作為第三個尺度矩陣。

法向量、點密度和表面曲率都具有旋轉平移不變性，描述的局部特征是唯一的［17］。因此，可以將上文三個尺度矩陣組合，形成特征點的多尺度矩陣描述符（Multi-scale Matrix Descriptor，MSM），如公式（11）所示。

公式（11）中，N是特征點數量，VMSM是N×60維的矩陣，j為鄰域點數量，·為四舍五入取整，|·|為絕對值，r為鄰域半徑，為曲率的平均值。vote(x)為定義的函數，表示將第x個直方圖記為1。

2.3 獲取匹配點對

為了提高特征匹配的識別度，提出以下方法：

（1）設提取特征點的點集為Pa和Qa，兩片點云的多尺度矩陣描述符為F(Pa)，F(Qa)。建立k維樹，搜索F(Pa)和F(Qa)歐氏距離最小的子向量，把向量所對應的點，初步建立起匹配關系，記為：

式中，Ni為初步匹配點對數量。通過k維樹初步建立匹配關系，其替代了常用方法中需要遍歷全部點過程，提高了效率，同時也替代了設定歐氏距離閾值獲取匹配點對的方式，解決了因歐氏距離閾值設置不合理而導致刪除正確匹配點對的問題［9］。

（2）對于點對集K1，內任意點pi對應的法向量夾角、點密度、曲率，分別為F1(pj)、F2(pj)、F3(pj)。內任意點qi對應的法向量夾角、點密度、曲率，分別為F1(qj)、F2(qj)、F3(qj)。滿足公式（13）則認為對應關系正確，建立如公式（14）所示的匹配點對集K2。

式 中，|·|為 絕 對 值，σ1=σ2=σ3=0.1為 設 定閾值［18］。

式中，Nj為匹配點對數量。

（3）使用剛性距離約束檢驗點對集K2，在K2任選兩個匹配點對(pj，qj)，(pk，qk)，根據文獻［4］設定閾值τ=0.005，若匹配點對不滿足公式（15），則從集合K3中剔除，直到遍歷K2，得到精煉的匹配點對集K3。

式中，||·||表示歐氏距離。

式中，Nk為經過兩種配對參數精煉后獲得的正確匹配點對的數量。

最后對匹配點對集K3用單位四元數算法［18］，獲取旋轉矩陣Ra和平移向量Ta。

3 精確配準

對于源點云Q，根據初始配準獲取的旋轉矩陣Ra和平移向量Ta，用公式（17）變換得到Q'，將P和Q'作為輸入點云。

為了得到精度更高的配準結果，本文使用KICP算法進行精確配準，具體步驟如下：

（1）輸入待配準的點云P和Q'，設定終止閾值η［5］和最大迭代次數；

（2）對于Q'中一點，用k維樹在P中尋找歐氏距離最近的對應點，得到對應點對，直到遍歷Q'；

（3）用SVD法計算對應點對的旋轉矩陣Rk和平移向量Tk，然后計算公式（18）的目標函數；

（4）將步驟（3）的旋轉矩陣Rk和平移向量Tk應用于Q'，計算Q'k+1=RK·Q'k+TK，用Q'k+1替代Q'；

（5）計算步驟（3）中相鄰兩次目標函數的變化，若Ek+1-Ek＜η或達到設定的最大迭代次數時，則終止迭代，否則轉向步驟（2）。

本文的粗配準和精確配準總流程可歸結為圖5所示。

圖5 算法總流程圖Fig.5 Algorithm flow chart

4 實驗結果及分析

本部分實驗均在配置為Intel core i5 2.3GHz CPU、8GB內存的計算機上完成。

4.1 特征點提取結果的比較

為了驗證提取CME特征點的必要性，采用斯坦福大學的Bunny（35947）和Dragon（56053）點云模型，將CME約束提取特征點與曲率變化（Curvature Change，CC）、測量角度（Measuring Angle，MA）、特征值性質（Eigenvalue Property，EP）以及兩兩組合的方法比較，分別討論7類提取方法獲取的特征點數量、對粗配準效率和精度的影響，如圖6所示。人為選定點云第一個點為初始查詢點。

由圖6（a）可以看出，使用單一特征提取的特征點數量多于兩個特征組合，而CME約束提取特征點數量最少，證明特征約束越多，提取特征點數量越少。由圖6（b）可知，特征點數量越多，在計算描述符和獲取匹配點對時越耗時，CME特征點數量最少，所以粗配準效率最高。由圖6（c）可知，由于使用單一特征提取的點集內存在大量異常點，故配準誤差大，精度低；而使用兩個特征組合配準精度明顯提升；CME特征點提取效果最好，因此粗配準誤差最小。綜上證明，使用CME特征點會讓兩組待配準的點云有更好的初始位置，且配準效率最高。

圖6 特征點粗配準結果對比Fig.6 Influence of feature points on initial registration result

4.2 描述符描述結果的比較

為了證明構造多尺度矩陣描述符（MSM）的必要性，將MSM與法向量（Normal Vector，NV）矩陣、點密度（Point Density，PD）矩陣、曲率（Curvature，CU）矩陣以及兩兩組合構成的描述符進行對比，如圖7所示，并討論7種描述符計算效率以及對粗配準結果的影響。

圖7 描述符對比實驗Fig.7 Descriptor comparison experiment

由圖7（a）可以看出，相同模型下7種描述粗配準時間最大與最小相差都約1s，證明描述符的三個子尺度矩陣計算效率都很高，提出的MSM描述符時間計算復雜度低。由圖7（b）可知，對于簡單的Bunny模型，使用MSM描述符精度雖然最小，但相比于其他描述符，配準精度提升不大，所以提出的三個子尺度矩陣及其組合的描述符，都能完成簡單點云模型的配準，證明了三個子特征的有效性。而較復雜的Dragon模型配準，使用一個子尺度矩陣描述特征點描述性有限，所以粗配準誤差最高；使用兩個子尺度矩陣粗配準精度明顯提升；而MSM粗配準精度最高，證明MSM的描述性最強，獲取的匹配點對幾何特征最接近。

4.3 物體表面點云配準比較

本節(jié)使用的物體點云是Leica P50掃描儀在實際環(huán)境不同視角下對實驗室某機器人（Robot）進行掃描獲取的點云數據［19］，點數為296 356、269 132。評價指標為配準總耗時和配準誤差（ξRMSE），ξRMSE定義［9］見公式（19）。

為驗證提出的配準方法有更高的配準效率和精度，將本文的NDM-KICP算法與經典ICP算法、粗配準不同但精確配準相同的基于LPFH的KICP算法（LPFH-KICP）［3］、粗配準和精確配準都不同的基于NCCP的雙向k維樹ICP（Bidirectional k-d tree ICP，BKICP）算 法（NCCPBKICP）［4］、基于半正定的隨機算法（Semidefinite-Based Randomized Approach，SDRSAC）［20］、采樣一致性無損檢測（Sampling Consensus and Nondestructive Testing，SAC-IA NDT）算 法［21］做對比，配準結果如圖8，配準數據見表1。

圖8 Robot配準效果Fig.8 Registration result of Robot

圖8整體來看，經典ICP無法完成有效的點云配準，SAC-IA NDT出現明顯的錯配，其他4種配準算法都能完成較好的配準；觀察左臂、胸口屏幕和底座三處細節(jié)進行對比，LPFH-KICP和NCCP-BKICP在左臂和胸口屏幕按鈕處配準有微小的偏差，在底座拼接處有明顯的錯位；SDRSAC在左臂處有明顯的錯位，胸口屏幕按鈕處有微小偏差；而本文算法在這三個部分的配準線條更加流暢，配準效果最好。由于Robot點云數量大且重疊率低，即使配準框架的計算復雜度低，配準也會非常耗時。通過表1定量分析，本文算法比經典ICP的配準精度高2個量級。與LPFH-KICP算法相比，本文算法配準效率提高了40%，配準精度提高了35%。與NCCP-BKICP算法相比，本文算法配準效率提高了51%，配準精度提高了40%。與SDRSAC算法相比，本文算法配準效率提高了58%，配準精度提高了29%。與SAC-IA NDT算法相比，本文算法配準效率提高了54%，配準精度提高2個量級。所以，本文算法在配準效率和精度上有較大優(yōu)勢，證明本文三重約束提取特征點和多特征低維度描述符組合，又兩次剔除錯誤匹配點對的方法，使配準效果更好，效率更高。

表1 Robot模型配準數據Tab.1 Registration data of Robot model

綜合上述實驗，本文算法可以完成具有較大初始位置差異和部分區(qū)域重疊的實際物體表面點云的配準，并且有較高的配準效率和精度。

4.4 數據庫模型點云配準實驗

為了驗證本文NDM-KICP算法在其他點云模型上也具有較高的配準效率和精度，采用斯坦福大學模型設置兩組實驗。第一組實驗考慮實際采集的點云有殘缺，源點云為Bunny（35947）模型，目標點云為源點云在三維空間做（Re，te）的大幅度旋轉平移變換，對目標點云隨機剔除5%，10%，15%，20%的點數。第二組實驗考慮實際采集點云數據帶有噪聲，源點云為Dragon（52368）模型，目標點云為源點云做（Re，te）變換，對目標點云分別加入0.1ˉr、0.2ˉr、0.3ˉr、0.5ˉr的高斯白噪聲，其中ˉr為平均密度。配準耗時是粗、精配準總時間，配準精度使用變換矩陣相對誤差［8］，如式（20）所示。Bunny模型配準結果圖和配準數據見圖9和表2，Dragon模型配準結果圖和配準數據見圖10和表3。

表2 Bunny模型配準數據Tab.2 Registration data of Bunny model

圖9 Bunny配準效果Fig.9 Registration result of Bunny

式中，ξR是相對旋轉誤差，ξt是相對平移誤差，（Rk，tk）是精確配準求解的旋轉平移矩陣，（Re，te）是設定的旋轉平移矩陣，‖·‖F是F范數，‖·‖2是2范數。

由于斯坦福模型的結構相對簡單，點云表面多處特征相似，若提取的特征點代表性不強，描述符描述性差，會導致錯誤匹配點對急劇增加，從而影響精確配準的結果和效率，并且對于數據量小的點云，粗配準提取特征點的耗時較少，此時描述符的耗時會直接影響配準效率。

由圖9可以看出，對于不同數據殘缺（剔除點云比例）的Bunny模型點云，經典ICP無法完成有效的配準，其他算法配準效果良好。根據表2定量分析，6種算法的相對平移誤差相差不大。相比于經典ICP，本文算法的相對旋轉誤差降低3～4個量級。與LPFH-KICP算法相比，本文算法的配準耗時降低了26%～36%，相對旋轉誤差減少了8%～66%。與NCCP-BKICP算法相比，本文算法的配準耗時降低55%～64%；在相對旋轉誤差上，剔除15%～20%點云時，本文算法減少43%～66%。與SDRSAC算法相比，本文算法的配準耗時降低60%～64%，相對旋轉誤差降低82%～99%。與SAC-IA NDT算法相比，本文算法的配準耗時降低92%～94%，相對旋轉誤差降低89%～99%?？梢钥闯?，本文算法在配準效率和精度上均有優(yōu)勢，并且，隨著剔除點云比例的增加，本文算法能保持較穩(wěn)定的配準精度，而LPFH-KICP和NCCP-BKICP的相對旋轉誤差不斷變大，SDRSAC和SAC-IA NDT的相對旋轉誤差起伏較大，由此可以證明，本文算法有更強的魯棒性。

根據圖10和表3，對帶有不同噪聲的Dragon點云，經典ICP仍無法完成有效的配準，配準精度較本文差4～5個量級。且6種算法相對平移誤差相同。相比于LPFH-KICP算法，本文算法的配準耗時減少了3%～18%，相對旋轉誤差減少了1%～33%。相比于NCCP-BKICP算法，本文算法的配準耗時減少17%～30%，相對旋轉誤差減少了1%～19%。相比于SDRSAC算法，本文算法的配準耗時降低16%～53%，相對旋轉誤差降低14%～77%。相比于SAC-IA NDT算法相比，本文算法的配準耗時降低83%～88%，相對旋轉誤差降低2～3個量級。由此可以證明，對于帶噪聲的點云，本文算法的配準效率和精度仍有優(yōu)勢，這是由于本文多重幾何特征約束的粗配準方法，在噪聲環(huán)境下有更好的魯棒性，能給精確配準提供更好的初始點云。

表3 Dragon模型配準數據Tab.3 Registration data of Dragon model

圖10 Dragon配準效果Fig.10 Registration result of Dragon

通過對斯坦福點云模型的配準結果分析，可以確定在帶有數據缺失和噪聲環(huán)境下，對于不同的點云模型，本文算法仍有較高的配準效率和精度，并且有較好的魯棒性。

5 結 論

為了提高點云配準技術的配準精度和效率，對3D點云數據局部坐標的幾何信息深入分析，本文提出了應用鄰域點信息描述與匹配的點云配準方法。首先，針對點云表面深度變化，提取特征點；其次，改進法向量夾角，點密度和曲率三個幾何特征，提出了一個多尺度矩陣描述符；再次，提出用k維樹初步建立匹配關系，并結合幾個特征和剛性距離約束，剔除錯誤匹配點對，實現粗配準，最后，在粗配準基礎上，用k維樹ICP算法完成精確配準。對于實物點云配準，本文算法比經典ICP配準精度提高2個量級，相比于其他算法，配準精度、效率至少提升29%、54%，證明本文算法有較高的配準精度和效率。算法的優(yōu)勢還體現在：一是對于不同數據缺失和帶有不同噪聲的點云，精度和效率優(yōu)勢不受影響；二是在模型數據量較大時，配準效率優(yōu)勢更加明顯。因此本文算法可以應用于3D物體表面點云的配準中。進一步優(yōu)化配準方法，并將本文算法應用于其他環(huán)境中，如生物醫(yī)學中顱骨、膝蓋等模型的配準問題，是下一步研究的方向。