侯東興 李宗春 鄧 勇

1 信息工程大學(xué)地理空間信息學(xué)院，鄭州市科學(xué)大道62號(hào)，450001

地面三維激光掃描技術(shù)是一種高效獲取海量點(diǎn)云并運(yùn)用點(diǎn)云數(shù)據(jù)處理軟件進(jìn)行數(shù)據(jù)預(yù)處理和三維建模的技術(shù)。點(diǎn)云拼接是三維激光掃描儀數(shù)據(jù)采集和處理中的關(guān)鍵環(huán)節(jié)，影響數(shù)據(jù)預(yù)處理和建模的精度，是激光掃描技術(shù)數(shù)據(jù)處理中的熱點(diǎn)和難點(diǎn)［1－2］。點(diǎn)云拼接技術(shù)實(shí)質(zhì)是坐標(biāo)轉(zhuǎn)換，而坐標(biāo)轉(zhuǎn)換常用的方法是通過尋找3對(duì)或者3對(duì)以上的公共點(diǎn)來計(jì)算坐標(biāo)轉(zhuǎn)換參數(shù)，完成坐標(biāo)系的統(tǒng)一。在三維激光掃描模式下，公共點(diǎn)轉(zhuǎn)換方法依然可行，但是尋求公共點(diǎn)并精確計(jì)算點(diǎn)的坐標(biāo)費(fèi)時(shí)費(fèi)力，研究基于公共特征的坐標(biāo)轉(zhuǎn)換方法勢(shì)在必行［3］。

掃描場(chǎng)景中廣泛存在著平面特征。王力等［3］針對(duì)多平面特征的場(chǎng)景（平面數(shù)至少為3個(gè)），提出基于平面的坐標(biāo)轉(zhuǎn)換模型。目前，市場(chǎng)上掃描儀一般都具有整平和補(bǔ)償功能，掃描測(cè)量中不同測(cè)站坐標(biāo)系Z軸方向與鉛垂線方向相反，即不同測(cè)站坐標(biāo)系繞Z軸旋轉(zhuǎn)即能滿足三軸指向概略相同?；诖?，本文提出基于共有單平面的點(diǎn)云拼接技術(shù)：首先進(jìn)行平面的提取、擬合，精確計(jì)算平面的法向量；再運(yùn)用平面的法向量信息，利用旋轉(zhuǎn)角構(gòu)造旋轉(zhuǎn)矩陣；最后通過平面的中心點(diǎn)（視為公共點(diǎn)）進(jìn)行坐標(biāo)平移，完成點(diǎn)云拼接。

1 平面特征的提取

提取平面特征的方法一般分為人機(jī)交互法和平面特征自動(dòng)提取算法等［4－6］。人機(jī)交互法是操作員通過計(jì)算機(jī)顯示的掃描點(diǎn)云數(shù)據(jù)，將共有的平面點(diǎn)云提取出來，或者在掃描階段對(duì)共有的平面進(jìn)行單獨(dú)掃描，直接將平面點(diǎn)云提取出來。自動(dòng)提取算法包括以下兩種。

1）基于法向量聚類的平面點(diǎn)云提取。平面特征自動(dòng)提取需要提前獲取平面的大致個(gè)數(shù)。平面點(diǎn)云中的點(diǎn)存在一個(gè)共同特性，即點(diǎn)的單位法向量基本相同，它們的單位法向量與平面的單位法向量一致。因此，可以采用聚類的算法（將單位法向量是否相同作為是否歸為一類的評(píng)判標(biāo)準(zhǔn)）將平面提取出來。平面特征自動(dòng)提取的步驟為：①對(duì)原始點(diǎn)云進(jìn)行采樣；②計(jì)算點(diǎn)云的單位法向量［7］；③對(duì)點(diǎn)云的法向量進(jìn)行聚類運(yùn)算；④編號(hào)存儲(chǔ)法向量相同的點(diǎn)云。

2）基于回光強(qiáng)度的平面點(diǎn)云提取［1－2］。三維激光掃描儀不僅記錄目標(biāo)點(diǎn)的三維坐標(biāo)，還記錄照射材料的回光強(qiáng)度。目前對(duì)回光強(qiáng)度的利用相對(duì)較少，大多數(shù)情況下僅僅用于點(diǎn)云的可視化。在不同于其他材質(zhì)的平面特征中，人工標(biāo)志的反射強(qiáng)度明顯高于一般地物，可以利用點(diǎn)的回光強(qiáng)度值進(jìn)行聚類，從而將平面自動(dòng)提取出來。

