蒼桂華 李明峰 岳建平

1)南京工業(yè)大學(xué)測(cè)繪學(xué)院，南京 210009

2)河海大學(xué)地球科學(xué)與工程學(xué)院，南京 210098

以入射角定權(quán)的點(diǎn)云數(shù)據(jù)加權(quán)總體最小二乘平面擬合研究*

蒼桂華1)李明峰1)岳建平2)

1)南京工業(yè)大學(xué)測(cè)繪學(xué)院，南京 210009

2)河海大學(xué)地球科學(xué)與工程學(xué)院，南京 210098

顧及入射角對(duì)點(diǎn)云數(shù)據(jù)點(diǎn)位精度的影響，提出入射角定權(quán)的點(diǎn)云數(shù)據(jù)加權(quán)總體最小二乘平面擬合方法。該方法以加權(quán)總體最小二乘方法為基礎(chǔ)，利用入射角信息確定點(diǎn)云數(shù)據(jù)中各點(diǎn)的平面擬合權(quán)重，并根據(jù)系數(shù)矩陣結(jié)構(gòu)特點(diǎn)定義其權(quán)陣。通過(guò)3種不同材質(zhì)平面樣本數(shù)據(jù)的計(jì)算分析表明，與最小二乘法、總體最小二乘法相比，該方法可得到更高的平面擬合精度。

點(diǎn)云數(shù)據(jù);平面擬合;加權(quán)總體最小二乘;入射角定權(quán);系數(shù)矩陣

地面三維激光現(xiàn)實(shí)掃描場(chǎng)景中含有大量平面特征，可用于前期數(shù)據(jù)預(yù)處理或后期建模的數(shù)據(jù)簡(jiǎn)化，還可用于多視點(diǎn)云數(shù)據(jù)配準(zhǔn)，實(shí)現(xiàn)建模。平面特征的提取是利用點(diǎn)云數(shù)據(jù)進(jìn)行3D建模的一項(xiàng)重要工作，提取精度直接關(guān)系到后續(xù)應(yīng)用成果的精度。以往提取平面特征常用的方法為最小二乘平面擬合法。此方法根據(jù)點(diǎn)的三維坐標(biāo)(xi，yi，zi)建立平面方程z=ax+by+c，假設(shè)誤差e只存在于觀測(cè)向量Z 中，建立高斯-馬爾科夫(Gauss-Markov，G-M)模型，采用最小二乘(least squares，LS)方法擬合平面，獲取平面擬合參數(shù)。然而，由于儀器設(shè)備、外界環(huán)境、地物特性等因素影響，使得實(shí)際點(diǎn)云數(shù)據(jù)的xi、yi、zi三個(gè)觀測(cè)值中均存在誤差，從而使包含變量xi、yi的系數(shù)矩陣A也含有誤差。因此，利用最小二乘法進(jìn)行平面數(shù)據(jù)擬合的結(jié)果不是最優(yōu)，而是有偏的［1］。針對(duì)這種觀測(cè)向量和系數(shù)矩陣均包含誤差的模型，即所謂的變量中的誤差(error-in-variables，EIV)模型，Golub等［2］最早提出了總體最小二乘(total least squares，TLS)估計(jì)方法。但總體最小二乘估計(jì)僅在設(shè)計(jì)矩陣和殘差元素均服從獨(dú)立等精度分布時(shí)才是最優(yōu)估計(jì)［3］，而實(shí)際點(diǎn)云數(shù)據(jù)各點(diǎn)的精度是不等的，因此簡(jiǎn)單的總體最小二乘方法并非最優(yōu)估計(jì)［4］。為了解決矩陣和殘差的不等精度估計(jì)問(wèn)題，Markovsky等［5］提出了加權(quán)總體最小二乘(weighted total least squares，WTLS)方法。Schaffrin等［6］則進(jìn)一步擴(kuò)展了WTLS方法。國(guó)內(nèi)學(xué)者也從理論、應(yīng)用等方面對(duì) WTLS方法進(jìn)行了深入研究［7－10］。本文根據(jù)入射角對(duì)點(diǎn)云數(shù)據(jù)點(diǎn)位精度的影響，提出利用入射角定權(quán)的加權(quán)總體最小二乘點(diǎn)云數(shù)據(jù)平面擬合方法，簡(jiǎn)稱角度加權(quán)總體最小二乘(angle weighted total least squares，AWTLS)方法。

