加權(quán)總體最小二乘點(diǎn)云平面擬合定權(quán)方法探討

李明峰 歐江霞 檀 丁 楊輝輝

1 南京工業(yè)大學(xué)地球空間信息研究中心，南京市浦珠南路30號(hào)，210009

2 武漢大學(xué)測(cè)繪學(xué)院，武漢市珞喻路129號(hào)，430079

受觀測(cè)環(huán)境、系統(tǒng)誤差等因素影響，三維激光掃描儀獲取的平面三維坐標(biāo)點(diǎn)云數(shù)據(jù)x、y、z三個(gè)方向上均含有誤差且點(diǎn)位精度各不相同，若采用最小二乘（least squares，LS）法及總體最小二乘（total least squares，TLS）法［1－3］對(duì)平面點(diǎn)云數(shù)據(jù)進(jìn)行擬合，所得參數(shù)解并非平面參數(shù)的最或然值。此時(shí)，可基于各點(diǎn)的點(diǎn)位精度確定其相對(duì)權(quán)值，引入加權(quán)總體最小二乘（weighted total least squares，WTLS）法［4－6］解算平面參數(shù)。在參數(shù)迭代解算過程中，權(quán)陣的設(shè)定將對(duì)解算結(jié)果產(chǎn)生較大影響。而由于缺少統(tǒng)一的定權(quán)準(zhǔn)則，無法確保解算結(jié)果的精度，其可信度較低。對(duì)此，本文在原有以強(qiáng)度值定權(quán)及后驗(yàn)入射角定權(quán)的基礎(chǔ)上，提出以先驗(yàn)入射角定權(quán)及距離定權(quán)兩種新的定權(quán)方法，同時(shí)結(jié)合穩(wěn)健估計(jì)，分別以強(qiáng)度值標(biāo)準(zhǔn)差、距離標(biāo)準(zhǔn)差為依據(jù)，設(shè)定相應(yīng)閾值，在剔除異常數(shù)據(jù)的基礎(chǔ)上，構(gòu)造了以強(qiáng)度值定權(quán)的穩(wěn)健加權(quán)總體最小二乘（robust weighted total least squares based on intensity，RWTLSI）法、以先驗(yàn)入射角定權(quán)的加權(quán)總體最小二乘（weighted total least squares based on prior incident angle，WTLSPIA－1）法、以距離定權(quán)的加權(quán)總體最小二乘（weighted total least squares based on distance，WTLSD）法、以強(qiáng)度值定權(quán)的加權(quán)最小二乘（weighted total least squares based on intensity，WTLSI）法、以后驗(yàn)入射角定權(quán)的加權(quán)最小二乘（weighted total least squares based on posteriori incident angle，WTLSPIA－2）法［8－9］，應(yīng)用于擬合掃描不同反射材質(zhì)獲取的平面點(diǎn)云數(shù)據(jù)，驗(yàn)證各WTLS法的有效性及適用性。

1 加權(quán)總體最小二乘基本原理

點(diǎn)云數(shù)據(jù)平面擬合函數(shù)模型為：

式中，a、b、c為待求平面擬合參數(shù)。顧及觀測(cè)向量誤差與系數(shù)矩陣誤差的EIV（errors－in－variable）模型為：

式中，Z為含有隨機(jī)誤差ez的n×1維觀測(cè)向量，A是含有隨機(jī)誤差EA的n×m維系數(shù)矩陣，ξ為待估參數(shù)：

隨機(jī)誤差ez與EA的統(tǒng)計(jì)性質(zhì)如下：

式中，“?”為Kronecker積，vec（）為矩陣?yán)弊儞Q，為未知方差分量，QZ、QA為ez與eA的對(duì)稱、非奇異協(xié)因數(shù)陣，且有：

式中，PZ為觀測(cè)值權(quán)陣，PA為系數(shù)陣A的權(quán)陣，P0、PXY分別為系數(shù)陣A的列向量權(quán)陣及行向量權(quán)陣。

加權(quán)總體最小二乘的參數(shù)估計(jì)準(zhǔn)則為：