2 平面擬合

2.1 求解平面參數(shù)［8－10］

平面擬合實(shí)質(zhì)上是求解平面的最優(yōu)參數(shù)。三維空間中平面的數(shù)學(xué)表達(dá)式一般為：ax＋by＋cz＋d＝0。該表達(dá)式唯一確定了一個(gè)平面，但是平面方程并不是只有一個(gè)這樣的表達(dá)式（如2ax＋2by＋2cz＋2d＝0）。為了能讓平面有唯一的表達(dá)式，定義平面方程為：

式中，a、b、c為平面擬合的待估參數(shù)（此平面不與XOY面平行）。

掃描儀獲取的點(diǎn)坐標(biāo)在x、y、z3個(gè)方向均存在誤差，誤差方程式為［11］：

式中，vxi、vyi、vzi為點(diǎn)坐標(biāo)x、y、z3個(gè)方向的誤差改正數(shù)。整理成誤差方程向量形式：

取方程的最小二乘解作為參數(shù)的估值。考慮到系數(shù)矩陣和觀測(cè)值中均含有誤差，采用矩陣的SVD（奇異值分解）［8］法求解參數(shù)最小二乘解。步驟為：①組成增廣矩陣，并進(jìn)行奇異值分解，即，其中V＝；②求解參數(shù)估值，如果V22非奇異，則；③進(jìn)行精度評(píng)定。

2.2 優(yōu)化平面參數(shù)

由于各種因素的影響，三維激光掃描儀獲取的平面點(diǎn)云必然會(huì)存在粗差點(diǎn)，特別是在平面的邊界區(qū)域。因此，為了獲取較高精度的平面參數(shù)，需要對(duì)這些粗差點(diǎn)進(jìn)行剔除，再進(jìn)行平面參數(shù)的重新估值。粗差點(diǎn)的剔除主要根據(jù)點(diǎn)到擬合平面的距離進(jìn)行判定。平面參數(shù)優(yōu)化的具體步驟為。

1）計(jì)算點(diǎn)到平面的距離。引入點(diǎn)到平面的距離公式：

式中，A＝a，B＝b，C＝－1，D＝c。

2）應(yīng)用貝塞爾公式計(jì)算距離中誤差［11］：

3）粗差點(diǎn)剔除。由于三維掃描儀獲取的點(diǎn)云數(shù)據(jù)量大，粗差點(diǎn)數(shù)量不易估計(jì)，粗差點(diǎn)對(duì)初始平面參數(shù)的影響也不易估計(jì)。設(shè)置3md為極限誤差來剔除粗差點(diǎn)，若di＞3md，此點(diǎn)作為粗差點(diǎn)剔除。

4）重新計(jì)算平面參數(shù)。利用保留下來的平面點(diǎn)云和前文方法重新計(jì)算平面參數(shù)a、b、c。

5）迭代優(yōu)化平面參數(shù)。設(shè)置3md為極限誤差來剔除粗差點(diǎn)，將2md設(shè)置為一般誤差來優(yōu)化點(diǎn)云。利用優(yōu)化的點(diǎn)云重新計(jì)算平面參數(shù)，將剔除粗差后的點(diǎn)云作為新的平面點(diǎn)云繼續(xù)迭代優(yōu)化計(jì)算。通過設(shè)置迭代次數(shù)閾值或者平面參數(shù)變化閾值，來完成平面參數(shù)的優(yōu)化。

3 解算轉(zhuǎn)換矩陣

3.1 計(jì)算旋轉(zhuǎn)矩陣

平面的法向量可由平面參數(shù)來表示。平面法向量可表示為（a，b，－1）或（－a，－b，1），兩種法向量方向相反。

