劉支亮李明峰陸建華

（1.蘇州市測(cè)繪院有限責(zé)任公司，江蘇 蘇州215006;2.南京工業(yè)大學(xué)測(cè)繪科學(xué)與技術(shù)學(xué)院，江蘇 南京211816）

0 引 言

在逆向工程、工業(yè)檢測(cè)和工程測(cè)量中，圓柱面的擬合問(wèn)題十分常見(jiàn)[1]。三維激光掃描技術(shù)突破了傳統(tǒng)單點(diǎn)測(cè)量模式[2]，與傳統(tǒng)圓柱面擬合數(shù)據(jù)獲取手段相比，三維激光掃描能夠快速獲取大量點(diǎn)云數(shù)據(jù)。通過(guò)圓柱面點(diǎn)云坐標(biāo)數(shù)據(jù)的擬合處理，可得圓柱面在空間直角坐標(biāo)系中的幾何方程。根據(jù)圓柱面的幾何關(guān)系，僅考慮觀測(cè)向量的隨機(jī)誤差，采用最小二乘（Least Squares，LS）原理，構(gòu)建經(jīng)典的Gauss-Marko模型計(jì)算幾何參數(shù)是最常用的擬合算法[3-5]。實(shí)際上，受儀器、掃描環(huán)境和人員等因素影響，圓柱面擬合平差模型中所有觀測(cè)值都含有隨機(jī)誤差，即觀測(cè)向量與系數(shù)矩陣均存在誤差，這類模型稱為EIV（Errors-In-Variables）模型[6]?？傮w最小二乘（Total Least Squares，TLS）法除了考慮觀測(cè)值誤差，還顧及系數(shù)矩陣的誤差，能夠得到EIV模型的最優(yōu)估計(jì)[7-9]。加權(quán)總體最小二乘（Weighted TLS，WTLS）法[10]進(jìn)一步考慮了觀測(cè)值點(diǎn)位精度的差異，并根據(jù)一定準(zhǔn)則賦予權(quán)值參與求解，從而獲得更為可信的參數(shù)解。袁慶等[11]基于點(diǎn)位精度賦予權(quán)重，通過(guò)WTLS方法解算點(diǎn)云數(shù)據(jù)平面擬合參數(shù)。陳瑋嫻等[12]將WTLS方法應(yīng)用于激光掃描球型標(biāo)靶的擬合，根據(jù)點(diǎn)云激光回波強(qiáng)度特性確定先驗(yàn)權(quán)，提高了參數(shù)解的精度。本文根據(jù)點(diǎn)到圓柱面的距離與點(diǎn)位精度的相關(guān)性，提出顧及距離的WTLS算法，即在迭代時(shí)自適應(yīng)修正觀測(cè)值權(quán)，使得權(quán)值的確定更加可靠。同時(shí)，由于獲取的點(diǎn)云數(shù)據(jù)中不可避免存在粗差點(diǎn)，故以3倍距離標(biāo)準(zhǔn)差為閾值，剔除粗差點(diǎn)，構(gòu)造穩(wěn)健的加權(quán)總體最小二乘（Robust Weighted TLS，RWTLS）圓柱面擬合算法，并通過(guò)實(shí)例數(shù)據(jù)驗(yàn)證算法的有效性。

1 圓柱面擬合數(shù)學(xué)模型

圓柱擬合參數(shù)包括大小參數(shù)和定位參數(shù)。大小參數(shù)由圓柱的半徑?jīng)Q定，所以圓柱面的標(biāo)準(zhǔn)方程x2+y2=r2只有一個(gè)半徑參數(shù);定位參數(shù)由圓柱中軸線的位置和姿態(tài)參數(shù)兩部分構(gòu)成。結(jié)合圓柱的標(biāo)準(zhǔn)坐標(biāo)系與測(cè)量坐標(biāo)系之間的關(guān)系，圓柱面的定位參數(shù)實(shí)質(zhì)上就是兩個(gè)坐標(biāo)系的坐標(biāo)轉(zhuǎn)換參數(shù)，中軸線位置參數(shù)即坐標(biāo)系的平移參數(shù)，姿態(tài)參數(shù)即坐標(biāo)系的旋轉(zhuǎn)參數(shù)。

圓柱標(biāo)準(zhǔn)坐標(biāo)系與測(cè)量坐標(biāo)系之間的轉(zhuǎn)換關(guān)系為:

由于標(biāo)準(zhǔn)坐標(biāo)系在xoy平面內(nèi)投影是一個(gè)圓，x軸與y軸方向不需要嚴(yán)格指定，只需兩個(gè)旋轉(zhuǎn)參數(shù)εX、εY就能實(shí)現(xiàn)姿態(tài)參數(shù)的確定，此時(shí)εZ=0，z軸指向天頂方向。另外，標(biāo)準(zhǔn)坐標(biāo)系的原點(diǎn)可為中軸線上任意一點(diǎn)，即只需兩個(gè)平移參數(shù)就可以確定位置，文中令。那么，原來(lái)所需求取的7個(gè)參數(shù)，實(shí)際只需要求取ΔX、ΔY、εX、εY、r共5個(gè)參數(shù)。

結(jié)合式（1），可構(gòu)造目標(biāo)函數(shù)為:

式（2）中，r為 圓 柱半 徑 （r＞ 0），R=，且

將式（2）在近似值（ΔX0、ΔY0、（εX）0、（εY）0、r0）處以泰勒級(jí)數(shù)展開(kāi)，舍去二階及以上項(xiàng)，有:

則根據(jù)最小二乘原理，可列誤差方程為:

式（4）中，ξ表示參數(shù)向量的改正值，ξ=，V L表示觀測(cè)值誤差向量，B為系數(shù)矩陣，B=，L為觀測(cè)向量，

選取合適的參數(shù)初始值，求解參數(shù)向量改正值ξ，利用上次得到的參數(shù)近似值，加上本次迭代的ξ，作為新的近似值，重新計(jì)算ξ。經(jīng)迭代收斂，最后得到5個(gè)參數(shù)的最佳估值，進(jìn)而可得圓柱方程。

2 RWTLS算法擬合圓柱面

2.1 WTLS算法基本原理

顧及系數(shù)矩陣的隨機(jī)誤差，建立系數(shù)矩陣與觀測(cè)向量均含有誤差的EIV模型，誤差方程為:

設(shè)觀測(cè)值誤差向量和系數(shù)矩陣誤差向量服從的如下分布:

式（6）中，vec（·）表示矩陣的拉直變換，V L和vec（V B）分別為觀測(cè)值和系數(shù)矩陣的誤差向量，Q L、Q B為對(duì)稱陣，Q L為觀測(cè)值的協(xié)方差陣，且Q L，系數(shù)矩陣的協(xié)因數(shù)陣Q B可以用行矩陣和列矩陣表示:

加權(quán)總體最小二乘估計(jì)準(zhǔn)則為:

其中，v B=vec（V B）。

2.2 權(quán)陣的定義

激光點(diǎn)與對(duì)象越貼近，表明測(cè)量精度越高，兩者的相關(guān)性越強(qiáng);反之，則說(shuō)明受其他因素干擾較多，兩者相關(guān)性較弱[13]。因此，可通過(guò)激光點(diǎn)到圓柱面的距離確定該點(diǎn)的權(quán)重。根據(jù)標(biāo)準(zhǔn)圓柱方程，由式（9）計(jì)算各點(diǎn)到圓柱面的距離。再按式（10）將距離轉(zhuǎn)換成[0，1]的數(shù)值，并將其作為各點(diǎn)對(duì)應(yīng)的權(quán)值P i。

式（9）—（10）中，D i為點(diǎn)到圓柱面的距離，max（D i）為點(diǎn)到圓柱面距離的最大值，δ為預(yù)設(shè)的一個(gè)小數(shù)值，通常取10-6，使P i不為0。那么，考慮每個(gè)點(diǎn)的點(diǎn)位精度差異，有系數(shù)矩陣B的行向量和觀測(cè)向量L權(quán)陣為。假設(shè)待求圓柱面擬合參數(shù)相互獨(dú)立，則系數(shù)矩陣的列向量權(quán)陣為單位陣，即P0=I5。

2.3 RWTLS求解及精度評(píng)定

利用Lagrange乘數(shù)法，可構(gòu)造Lagrange目標(biāo)函數(shù)為:

式（11）中，λ為L(zhǎng)agrange乘數(shù)子列向量?；谶x權(quán)迭代[14]的思想，利用Schaffrin B等[10]，[15]的求解算法，以3倍標(biāo)準(zhǔn)差作為閾值剔除粗差點(diǎn)，構(gòu)造一種RWTLS算法。算法迭代計(jì)算步驟如下:

（1）由TLS解作為圓柱參數(shù)估值，計(jì)算點(diǎn)到圓柱面的距離，得到系數(shù)矩陣行向量和觀測(cè)向量的初始權(quán)陣

（2）計(jì)算RWTLS迭代初始值:

（3）利用前一次得到的參數(shù)近似值，加上本次的ξ（i）作為新的擬合參數(shù)近似值，更新誤差方程系數(shù)矩陣B和觀測(cè)向量L，依據(jù)2.2節(jié)中權(quán)陣定義，更新權(quán)陣和

（4）迭代更新相關(guān)值:

（5）計(jì)算ξ（i+1）:

（6）由式（15）計(jì)算點(diǎn)到圓柱中軸線的距離與圓柱半徑的較差d i、均值和標(biāo)準(zhǔn)差:

此時(shí)，圓柱面擬合單位權(quán)中誤差估值和擬合精度為:

3 實(shí)例分析

利用Leica ScanStation C10三維激光掃描儀對(duì)蘇州市太平天國(guó)忠王府某圓形立柱進(jìn)行掃描，經(jīng)配準(zhǔn)、去噪等預(yù)處理后，得到完整的圓柱面點(diǎn)云。通過(guò)簡(jiǎn)化降低點(diǎn)云數(shù)據(jù)量，獲得實(shí)驗(yàn)的原始點(diǎn)云數(shù)據(jù)（圖1）。分別利用LS、TLS、WTLS、RWTLS算法對(duì)圓柱面點(diǎn)云數(shù)據(jù)擬合，計(jì)算圓柱面擬合精度（表1）。

圖1 圓柱面點(diǎn)云數(shù)據(jù)分布圖

表1 點(diǎn)云數(shù)據(jù)擬合精度結(jié)果

由表1可知，幾種算法得到的最終擬合參數(shù)結(jié)果較為接近。RWTLS算法的單位權(quán)中誤差明顯優(yōu)于LS、TLS和WLTS算法，分別相對(duì)提高了45%、45%和27%。在圓柱面擬合精度方面，RWTLS算法比WTLS、TLS和LS算法精度更高，擬合效果最好。圓柱面擬合時(shí)，迭代過(guò)程直接解算的是參數(shù)向量改正值，數(shù)值上非常小，所以TLS算法與LS算法的單位權(quán)中誤差幾乎一致。WTLS算法顧及距離自適應(yīng)調(diào)整權(quán)陣，具備一定抗差能力，單位權(quán)中誤差優(yōu)于TLS算法。由于點(diǎn)云數(shù)據(jù)中粗差的存在，自適應(yīng)權(quán)陣無(wú)法反映在計(jì)算擬合精度中，使得WTLS解算的擬合參數(shù)精度反而低于TLS算法。而RWTLS算法除了能自適應(yīng)修正權(quán)值，還能探測(cè)剔除粗差，可獲得更可靠的參數(shù)解，其擬合得到圓柱面的標(biāo)準(zhǔn)方程為x2+y2=0.023 313 7。

為說(shuō)明本文算法的正確性，利用田曉等[5]的全站儀觀測(cè)數(shù)據(jù)進(jìn)行驗(yàn)證，得到ΔX=107.176 75 m，ΔY=225.396 40 m，εX=0°7'16.597 67″，εX=0°5'28.252 83″，r=0.124 08 m，擬合精度（即圓度的標(biāo)準(zhǔn)差[16]）為0.53 mm，遠(yuǎn)優(yōu)于田曉等得到的0.92 mm。結(jié)果表明，本算法是有效的，且優(yōu)于傳統(tǒng)算法。

4 結(jié) 語(yǔ)

（1）針對(duì)圓柱面點(diǎn)云數(shù)據(jù)的特點(diǎn)，根據(jù)圓柱面擬合的數(shù)學(xué)模型，提出了RWTLS圓柱面擬合算法，較好地顧及了不同掃描點(diǎn)精度的微小差異，削弱了異常值對(duì)圓柱面點(diǎn)云數(shù)據(jù)擬合的干擾，且可任意選取參數(shù)初始值。

（2）經(jīng)過(guò)工程實(shí)例數(shù)據(jù)計(jì)算，證明了在圓柱面點(diǎn)云數(shù)據(jù)擬合中RWTLS算法能有效避免粗差影響，得到準(zhǔn)確可靠的擬合參數(shù)，擬合效果優(yōu)于傳統(tǒng)算法。盡管實(shí)例的圓形柱體近似垂直于地面，但算法亦適用于大旋轉(zhuǎn)角情形的柱體擬合，具有較強(qiáng)的工程應(yīng)用價(jià)值。

（3）RWTLS解算的迭代過(guò)程較為繁瑣，且收斂速度較慢，如何優(yōu)化先驗(yàn)定權(quán)的規(guī)則提高運(yùn)算效率有待進(jìn)一步研究。