2 定權(quán)方法探討

2.1 強(qiáng)度值定權(quán)

對(duì)于相同材質(zhì)的反射面，激光反射強(qiáng)度與點(diǎn)位精度相對(duì)應(yīng)，即點(diǎn)位的激光反射強(qiáng)度越大，點(diǎn)位精度越高，其擬合權(quán)重越大；反之，反射強(qiáng)度越小，點(diǎn)位精度越低，其擬合權(quán)重越小。對(duì)于范圍在［－2 047，2 048］的強(qiáng)度值，可利用式（6）轉(zhuǎn)換成［0，1］之間的數(shù)值，即可確定為該點(diǎn)的權(quán)值［6－7］：

以強(qiáng)度值定權(quán)可構(gòu)造RWTLSI法。本文顧及異常強(qiáng)度值對(duì)擬合精度的影響，通過式（7）計(jì)算觀測(cè)數(shù)據(jù)的強(qiáng)度值標(biāo)準(zhǔn)差，剔除的點(diǎn)，提高擬合數(shù)據(jù)精度，之后以強(qiáng)度值定權(quán)，構(gòu)造了RWTLSI法：

式中，Ii為點(diǎn)位強(qiáng)度值，為平均強(qiáng)度值，σI為強(qiáng)度值標(biāo)準(zhǔn)差。

2.2 入射角定權(quán)

激光反射強(qiáng)度與入射角存在如下關(guān)系：

式中，I′i表示入射角為θi時(shí)的點(diǎn)位強(qiáng)度值，I0為垂直入射時(shí)的點(diǎn)位強(qiáng)度值。由此可知，點(diǎn)位精度與入射角大小成反比，即點(diǎn)位入射角越小，該點(diǎn)的點(diǎn)位精度越高；反之，精度越低。因此，可將點(diǎn)位的入射角余弦值作為該點(diǎn)的權(quán)值［10－11］：

式中，（xi，yi，zi）為點(diǎn)i的三維坐標(biāo)，（a，b，c）為平面參數(shù)。

以入射角定權(quán)，可構(gòu)造WTLSPIA－1法與WTLSPIA－2法。對(duì)于WTLSPIA－1法，其定權(quán)方法如下：

1）選擇距離M（xi，yi，zi）最近的20個(gè)點(diǎn)，利用TLS法進(jìn)行平面擬合，求取平面參數(shù)；

2）根據(jù)式（8），點(diǎn)M（xi，yi，zi）的入射角余弦值為cosθi，令Pi＝cosθi；

3）重復(fù)步驟1）～2），直至得到所有點(diǎn)的入射角余弦值。

對(duì)于WTLSPIA－2法，可先利用TLS對(duì)觀測(cè)數(shù)據(jù)進(jìn)行擬合，之后計(jì)算各點(diǎn)的入射角余弦值，并以此作為迭代算法中各點(diǎn)權(quán)值的初始值，即Pi0＝cosθi［8－9］。

2.3 距離定權(quán)

點(diǎn)到平面的距離反映了點(diǎn)與平面的相關(guān)性，即點(diǎn)距離平面越近，其與平面的相關(guān)性越強(qiáng)，在擬合過程中所占權(quán)重應(yīng)越大；反之，點(diǎn)距離平面越遠(yuǎn)，其與平面的相關(guān)性越弱，在擬合過程中所占的權(quán)重應(yīng)越小。因此，可根據(jù)點(diǎn)到平面的距離確定該點(diǎn)的擬合權(quán)值。

本文以距離定權(quán)構(gòu)造了WTLSD 法，其先利用TLS對(duì)觀測(cè)數(shù)據(jù)進(jìn)行擬合，之后計(jì)算各點(diǎn)相對(duì)于擬合平面的距離，根據(jù)式（10）將點(diǎn)到平面距離的倒數(shù)轉(zhuǎn)換成［0，1］之間的數(shù)值，并以此作為迭代算法中各點(diǎn)權(quán)值的初始值，即Pi0＝d′i：