不同測(cè)站下掃描獲取的同一平面點(diǎn)云的法向量是不同的。假設(shè)測(cè)站1與測(cè)站2之間的旋轉(zhuǎn)矩陣為R21（R21為3×3矩陣），那么一定滿足（a1，b1，－1）＝R21（a2，b2，－1）（設(shè)（a1，b1，－1）與（a2，b2，－1）同為背離測(cè)站方向或同為指向測(cè)站方向）。旋轉(zhuǎn)矩陣具有3個(gè)自由度，如果只將法向量視為一個(gè)公共點(diǎn)坐標(biāo)，是不能解出旋轉(zhuǎn)矩陣的。平面的法向量由平面點(diǎn)云計(jì)算得出，所以不能只將法向量視為一個(gè)點(diǎn)坐標(biāo)。

運(yùn)用兩個(gè)法向量求解一個(gè)旋轉(zhuǎn)矩陣，使法向量重合有較多方法［12］，其中較經(jīng)典的是基于羅德里格矩陣求解旋轉(zhuǎn)矩陣。其求解公式為：

式中，I為單位陣，r為向量叉乘計(jì)算的單位向量，θ為向量之間的夾角。

此旋轉(zhuǎn)矩陣R滿足（a1，b1，－1）＝R（a2，b2，－1）。但此旋轉(zhuǎn)矩陣只是將平面的法向量方向重合，而沒有將平面上的點(diǎn)云重合，不是歐氏空間轉(zhuǎn)換的旋轉(zhuǎn)矩陣。歐氏空間轉(zhuǎn)換的旋轉(zhuǎn)矩陣可由歐拉角、單位四元數(shù)、Kardan旋轉(zhuǎn)（依次繞坐標(biāo)軸Z、Y、X旋轉(zhuǎn)）等方法得到［1］。已知同一平面在不同測(cè)站的法向量，可根據(jù)Kardan 旋轉(zhuǎn)方法求取旋轉(zhuǎn)矩陣。具體步驟如下。

1）統(tǒng)一兩組法線方向（同為指向測(cè)站方向或同為背離測(cè)站方向），方法是構(gòu)造一個(gè)向量來約束平面法線方向。向量的構(gòu)造過程為：任取共有平面上一點(diǎn)P（也可取平面重心點(diǎn)或中心點(diǎn)），設(shè)測(cè)站坐標(biāo)原點(diǎn)為S，連接P、S兩點(diǎn)，方向指向測(cè)站方向（圖1）。將構(gòu)造的向量點(diǎn)乘平面法向量，如果不同測(cè)站點(diǎn)乘的結(jié)果符號(hào)相同，則將此時(shí)的法向量作為平面法向量；如果不同測(cè)站點(diǎn)乘的結(jié)果符號(hào)不同，則將其中一個(gè)測(cè)站的法向量的反向量作為平面的法向量，另一個(gè)測(cè)站的法向量不變，最終確保兩個(gè)平面法向量相對(duì)于測(cè)站具有統(tǒng)一的指向。

圖1 構(gòu)造的向量與平面法線的關(guān)系Fig.1 The relationship of constructed vector and plane normal

2）將不同測(cè)站共有平面的法向量繞Z軸旋轉(zhuǎn)到XOZ面的第一象限，分別計(jì)算旋轉(zhuǎn)后法向量與Z軸正方向的夾角α1、α2，計(jì)算待拼接測(cè)站坐標(biāo)系繞Y軸旋轉(zhuǎn)的角度Δα＝｜α1－α2｜。

3）將不同測(cè)站的共有平面的法向量繞Z軸旋轉(zhuǎn)到Y(jié)OZ面的第一象限，分別計(jì)算旋轉(zhuǎn)后法向量與Z軸正方向的夾角β1、β2，計(jì)算待拼接測(cè)站坐標(biāo)系繞X軸旋轉(zhuǎn)的角度Δβ＝｜β1－β2｜。

4）若Δα＜Δβ，分別計(jì)算不同測(cè)站的共有平面法向量在XOY面的投影與X軸正方向的夾角θ1、θ2；若Δα＞Δβ，則分別計(jì)算不同測(cè)站的共有平面法向量在XOY面的投影與Y軸正方向的夾角θ1、θ2。

