穩(wěn)健估計(jì)方法在點(diǎn)云圓柱擬合中的應(yīng)用

劉東林王解先

1同濟(jì)大學(xué)測(cè)繪與地理信息學(xué)院，上海，200092

與傳統(tǒng)的全站儀測(cè)量相比，三維激光掃描技術(shù)的高效、快捷、精確、方便等優(yōu)點(diǎn)，使其在城市三維建模、礦山地質(zhì)測(cè)繪、變形監(jiān)測(cè)和古跡測(cè)量等領(lǐng)域具有廣泛的應(yīng)用［1］。在工業(yè)測(cè)量中，需要根據(jù)實(shí)際情況擬合各種圓柱形物體，進(jìn)行圓柱擬合的目的就是計(jì)算得到圓柱的圓心坐標(biāo)、方向和半徑［2］。

目前，常用的圓柱擬合方法中：高斯圖法、遺傳算法和特征值法這3種方法的理論模型較為復(fù)雜，不易理解和實(shí)現(xiàn)，由于參數(shù)之間的相關(guān)性，需要選取較為準(zhǔn)確的初值且不具備抗差性；幾何特性法、圓度判別法和坐標(biāo)轉(zhuǎn)換法這3種方法較容易理解和實(shí)現(xiàn)，但是對(duì)初值也有較強(qiáng)的依賴性，同時(shí)也不具備抗差性［3-5］。文獻(xiàn)［6］通過變換基于坐標(biāo)轉(zhuǎn)換的誤差方程來降低參數(shù)之間的相關(guān)性。文獻(xiàn)［7］主要利用三倍中誤差探測(cè)法剔除粗差后實(shí)現(xiàn)了坐標(biāo)轉(zhuǎn)換法擬合圓柱。文獻(xiàn)［8］主要介紹了再生權(quán)最小二乘法（selfborn weighted least squares，SBWLS）和其他穩(wěn)健估計(jì)方法的對(duì)比分析。本文將介紹一種有效可行的圓柱擬合方法，該方法基于坐標(biāo)轉(zhuǎn)換擬合法和再生權(quán)最小二乘穩(wěn)健估計(jì)法，在降低參數(shù)對(duì)初值選取的要求的同時(shí)，也不用對(duì)觀測(cè)值進(jìn)行粗差剔除處理，還可以擬合任意傾斜角度的圓柱。

1 擬合模型與參數(shù)估計(jì)方法

1.1 坐標(biāo)轉(zhuǎn)換擬合法獲取初值

為了便于求解圓柱的幾何參數(shù)，以圓柱中心軸線與測(cè)量坐標(biāo)系XOY平面的交點(diǎn)為坐標(biāo)原點(diǎn)，圓柱中軸線為Z軸，建立圓柱標(biāo)準(zhǔn)坐標(biāo)系，將測(cè)量坐標(biāo)系X通過坐標(biāo)轉(zhuǎn)換到標(biāo)準(zhǔn)坐標(biāo)系Y［4］：

式中，X0=(x0y0z0)T為平移量；R=R1(α)R2(β)R3(γ)為旋轉(zhuǎn)矩陣；α、β、γ分別是繞3個(gè)坐標(biāo)軸進(jìn)行旋轉(zhuǎn)的角度。由標(biāo)準(zhǔn)坐標(biāo)系的定義可知，標(biāo)準(zhǔn)坐標(biāo)系與測(cè)量坐標(biāo)系的z軸平行，因此有γ=0且z0可以任意選取，本文取z0=0。這樣，只需要求取x0、y0、α、β和圓柱半徑r這5個(gè)參數(shù)，就可以唯一確定圓柱的姿態(tài)和大小。

標(biāo)準(zhǔn)坐標(biāo)系中圓柱面可表示為：

得到以測(cè)量點(diǎn)到圓柱表面的距離為殘差的誤差方程：

式中，x和y是標(biāo)準(zhǔn)坐標(biāo)系中的點(diǎn)坐標(biāo)；r為圓柱半徑。

將誤差方程轉(zhuǎn)換為矩陣的形式：

式中，Λ=diag(1 1 0)。且其中z0=0，γ=0為已知值，不參與參數(shù)求解。