式中，di為點(diǎn)到平面的距離，min（di）為點(diǎn)到平面的最小距離，max（di）點(diǎn)到平面的最大距離，ε為極小值（本文取ε＝10－6mm）。

若點(diǎn)云數(shù)據(jù)在x、y、z三個(gè)方向等精度獲取，則依據(jù)以上定權(quán)方法，觀測(cè)值權(quán)陣PZ及系數(shù)陣A的行向量權(quán)陣PXY定義如下：

同時(shí)，由于系數(shù)矩陣A第3列為常數(shù)，不需要修正，則A的列向量權(quán)陣為：

3 平面參數(shù)迭代解算及精度評(píng)定

參照非線性最小二乘的牛頓－高斯迭代算法，式（2）可改寫成如下形式：

式中，λ為n×1維拉格朗日乘數(shù)，通過對(duì)式（14）求導(dǎo)、求極值可計(jì)算出λ及ξ的值。

對(duì)于WTLSPIA－2法與WTLSD 法，模型參數(shù)在迭代解算過程中均可自適應(yīng)地修正觀測(cè)向量權(quán)陣與系數(shù)矩陣權(quán)陣，同時(shí)以3倍距離標(biāo)準(zhǔn)差為閾值，剔除異常數(shù)據(jù)，其參數(shù)解算過程如下［8］：

1）在利用TLS法求得平面參數(shù)估值的基礎(chǔ)上，依據(jù)入射角定權(quán)與距離定權(quán)及式（11）、（12）設(shè)定A的列向量權(quán)陣P0及A的行向量初始權(quán)陣PXY（0）、觀測(cè)值初始權(quán)陣PZ（0）；

2）求取參數(shù)ξ＝［a，b，c］T的迭代初始值：

式中，Q0為P0的廣義逆；

6）利用式（18）計(jì)算點(diǎn)至平面的距離di、距離標(biāo)準(zhǔn)差σd：

式中，n為總觀測(cè)數(shù)，若di＞3σd，則認(rèn)為該點(diǎn)為異常點(diǎn)，刪除；反之，保留；

對(duì)于以強(qiáng)度值為權(quán)值的RWTLSI 法與WTLSPIA－1法，除無需步驟4）與6）外，解算方法與WTLSPIA－2法一致。

經(jīng)迭代運(yùn)算得到擬合平面的參數(shù)解之后，可利用式（19）計(jì)算單位權(quán)中誤差及平面擬合精度，進(jìn)行精度評(píng)定：

4 算例分析

為對(duì)本文4種定權(quán)方法與由其構(gòu)造的4種平面擬合算法的可行性及適用性進(jìn)行驗(yàn)證，利用徠卡C10三維激光掃描儀分別掃描距其30m 的標(biāo)準(zhǔn)反射板、普通木板，獲得兩組平面點(diǎn)云數(shù)據(jù)（圖1）。

圖1 實(shí)驗(yàn)樣本的點(diǎn)云數(shù)據(jù)Fig.1 Point clouds in experiments

分別利用LS法、TLS法、WTLSI法、RWTLSI法、WTLSPIA－1 法、WTLSPIA－2 法、WTLSD 法對(duì)兩組平面點(diǎn)云數(shù)據(jù)進(jìn)行擬合。本文采用奇異值分解（SVD）求解TLS參數(shù)估值［2－3］；依據(jù)式（11）、（12）設(shè)定WTLSI法、RWTLSI 法 與WTLSPIA－1法中A的列向量權(quán)陣P0及A的行向量權(quán)陣PXY、觀測(cè)值權(quán)陣PZ，并將以上權(quán)陣作為WTLSPIA－2 法與WTLSD 法對(duì)應(yīng)權(quán)陣的初始值。各方法計(jì)算得到的平面參數(shù)及精度評(píng)定指標(biāo)如表1、2所示。

表1 標(biāo)準(zhǔn)反射板平面參數(shù)及擬合精度Tab.1 Plane parameters and fitting precision of standard reflecting board

表2 普通木板平面參數(shù)及擬合精度Tab.2 Plane parameters and fitting precision of common board