1 LS、TLS、WTLS 三種方法及其平差模型

點(diǎn)云數(shù)據(jù)經(jīng)過(guò)粗差剔除預(yù)處理后得到一組平面點(diǎn)云的三維坐標(biāo)(xi，yi，zi)(i=1，2，…，n)。設(shè)空間平面方程式為:

式中，a、b、c為待求的平面擬合參數(shù)。將式(1)寫(xiě)成矩陣形式為:

若同時(shí)考慮觀測(cè)向量Y和系數(shù)矩陣A的誤差，則建立EIV函數(shù)模型［6］:

式中，eY和EA分別表示觀測(cè)向量Y和系數(shù)矩陣A的隨機(jī)誤差矩陣;vec是指將矩陣按列拉直所得到的列向量，排列順序?yàn)閺淖蟮接?eA=vec(EA)∈R3n×1表示矩陣EA按列向量化后得到的向量;σ20為未知的單位權(quán)方差;QY=P－1Y，QA=P－1A，分別表示觀測(cè)向量和系數(shù)陣列向量化向量的協(xié)因素陣。

根據(jù)eA及其相關(guān)權(quán)陣的設(shè)置，采用不同的估計(jì)方法進(jìn)行平面擬合參數(shù)的求解。

1)eA≡0。此時(shí)只考慮觀測(cè)向量Y的誤差，則EIV模型變?yōu)榻?jīng)典測(cè)量平差的G-M模型:

應(yīng)采用最小二乘法進(jìn)行參數(shù)估計(jì)。

2)eA≠0，但 QY=In，QA=Inm。此時(shí)假設(shè)點(diǎn)云數(shù)據(jù)坐標(biāo)獨(dú)立且等精度，即觀測(cè)向量權(quán)陣PY和系數(shù)矩陣權(quán)陣PA均為單位陣，采用總體最小二乘法進(jìn)行參數(shù)估算。

3)eA≠0，QY、QA均不為單位陣。此時(shí)，假設(shè)點(diǎn)云數(shù)據(jù)坐標(biāo)為不等精度，采用加權(quán)總體最小二乘法進(jìn)行參數(shù)估計(jì)。根據(jù)系數(shù)矩陣的結(jié)構(gòu)特點(diǎn)設(shè)計(jì)權(quán)陣P0、Px、PY。P0是3×3矩陣，代表系數(shù)矩陣 A 的列向量權(quán)陣;PX是n×n矩陣，代表系數(shù)矩陣A的行向量權(quán)陣;PY是n×n矩陣，代表矩陣Y的權(quán)陣。P0、PX、PY相應(yīng)的協(xié)因素矩陣為 Q0、QX、QY，即

式中，?表示“kronecker積”。

3 種估算方法的估計(jì)準(zhǔn)則和單位權(quán)方差估計(jì)如表1所示。

表1 3種方法的估算準(zhǔn)則和單位權(quán)方差估計(jì)Tab.1 Estimation criterions and unit weight estimation variances with three methods

2 AWTLS方法權(quán)陣的確定及其步驟

2.1 權(quán)陣的確定