5）確定解算旋轉(zhuǎn)矩陣。根據(jù)掃描儀特性，Δα、Δβ為較小的值。若Δα＜Δβ，則將兩個(gè)測(cè)站的坐標(biāo)系分別繞Z軸旋轉(zhuǎn)－θ1、－θ2，將平面法向量旋轉(zhuǎn)到XOZ面內(nèi)，再分別計(jì)算不同測(cè)站的共有平面法向量與X軸正方向的夾角，將兩個(gè)測(cè)站的坐標(biāo)系分別繞Y軸旋轉(zhuǎn)將平面法向量旋轉(zhuǎn)到與X軸平行，最后將待拼接測(cè)站坐標(biāo)系繞X軸旋轉(zhuǎn)Δβ＝β1－β2。則：

3.2 計(jì)算平移矩陣

共有平面在待拼接測(cè)站下的平面點(diǎn)云經(jīng)過旋轉(zhuǎn)后，與這個(gè)平面在固定測(cè)站下的平面點(diǎn)云平行，再經(jīng)過平移即能完成平面點(diǎn)云的重合。平移矩陣需要一對(duì)公共點(diǎn)解算：

尋找公共點(diǎn)有多種方法，如將兩個(gè)平面的中心或重心作為公共點(diǎn)，在掃描時(shí)測(cè)得一組公共點(diǎn)，或通過人機(jī)交互直接確定一組或幾組公共點(diǎn)。公共點(diǎn)的坐標(biāo)精度將直接影響點(diǎn)云拼接的精度。

4 基于共有單平面特征的點(diǎn)云拼接實(shí)驗(yàn)

4.1 仿真數(shù)據(jù)實(shí)驗(yàn)

利用VS2010編程，獲取一個(gè)平面點(diǎn)云和一個(gè)半球形點(diǎn)云，視為第一測(cè)站點(diǎn)云數(shù)據(jù)。將第一測(cè)站點(diǎn)云數(shù)據(jù)經(jīng)過坐標(biāo)轉(zhuǎn)換，作為第二測(cè)站點(diǎn)云數(shù)據(jù)。模擬兩站點(diǎn)云數(shù)據(jù)如圖2（a）所示。

第一測(cè)站點(diǎn)云經(jīng)轉(zhuǎn)換矩陣后形成第二測(cè)站點(diǎn)云，滿足：

式中，P為第一測(cè)站點(diǎn)云，Q為模擬的第二測(cè)站點(diǎn)云。坐標(biāo)轉(zhuǎn)換矩陣如表1所示。

表1 生成第二測(cè)站點(diǎn)云的轉(zhuǎn)換矩陣Tab.1 The transformation matrix of generating second station point cloud

經(jīng)過平面提取、擬合，解算兩站點(diǎn)云的旋轉(zhuǎn)矩陣。將點(diǎn)云的重心作為公共點(diǎn)求解平移矩陣。利用求解出的轉(zhuǎn)換矩陣，對(duì)第二測(cè)站點(diǎn)云進(jìn)行坐標(biāo)轉(zhuǎn)換，轉(zhuǎn)換后點(diǎn)云如圖2（b）所示。

圖2 模擬數(shù)據(jù)拼接實(shí)驗(yàn)Fig.2 Registration experiment by simulated data

平面導(dǎo)出的轉(zhuǎn)換矩陣應(yīng)滿足：

式中，Q為模擬的第二測(cè)站點(diǎn)云，P2為反求的第一測(cè)站點(diǎn)云。單平面求得的轉(zhuǎn)換矩陣如表2所示。

表2 單平面導(dǎo)出的轉(zhuǎn)換矩陣Tab.2 The transformation matrix calculate by single plane

由式（10）和式（11）可導(dǎo)出，兩組轉(zhuǎn)換矩陣應(yīng)滿足：

式中，I為單位陣。將兩組轉(zhuǎn)換矩陣代入式（12），求得的結(jié)果如表3所示。

表3 R2R1 和R2T1＋T2 的值Tab.3 The values of R2R1and R2T1＋T2

從圖2（b）可知，兩站點(diǎn)云較好地拼接在一起。從表3可以推出，反求的轉(zhuǎn)換矩陣是正確的。經(jīng)過坐標(biāo)差計(jì)算，兩站模擬數(shù)據(jù)坐標(biāo)差基本為0，證明了此方法的可行性。

4.2 實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)