表1中，由于標(biāo)準(zhǔn)反射板的反射效果較好，各算法所得平面參數(shù)基本一致，單位權(quán)中誤差、平面擬合精度及點(diǎn)距平面的最大距離均較小，擬合精度較高，說明各算法對(duì)于均質(zhì)性較好的反射板都具有較好的適用性。其中，LS法忽視了系數(shù)矩陣誤差，TLS法雖同時(shí)顧及了觀測(cè)向量與系數(shù)矩陣誤差，但未考慮各觀測(cè)值的點(diǎn)位精度，因此以上兩種算法的擬合精度評(píng)定指標(biāo)均大于各WTLS法，擬合效果相對(duì)較差；WTLSI法與RWTLSI法、WTLSPIA－1 法分別以點(diǎn)位強(qiáng)度值、點(diǎn)位入射角為權(quán)值，較合理地反映了各點(diǎn)的點(diǎn)位精度；WTLSPIA－2 法與WTLSD法在解算過程中分別以入射角、距離為定權(quán)準(zhǔn)則，自適應(yīng)地修正觀測(cè)向量權(quán)陣及系數(shù)矩陣權(quán)陣，在獲取最優(yōu)權(quán)值的同時(shí)又以3 倍距離標(biāo)準(zhǔn)差為閾值，剔除觀測(cè)數(shù)據(jù)中的異常點(diǎn)，通過迭代計(jì)算，獲得高精度參數(shù)解。

由表2 可看出，對(duì)于均質(zhì)性較差的普通木板，各方法的擬合效果差別較大。對(duì)比表1、2可知，LS法、TLS法雖擬合精度一般，但對(duì)于不同材質(zhì)的反射板，均能得到較準(zhǔn)確的平面參數(shù)解；由于普通模板的反射性差，點(diǎn)位的強(qiáng)度值已不能代表其點(diǎn)位精度，因而以強(qiáng)度值定權(quán)的WTLSI法、RWTLSI法的各擬合精度評(píng)定指標(biāo)均比較大，參數(shù)解不具備參考性。由表1、2 可知，兩種方法擬合效果一致，說明強(qiáng)度值之間不會(huì)因?yàn)榉瓷洳馁|(zhì)的不同表現(xiàn)出較大的差異性。因而，強(qiáng)度值是反射性能的體現(xiàn)，并不能較好地代表點(diǎn)位精度。WTLSPIA－1 法根據(jù)概率大小設(shè)定擬合權(quán)陣，權(quán)陣的隨機(jī)性及不確定性使其擬合效果較差。WTLSPIA－2法與WTLSD法可自適應(yīng)地修正擬合權(quán)陣，不斷地使權(quán)陣趨于合理，因此也取得了較好的擬合效果。相對(duì)LS法、TLS法、RWTLSI法及WTLSPIA－1法，以上兩種方法的單位權(quán)中誤差分別提高了99.4%、99.6%、99.4%、99.3% 與99.7%、99.8%、99.7%、99.6%，平面擬合精度分別減少了0.3、0.1、79.9、79.9mm，點(diǎn)距平面的最大距離分別減少了0.7、0.4、155、149mm。同時(shí)，對(duì)比表1、2可發(fā)現(xiàn)，WTLSD 法的單位權(quán)中誤差與平面擬合精度略優(yōu)于WTLSPIA－2法，表明該方法的擬合效果最好，平面參數(shù)解的可信度最高。

5 結(jié) 語(yǔ)

1）不論擬合何種反射材質(zhì)的平面點(diǎn)云數(shù)據(jù)，LS法與TLS由于具有良好的穩(wěn)健性，均能得到較合理的參數(shù)解。

2）WTLS法以觀測(cè)值的點(diǎn)位精度作為其擬合權(quán)重，當(dāng)選擇的定權(quán)方法可準(zhǔn)確表達(dá)各點(diǎn)的點(diǎn)位精度時(shí)，本文4 種加權(quán)平面擬合算法的擬合效果均優(yōu)于傳統(tǒng)的LS法及TLS法；若定權(quán)方法不合理，即當(dāng)點(diǎn)的擬合權(quán)重與實(shí)際權(quán)重存在偏差時(shí)，擬合效果較差，所得參數(shù)解精度較低，如WTLSI法、RWTLSI法及WTLSPIA－1法。