點(diǎn)云數(shù)據(jù)中各點(diǎn)的精度不等，進(jìn)行平面擬合時(shí)應(yīng)根據(jù)各點(diǎn)精度確定其擬合權(quán)重。點(diǎn)位精度越高，參與擬合的權(quán)重也越大，從而保證結(jié)果的精度。相關(guān)文獻(xiàn)研究了入射角對(duì)點(diǎn)云數(shù)據(jù)精度的影響，指出點(diǎn)云數(shù)據(jù)入射角越小則點(diǎn)位精度越高［11－12］。而入射角越小，其余弦值越大，因此可將各點(diǎn)入射角余弦值作為其擬合權(quán)重。平面點(diǎn)云數(shù)據(jù)的各點(diǎn)入射角θ定義為入射光線與平面法向量之間的夾角［13］。設(shè)擬合的平面方程如式(1)所示，則平面上i點(diǎn)的入射角余弦值或其擬合權(quán)值為:

式中，θi為 i點(diǎn)的入射角，(xi，yi，zi)為點(diǎn) i的三維坐標(biāo)，(a，b，－1)為擬合平面的法向量。

設(shè)點(diǎn)云在x、y、z三個(gè)方向等精度獲取，對(duì)于平面的系數(shù)陣列向量和觀測(cè)向量，有σx=σy=σz。結(jié)合系數(shù)矩陣A的特點(diǎn)，設(shè)置相應(yīng)權(quán)陣:

式中，P0的第三個(gè)對(duì)角元素為0，表示系數(shù)矩陣A中第三列不需要改正;其余對(duì)角線元素為1，表示系數(shù)矩陣A中的第一、二列中的元素是等精度獲取的;PX、PY與入射角有關(guān)。

2.2 方法步驟

由式(9)可以看出，權(quán)陣PX、PY與平面擬合參數(shù)相關(guān)，因此需要迭代求解。AWTLS方法的具體步驟如下:

1)將最小二乘平面擬合參數(shù)作為起始解，根據(jù)式(8)計(jì)算各點(diǎn)的入射角余弦值，組成相關(guān)矩陣P0、PX、PY。

5)重新計(jì)算每點(diǎn)的入射角余弦值，設(shè)置相關(guān)權(quán)陣 P0、PX、PY，計(jì)算(i+1):

3 算例分析

為驗(yàn)證AWTLS方法的可行性，采用3種不同材質(zhì)的平面點(diǎn)云樣本數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)的平面樣本有標(biāo)準(zhǔn)反射板(反射率為35%)、普通木板及一般的水泥建材模板。3種材質(zhì)樣本中，標(biāo)準(zhǔn)反射板均質(zhì)性最好，板面最平整;普通木板均質(zhì)性和平面度次之;建材水泥模板均質(zhì)性和平面度最差，但其平面偏差目測(cè)在1 mm以內(nèi)。利用徠卡HDS ScanStation2分別對(duì)樣本進(jìn)行掃描，獲取平面點(diǎn)云數(shù)據(jù)(圖1)。

根據(jù)點(diǎn)云數(shù)據(jù)特點(diǎn)確定平面形式［14］。分別利用最小二乘法、總體最小二乘法和AWTLS法對(duì)各樣本點(diǎn)云數(shù)據(jù)進(jìn)行平面擬合，獲得平面擬合參數(shù)、、以及單位權(quán)中誤差(0)。假設(shè)確定的平面形式如式(1)所示，由式(15)計(jì)算點(diǎn)云數(shù)據(jù)各點(diǎn)到擬合面的距離(di)以及平面擬合精度(P)，并得到點(diǎn)到擬合面的最大距離(max(di)):

式中，n為擬合平面上點(diǎn)的個(gè)數(shù)。

將單位權(quán)中誤差、點(diǎn)到擬合面的最大距離、平面擬合精度作為估算方法優(yōu)劣的評(píng)判指標(biāo)。用3種方法對(duì)3種材質(zhì)的樣本數(shù)據(jù)進(jìn)行平面擬合，結(jié)果如表2所示。