將誤差方程線性化，由于存在根式，導(dǎo)致各參數(shù)之間的相關(guān)性很強(qiáng)，所以線性化過程中所舍棄的二次項(xiàng)不可忽略，因此需要較為準(zhǔn)確的初值，若參數(shù)初值與真值偏差較大，在參數(shù)求解時(shí)很難拿收斂到正確的值，從而得不到正確的結(jié)果［3］。

為了降低各參數(shù)之間相關(guān)性，對(duì)誤差方程進(jìn)行變換。令誤差方程為：

將誤差方程轉(zhuǎn)換為矩陣的形式：

與誤差方程式（4）相比，方程式（6）中沒有根式，故參數(shù)之間的相關(guān)性大大降低，這樣線性化過程中的二次項(xiàng)就可以舍去，降低了對(duì)參數(shù)初值選取的要求，避免參數(shù)求解錯(cuò)誤。

根據(jù)間接平差迭代求取參數(shù)估值，但是該方法不是以圓度平方和最小為準(zhǔn)則［7］，得到的結(jié)果并不是參數(shù)的最佳估值。但此時(shí)參數(shù)估值已經(jīng)較為準(zhǔn)確，可以將其作為參數(shù)初值按式（5）的誤差方程再次迭代求解，從而得到參數(shù)的最佳估值?？捎墒剑?）得到圓柱面方程，由式（5）得到的殘差值即為圓柱的圓度值。

1.2 再生權(quán)最小二乘穩(wěn)健估計(jì)法

統(tǒng)計(jì)學(xué)家指出，觀測(cè)值中出現(xiàn)粗差的概率大概為1%～10%，這些粗差會(huì)降低參數(shù)估計(jì)的精度，因此，在參數(shù)估計(jì)的過程中，需要消除或減弱粗差對(duì)參數(shù)估計(jì)的影響。目前，對(duì)含粗差觀測(cè)值的處理方法主要有兩種，第一種是將含粗差觀測(cè)值視為期望異常，用統(tǒng)計(jì)檢驗(yàn)法剔除含粗差的觀測(cè)值之后再進(jìn)行參數(shù)估計(jì)；第二種是將含粗差觀測(cè)值視為方差異常，采用穩(wěn)健估計(jì)方法對(duì)其進(jìn)行處理［8］。通過三維激光掃描獲取的點(diǎn)云數(shù)據(jù)中，包含了很多異常數(shù)據(jù)，根據(jù)點(diǎn)云數(shù)據(jù)對(duì)圓柱進(jìn)行擬合時(shí)，在不進(jìn)行粗差剔除的情況下擬合得到較為精確的圓柱參數(shù)。本文采用了穩(wěn)健估計(jì)方法中的SBWLS來解決這個(gè)問題。

SBWLS是一種有效可行的穩(wěn)健估計(jì)方法。它充分利用觀測(cè)值的改正數(shù)（殘差）條件方程提供的有效信息，用觀測(cè)值改正數(shù)的多個(gè)估值來構(gòu)造觀測(cè)值的權(quán)，而不是直接用最小二乘法（least squares，LS）得到的觀測(cè)值改正數(shù)來計(jì)算觀測(cè)值的權(quán)。實(shí)驗(yàn)證明，與其他常用的穩(wěn)健估計(jì)方法相比，SBWLS法能更有效地消除或減弱粗差對(duì)參數(shù)估計(jì)的影響［8］。

列出誤差方程：

式中，l=-(d+BX0-L)。

將式（8）中的系數(shù)矩陣B進(jìn)行行變換，將其分為兩部分，其中Bt為t×t階的滿秩方陣，對(duì)V和l進(jìn)行同步變換，得到Vt和lt：

由式（9）可得

用同一個(gè)觀測(cè)值改正數(shù)（殘差）的多個(gè)估值，由多個(gè)估值計(jì)算觀測(cè)值的再生方差，根據(jù)再生方差計(jì)算得到觀測(cè)值的再生權(quán)［9-11］。其計(jì)算公式為：

將計(jì)算得到的再生權(quán)作為觀測(cè)值的權(quán)，按LS法迭代求解參數(shù)估值的方法稱為SBWLS9］。

2 圓柱擬合計(jì)算實(shí)例

2.1 數(shù)據(jù)準(zhǔn)備

掃描了某圓柱形物體，得到如表1所示的點(diǎn)云數(shù)據(jù)（人為的加入了一些粗差點(diǎn)）。點(diǎn)云三維散點(diǎn)圖如圖1所示，圖1中散布于集中點(diǎn)云周圍的點(diǎn)，即為粗差點(diǎn)。