3）WTLSPIA－2法及WTLSD法分別依據(jù)其定權(quán)準(zhǔn)則，在參數(shù)解算的過程中，可自適應(yīng)地修正觀測(cè)向量權(quán)陣及系數(shù)矩陣權(quán)陣，使權(quán)陣合理化。因此，用以上兩種方法擬合標(biāo)準(zhǔn)反射板、普通模板兩種材質(zhì)的點(diǎn)云數(shù)據(jù)，均得到了良好的擬合效果。

4）以距離定權(quán)的WTLSD 法利用點(diǎn)與平面之間的相關(guān)性，可使參與擬合的各點(diǎn)的權(quán)更準(zhǔn)確地表達(dá)其點(diǎn)位精度，因而擬合效果最好，各項(xiàng)精度評(píng)定指標(biāo)最優(yōu)，適用性最強(qiáng)。

［1］Burkhard S，Yaron A F.On the Multivariate Total Least－Squares Approach to Empirical Coordinate Transformations［J］.Journal of Geodesy，2008，82（7）：373－383

［2］魯鐵定，周世健.總體最小二乘的迭代解法［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2010，35（11）：1 351－1 354（Lu Tieding，Zhou Shijian.An Iterative Algorithm for Total Least Squares Estimation［J］.Journal of Geodesy and Geodynamics，2010，35（11）：1 351－1 354）

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

［4］Burkhard S，Andreas W.On Weighted Total Least－Squars Adjustment for Linear Regression［J］.Journal of Geodesy，2008，82（6）：415－421

［5］Shen Y Z，Li B F，F(xiàn)elus Y A.An Iterative Solution of Weighted Total Least－Squares Adjustment［J］.Journal of Geodesy，2010，85（4）：229－238

［6］袁慶，樓立志，陳瑋嫻.基于加權(quán)整體最小二乘的點(diǎn)云數(shù)據(jù)平面擬合法［J］.測(cè)繪通報(bào)，2011（3）：1－3（Yuan Qing，Lou Lizhi，Chen Weixian.Applying Weighted Total Least－Squares to the Plane Point Cloud Fitting of Terrestrial Laser Scanning［J］.Bulletin of Surveying and Mapping，2011（3）：1－3）

［7］蒼桂華，岳建平.基于加權(quán)總體最小二乘法的點(diǎn)云平面擬合［J］.激光技術(shù)，2014（5）：307－310（Cang Guihua，Yue Jianping.Plane Fitting of Point Clouds Based on Weighted Total Least Squares Adjustment［J］.Laser Technology，2014（5）：307－310）

［8］歐江霞，李明峰，王永明.基于穩(wěn)健加權(quán)總體最小二乘的點(diǎn)云數(shù)據(jù)平面擬合［J］.大地測(cè)量與地球動(dòng)力學(xué)，2014（3）：160－163（Ou Jiangxia，Li Mingfeng，Wang Yongming.Plane Fitting of Point Clouds Based on Robust Weighted Total Least Squares［J］.Journal of Geodesy and Geodynamics，2014（3）：160－163）

［9］蒼桂華，李明峰，岳建平.以入射角定權(quán)的點(diǎn)云數(shù)據(jù)加權(quán)總體最小二乘平面擬合研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，2014（3）：95－98（Cang Guihua，Li Mingfeng，Yue Jianping.Robust Plane Fitting of Point Clouds Based on TLS［J］.Journal of Geodesy and Geodynamics，2014（6）：95－98）

［10］Bucksch A，Lindenbergh R C，Van R J.Error Budget of Terrestrial Laser Scanning：Influence of the Intensity Remission on the Scan Quality［C］.III International Scientific Congress Geo－Siberia，Novosibirsk，2007

［11］Soudarissanane S，Lindenbergh R，Teunissen M M.Incidence Angle Influence on the Quality of Terrestrial Laser Scanning Points［C］.ISPRS Workshop Laserscanning，Paris，2009