從表2可以看出，不論何種材質(zhì)樣本，LS方法的單位權(quán)中誤差都遠(yuǎn)遠(yuǎn)大于TLS方法和AWTLS方法，而TLS和AWTLS方法的單位權(quán)中誤差為同一數(shù)量級(jí)，這是因?yàn)長(zhǎng)S方法與TLS方法以及AWTLS方法的單位權(quán)中誤差計(jì)算方法不同所致。另外，對(duì)于3種材質(zhì)的樣本數(shù)據(jù)，利用AWTLS方法獲得的3個(gè)精度判定指標(biāo)均要比LS方法和TLS方法小得多。以平面度和均質(zhì)性最差的水泥模板為例，利用AWTLS方法得到的單位權(quán)中誤差比TLS方法提高了65%，平面擬合精度比LS方法和TLS方法分別提高了63%和19%，點(diǎn)到擬合面的最大距離也由LS方法和TLS方法的9.1 mm和7.0 mm 降為5.3 mm。

圖1 實(shí)驗(yàn)樣本的平面點(diǎn)云數(shù)據(jù)Fig.1 Plane point clouds used in the experiments

表2 三種方法下平面點(diǎn)云數(shù)據(jù)擬合的相關(guān)計(jì)算結(jié)果Tab.2 Fitting results of plane samples

4 結(jié)論

1)AWTLS方法是在加權(quán)總體最小二乘方法基礎(chǔ)上，利用入射角對(duì)點(diǎn)云數(shù)據(jù)精度的影響確定各點(diǎn)的擬合權(quán)重。從實(shí)驗(yàn)結(jié)果看，該定權(quán)方法符合點(diǎn)云數(shù)據(jù)平面擬合的實(shí)際，結(jié)果也更加合理。

2)由于LS方法、TLS方法以及AWTLS方法的單位權(quán)中誤差計(jì)算方法不同，因此不能簡(jiǎn)單地以單位權(quán)中誤差的大小來(lái)判定某個(gè)估算方法的優(yōu)劣，應(yīng)結(jié)合其他精度評(píng)判指標(biāo)進(jìn)行評(píng)定。

3)AWTLS方法系數(shù)矩陣的行向量權(quán)陣PX和觀測(cè)向量權(quán)矩陣PY與擬合參數(shù)有關(guān)，因此在計(jì)算過(guò)程中需要迭代運(yùn)算。

利用WTLS方法進(jìn)行點(diǎn)云數(shù)據(jù)平面擬合的關(guān)鍵是如何確定各點(diǎn)的擬合權(quán)值。本文只進(jìn)行了利用入射角定權(quán)的WTLS方法在點(diǎn)云數(shù)據(jù)平面擬合中的適用性研究，后續(xù)將進(jìn)一步研究其他權(quán)值確定方法，并與本文的AWTLS方法進(jìn)行比較，研究各方法的優(yōu)劣。

1 邱衛(wèi)寧，等.測(cè)量數(shù)據(jù)處理理論與方法［M］.武漢:武漢大學(xué)出版社，2008.(Qiu Weining，et al.The theory and method of surveying data processing［M］.Wuhan:Wuhan University Press，2008)

2 Golub H G，Loan C F.An analysis of the total least squares problem［J］.SIAM Journal on Numerical Analysis，1980，17(6):883－893.

3 Van H S，Vandewalle J.The total least squares problem:computational aspects and analysis［M］.Philadelphia:SIAM，1991.

4 周擁軍，鄧才華.加權(quán)和不加權(quán)TLS方法及其在不等精度坐標(biāo)變換中的應(yīng)用［J］.武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2012，37(8):976 － 979.(Zhou Yongjun，Deng Caihua.Weighted and unweighted total least square methods and applications to heteroscedastic 3d coordinate transformation［J］.Geomatics and Information Science of Wuhan University，2012，37(8):976 －979)

5 Markovsky I，et al.The element-wise weighted total leastsquares problem ［J］.Computational Statistics and Data A-nalysis，2006，50:181 －209.

6 Schaffrin B，Wieser A.On weighted total least-square adjustment for linear regression［J］.Journal of Geodesy，2008，82(7):415－421.

