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

蒼桂華 李明峰 岳建平

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

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

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

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

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

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

顧及入射角對點(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)陣。通過3種不同材質(zhì)平面樣本數(shù)據(jù)的計算分析表明，與最小二乘法、總體最小二乘法相比，該方法可得到更高的平面擬合精度。

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

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

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

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

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

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

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

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

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

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

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

3)eA≠0，QY、QA均不為單位陣。此時，假設(shè)點(diǎn)云數(shù)據(jù)坐標(biāo)為不等精度，采用加權(quán)總體最小二乘法進(jìn)行參數(shù)估計。根據(jù)系數(shù)矩陣的結(jié)構(gòu)特點(diǎn)設(shè)計權(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 種估算方法的估計準(zhǔn)則和單位權(quán)方差估計如表1所示。

表1 3種方法的估算準(zhǔn)則和單位權(quán)方差估計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)行平面擬合時應(yīng)根據(jù)各點(diǎn)精度確定其擬合權(quán)重。點(diǎn)位精度越高，參與擬合的權(quán)重也越大，從而保證結(jié)果的精度。相關(guān)文獻(xiàn)研究了入射角對點(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三個方向等精度獲取，對于平面的系數(shù)陣列向量和觀測向量，有σx=σy=σz。結(jié)合系數(shù)矩陣A的特點(diǎn)，設(shè)置相應(yīng)權(quán)陣:

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

2.2 方法步驟

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

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

5)重新計算每點(diǎn)的入射角余弦值，設(shè)置相關(guān)權(quán)陣 P0、PX、PY，計算(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ì)性和平面度最差，但其平面偏差目測在1 mm以內(nèi)。利用徠卡HDS ScanStation2分別對樣本進(jìn)行掃描，獲取平面點(diǎn)云數(shù)據(jù)(圖1)。

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

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

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

從表2可以看出，不論何種材質(zhì)樣本，LS方法的單位權(quán)中誤差都遠(yuǎn)遠(yuǎn)大于TLS方法和AWTLS方法，而TLS和AWTLS方法的單位權(quán)中誤差為同一數(shù)量級，這是因?yàn)長S方法與TLS方法以及AWTLS方法的單位權(quán)中誤差計算方法不同所致。另外，對于3種材質(zhì)的樣本數(shù)據(jù)，利用AWTLS方法獲得的3個精度判定指標(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)計算結(jié)果Tab.2 Fitting results of plane samples

4 結(jié)論

1)AWTLS方法是在加權(quán)總體最小二乘方法基礎(chǔ)上，利用入射角對點(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)中誤差計算方法不同，因此不能簡單地以單位權(quán)中誤差的大小來判定某個估算方法的優(yōu)劣，應(yīng)結(jié)合其他精度評判指標(biāo)進(jìn)行評定。

3)AWTLS方法系數(shù)矩陣的行向量權(quán)陣PX和觀測向量權(quán)矩陣PY與擬合參數(shù)有關(guān)，因此在計算過程中需要迭代運(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)寧，等.測量數(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é)報:信息科學(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 王樂洋，許才軍，魯鐵定.病態(tài)加權(quán)總體最小二乘平差的嶺估計解法［J］.武漢大學(xué)學(xué)報:信息科學(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］.大地測量與地球動力學(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］.測繪學(xué)報，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］.大地測量學(xué)與地球動力學(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è)測量擬合［M］.北京:測繪出版社，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

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

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