表1 點(diǎn)云坐標(biāo)數(shù)據(jù)Tab.1 Point Cloud Data

圖1 點(diǎn)云數(shù)據(jù)散點(diǎn)圖Fig.1 Coordinate 3D Scatter Images

2.2 圓柱擬合

基于上述方法擬合此點(diǎn)云柱面參數(shù)，將五參數(shù)初值設(shè)置為x0=0，y0=0，α=0，β=0，r=1。首先，使用初值任意選取的方法迭代求取參數(shù)，將此作為初值，分別使用LS（不抗差）和SBWLS（抗差）進(jìn)行參數(shù)估計(jì)，比較分析兩種方法的平差結(jié)果。

從表2（表中參數(shù)含義見§1.1）可以得出，LS和SBWLS-1兩種方法平差參數(shù)結(jié)果接近，SBWLS-1法平差的單位權(quán)中誤差小于LS法，證明了SBWLS法的抗差性。圖2為圓柱擬合結(jié)果圖。

圖2 圓柱三維擬合圖Fig.2 Cylinder 3D Fitting Images

為了驗(yàn)證該方法的通用性，將表1中的點(diǎn)云數(shù)據(jù)進(jìn)行一定數(shù)值的坐標(biāo)平移和旋轉(zhuǎn)變換，然后按照SBWLS法求解圓柱面參數(shù)，由表2中的SBWLS-2可得出，擬合結(jié)果中圓柱的半徑r不變，說明該方法能成功擬合各種情況下的圓柱。

為了進(jìn)一步驗(yàn)證該方法的正確性，本文采用文獻(xiàn)［6］中的實(shí)例數(shù)據(jù)（不含粗差），用SBWLS法的擬合結(jié)果與文獻(xiàn)［6］的擬合結(jié)果進(jìn)行比較。如表3所示（表中參數(shù)含義見§1.1），可得出兩種方法的擬合結(jié)果基本一致，擬合精度也相當(dāng)。

表3 兩種方法結(jié)果對(duì)比Tab.3 Comparison of Results of Two Methods

2.3 擬合結(jié)果評(píng)定

一般可以用中誤差、殘差、殘差分布特性等來評(píng)定圓柱面擬合結(jié)果［6］。SBWLS法計(jì)算得到的單位權(quán)中誤差為0.524 mm，與采集點(diǎn)云數(shù)據(jù)儀器的精度相吻合，說明此方法的擬合結(jié)果良好。另外，通過對(duì)計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)分析，繪制殘差的折線圖和分布圖，如圖3（a）和圖4（a）所示，發(fā)現(xiàn)粗差點(diǎn)的殘差非常大。對(duì)非粗差點(diǎn)的殘差進(jìn)行局部放大，如圖3（b）和圖4（b）所示，發(fā)現(xiàn)其分布呈現(xiàn)隨機(jī)性和正態(tài)性，符合偶然誤差的幾大特性。由此可以得出SBWLS法的擬合結(jié)果具有較高可靠性。

圖3 殘差（圓度）折線圖Fig.3 Residual Line Images

圖4 殘差（圓度）分布圖Fig.4 Residual Plot Images

3 結(jié)束語

在圓柱擬合的實(shí)際工程應(yīng)用中，很多情況下無法知道圓柱參數(shù)的近似值，導(dǎo)致誤差方程結(jié)構(gòu)不夠理想，同時(shí)觀測(cè)值中不可避免的存在一些粗差，因此僅僅依靠最小二乘法求解出正確的圓柱參數(shù)就比較困難。本文方法無需顧及初值選取、粗差探測(cè)。通過對(duì)含有粗差的不同方向的點(diǎn)云數(shù)據(jù)進(jìn)行擬合，發(fā)現(xiàn)能夠有效地消除或減弱粗差對(duì)參數(shù)估計(jì)的影響；通過與其他圓柱面擬合的方法比較，發(fā)現(xiàn)該方法所求解的圓柱參數(shù)結(jié)果和精度符合實(shí)際，適用于各種情況下的圓柱擬合；通過進(jìn)一步對(duì)參數(shù)結(jié)果進(jìn)行統(tǒng)計(jì)分析，發(fā)現(xiàn)圓度的分布隨機(jī)且符合正態(tài)特性，驗(yàn)證了該方法的正確性。