點(diǎn)云采用RIEGL VZ－400掃描儀獲取，配套軟件為RiSCAN－PRO。進(jìn)行兩組實(shí)驗(yàn)，實(shí)驗(yàn)1對(duì)方凳和儀器箱分別進(jìn)行掃描，實(shí)驗(yàn)2對(duì)雕像進(jìn)行掃描。平面點(diǎn)云提取方式為對(duì)共有的平面單獨(dú)掃描，直接將平面點(diǎn)云提取出來（點(diǎn)云為方凳平面點(diǎn)云）。首先對(duì)點(diǎn)云進(jìn)行均勻采樣，保證一定數(shù)量的點(diǎn)云個(gè)數(shù)（數(shù)量過大，奇異值分解計(jì)算無法進(jìn)行），求取最優(yōu)的平面參數(shù)。然后，將求取的平面參數(shù)運(yùn)用到原始點(diǎn)云中，進(jìn)行粗差點(diǎn)剔除，將優(yōu)化后的平面點(diǎn)云重心坐標(biāo)作為公共點(diǎn)。

利用最優(yōu)的平面參數(shù)和公共點(diǎn)解算轉(zhuǎn)換矩陣，將轉(zhuǎn)換矩陣應(yīng)用到待拼接的點(diǎn)云中完成點(diǎn)云拼接。兩站平面點(diǎn)云拼接前后分別如圖4（a）、（b）所示。

實(shí)驗(yàn)1中，利用共有平面（方凳）導(dǎo)出的轉(zhuǎn)換矩陣對(duì)儀器箱進(jìn)行拼接，如圖5所示。實(shí)驗(yàn)2中，利用共有平面（雕像基座）導(dǎo)出的轉(zhuǎn)換矩陣對(duì)人物雕像進(jìn)行拼接，如圖6所示。從圖4～6看出，不同站間點(diǎn)云比較好地拼接在一起，證明此方法正確可行。

圖3 平面點(diǎn)云優(yōu)化Fig.3 The optimization of plane point cloud

圖4 平面點(diǎn)云拼接實(shí)驗(yàn)Fig.4 The registration of plane point cloud

圖5 基于平面的儀器箱拼接實(shí)驗(yàn)Fig.5 The registration of instrument container based on plane

圖6 基于平面的雕像拼接實(shí)驗(yàn)Fig.6 The registration of statue based on plane

5 結(jié) 語

基于共有單平面特征的點(diǎn)云拼接具有較強(qiáng)的適用性，模擬數(shù)據(jù)和實(shí)驗(yàn)數(shù)據(jù)均驗(yàn)證了本文算法的正確性和可靠性。由于測(cè)量誤差的存在，提取的特征往往存在誤差，相鄰測(cè)站通過特征建立的關(guān)系模型也會(huì)存在模型誤差，因此通過特征法得到的轉(zhuǎn)換參數(shù)精度相對(duì)較低。在地形測(cè)量、工程測(cè)量、雕像建模、數(shù)字城市等領(lǐng)域，容易提取和應(yīng)用的特征相對(duì)較多，特征拼接容易滿足精度要求。但在文物保護(hù)、逆向工程等精度要求較高的領(lǐng)域，提取和使用的特征精度往往難以滿足要求。將特征法和ICP 算法相結(jié)合，先用特征法進(jìn)行粗拼接，作為ICP 的初始轉(zhuǎn)換參數(shù)，然后再由ICP 迭代提高拼接精度。在下一步的工作中，將設(shè)計(jì)平面度較高的平面圓形靶標(biāo)進(jìn)行多站點(diǎn)云拼接，用以提高測(cè)量工作的效率。在多平面場(chǎng)景中，研究平面的自動(dòng)匹配以及單平面與多平面坐標(biāo)轉(zhuǎn)換模型的組合求解。

［1］李廣云，李宗春.工業(yè)測(cè)量系統(tǒng)原理與應(yīng)用［M］.北京：測(cè)繪出版社，2011（Li Guangyun，Li Zongchun.The Principles and Applications of Industrial Measuring Systems［M］.Beijing：Survey and Mapping Press，2011）

［2］王力.基于人工標(biāo)志的激光掃描數(shù)據(jù)自動(dòng)拼接技術(shù)研究［D］.鄭州：信息工程大學(xué)，2010（Wang Li.Automatic Registration of Terrestrial LiDAR Data Using Planar Targets［D］.Zhengzhou：Information Engineering University，2010）