7 王樂(lè)洋，許才軍，魯鐵定.病態(tài)加權(quán)總體最小二乘平差的嶺估計(jì)解法［J］.武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2010，35(11):1 346 － 1 350.(Wang Leyang，Xu Caijun，Lu Tieding.Ridge estimation method in ill-posed weighted total least squares adjustment［J］.Geomatics and Information Science of Wuhan University，2010，35(11):1 346 －1 350)

8 趙輝.基于加權(quán)總體最小二乘的GPS數(shù)據(jù)高程擬合［J］.大地測(cè)量與地球動(dòng)力學(xué)，2011，31(5):88 －90，96.(Zhao Hui.GPS height fitting of weighted total least-squares adjustment［J］.Journal of Geodesy and Geodynamics，2011，31(5):88 －90，96)

9 周擁軍，朱建軍，鄧才華.附參數(shù)的條件平差與按行獨(dú)立的加權(quán)總體最小二乘法估計(jì)的一致性研究［J］.測(cè)繪學(xué)報(bào)，2012，41(1):48 － 53.(Zhou Yongjun，Zhu Jianjun，Deng Caihua.The consistency between row-wised weighted total least squares and condition adjustment with parameters［J］.Acta Geodaetica et Cartographica Sinica，2012，41(1):48－53)

10 陳瑋嫻，等.加權(quán)總體最小二乘在三維激光標(biāo)靶擬合中的應(yīng)用［J］.大地測(cè)量學(xué)與地球動(dòng)力學(xué)，2010，30(5):90－96.(Chen Weixian，et al.Application of weighted total least squares to target fitting of three-dimensional laser scanning［J］.Journal of Geodesy and Geodynamics，2010，30(5):90－96)

11 Soudarissanane S，Van R J，Bucksch A，et al.Error budget of terrestrial laser scanning:influence of the incidence angle on the scan quality［C］.3D-Nord Ost，Berlin，2007.

12 Soudarissanane S，et al.Incidence angle influence on the quality of terrestrial laser scanning points［C］.ISPRS Workshop Laserscanning，2009.

13 Soudarissanane S，et al.Scanning geometry:influence factor on the quality of terrestrial laser scanning points［J］.ISPRS Journal of Photogrammetry and Remote Sensing，2011，66(4):389－399.

14 王解先，季凱敏.工業(yè)測(cè)量擬合［M］.北京:測(cè)繪出版社，2008.(Wang Jiexian，Ji Kaimin.Industrial measurement fitting［M］.Beijing:Surveying and Mapping Press，2008)

STUDY ON POINT CLOUDS PLANE FITTING WITH WEIGHTED TOTAL LEAST SQUARES BASED ON INCIDENCE ANGLE WEIGHTING

Cang Guihua1)，Li Mingfeng1)and Yue Jianping2)
1)Department of Geomatics Engineering，Nanjing University of Technology，Nanjing 210009
2)School of Earth Sciences and Engineering，Hohai University，Nanjing210098

The incidence angle has an impact on point clouds accuracyin 3d terrestrial laser scanning.For improving the accuracy，a method of weighted total least squares(WTLS)incidence angle weighting for point clouds plane fitting is proposed.In the method，the weight value of each point is determined by incidence angle，and the weight matrix is defineson the basis of the structure characteristics of coefficient matrix.The experiments with point clouds data of the plane made of three different materials show that the method is more reasonable and more accurate than least squares(LS)and total least squares(TLS).

point clouds;plane fitting;weighted total least squares(WTLS);incidence angle;coefficient matrix

P207

A

1671-5942(2014)03-0095-04

2013-09-22

國(guó)家自然科學(xué)基金項(xiàng)目(41174002);江蘇省測(cè)繪科研項(xiàng)目(JSCHKY201303)。

蒼桂華，女，1971年生，講師，主要研究方向?yàn)閿z影測(cè)量與遙感、大地測(cè)量與測(cè)繪工程數(shù)據(jù)處理。E－mail:ghuac@163.com。