［3］王力，李廣云，張啟福，等 激光掃描中平面擬合及坐標(biāo)轉(zhuǎn)換模型構(gòu)建［J］.測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2011，29（2）：101－104（Wang Li，Li Guangyun，Zhang Qifu，et al.Plane Fitting and Transformation in Laser Scanning［J］.Journal of Geomatics Science and Technology，2011，29（2）：101－104

［4］潘國(guó)榮，秦世偉，蔡潤(rùn)彬，等.三維激光掃描擬合平面自動(dòng)提取算法［J］.同濟(jì)大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，37（9）：1 250－1 255（Pan Guorong，Qin Shiwei，Cai Runbin，et al.Fitted Plane Automatic Extraction Algorithm of 3D Laser Scanning［J］.Journal of Tongji University：Natural Science，2009，37（9）：1 250－1 255）

［5］Biosca J，Lerma J.Unsupervised Robust Planar Segmentation of Terrestrial Laser Scanner Point Clouds Based on Fuzzy Clustering Methods［J］.ISPRS Journal of Photogrammetry and Remote Sensing，2008，63（1）：84－98

［6］戴楠，李傳榮，蘇國(guó)中，等.激光點(diǎn)云提取建筑物平面目標(biāo)算法研究［J］.微計(jì)算機(jī)信息，2010，26（7）：205－207（Dai Nan，Li Chuanrong，Su Guozhong，et al.Algorithms Study of Building Planar Objects Detection on LiDAR Data［J］.Microcomputer Information，2010，26（7）：205－207）

［7］李澤宇，李德華，胡漢平，等.基于八叉樹的三維散亂數(shù)據(jù)點(diǎn)法矢的估計(jì)［J］.計(jì)算機(jī)與數(shù)字工程，2000，28（4）：62－65（Li Zeyu，Li Dehua，Hu Hanping，et al.Estimation of Normal Vector of 3DScattered Point Based on Octree［J］.Computer ＆Digital Engineering，2000，28（4）：62－65）

［8］魯鐵定，周世健，張立亭，等.基于整體最小二乘的地面激光掃描標(biāo)靶球定位方法［J］.大地測(cè)量與地球動(dòng)力學(xué)，2009，29（4）：102－105（Lu Tieding，Zhou Shijian，Zhang Liting，et al.Sphere Target Fixing of Point Cloud Data Based on TLS［J］.Journal of Geodesy and Geodynamics，2009，29（4）：102－105）

［9］官云蘭，程效軍，施貴剛.一種穩(wěn)健的點(diǎn)云數(shù)據(jù)平面擬合方法［J］.同 濟(jì) 大 學(xué) 學(xué) 報(bào)：自 然 科 學(xué) 版，2008，36（7）：981－984（Guan Yunlan，Cheng Xiaojun，Shi Guigang.A Robust Method for Fitting a Plane to Point Clouds［J］.Journal of Tongji University：Natural Science，2008，36（7）：981－984）

［10］官云蘭，劉紹堂，周世健，等.基于整體最小二乘的穩(wěn)健點(diǎn)云數(shù)據(jù)平面擬合［J］.大地測(cè)量與地球動(dòng)力學(xué)，2011，31（5）：80－83（Guan Yunlan，Liu Shaotang，Zhou Shijian，et al.Robust Plane Fitting of Point Clouds Based on TLS［J］.Journal of Geodesy and Geodynamics，2011，31（5）：80－83）

［11］隋立芬，宋力杰.誤差理論與測(cè)量平差基礎(chǔ)［M］.北京：解放軍出版社，2004（Sui Lifen，Song Lijie.Error Theory and Base of Surveying Adjustment［M］.Beijing：PLA Press，2004）

［12］楊凡，李廣云，王力.三維坐標(biāo)轉(zhuǎn)換方法研究［J］.測(cè)繪通報(bào)，2010（6）：5－7（Yang Fan，Li Guangyun，Wang Li.Research on the Methods of Calculating 3DCoordinate Transformation Parameters［J］.Bulletin of Survey and Mapping，2010（6）：5－7）