加權總體最小二乘點云平面擬合定權方法探討

李明峰 歐江霞 檀 丁 楊輝輝

1 南京工業大學地球空間信息研究中心，南京市浦珠南路30號，210009

2 武漢大學測繪學院，武漢市珞喻路129號，430079

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

1 加權總體最小二乘基本原理

點云數據平面擬合函數模型為：

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

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

隨機誤差ez與EA的統計性質如下：

式中，“?”為Kronecker積，vec（）為矩陣拉直變換，為未知方差分量，QZ、QA為ez與eA的對稱、非奇異協因數陣，且有：

式中，PZ為觀測值權陣，PA為系數陣A的權陣，P0、PXY分別為系數陣A的列向量權陣及行向量權陣。

加權總體最小二乘的參數估計準則為：

2 定權方法探討

2.1 強度值定權

對于相同材質的反射面，激光反射強度與點位精度相對應，即點位的激光反射強度越大，點位精度越高，其擬合權重越大；反之，反射強度越小，點位精度越低，其擬合權重越小。對于范圍在［－2 047，2 048］的強度值，可利用式（6）轉換成［0，1］之間的數值，即可確定為該點的權值［6－7］：

以強度值定權可構造RWTLSI法。本文顧及異常強度值對擬合精度的影響，通過式（7）計算觀測數據的強度值標準差，剔除的點，提高擬合數據精度，之后以強度值定權，構造了RWTLSI法：

式中，Ii為點位強度值，為平均強度值，σI為強度值標準差。

2.2 入射角定權

激光反射強度與入射角存在如下關系：

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

式中，（xi，yi，zi）為點i的三維坐標，（a，b，c）為平面參數。

以入射角定權，可構造WTLSPIA－1法與WTLSPIA－2法。對于WTLSPIA－1法，其定權方法如下：

1）選擇距離M（xi，yi，zi）最近的20個點，利用TLS法進行平面擬合，求取平面參數；

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

3）重復步驟1）～2），直至得到所有點的入射角余弦值。

對于WTLSPIA－2法，可先利用TLS對觀測數據進行擬合，之后計算各點的入射角余弦值，并以此作為迭代算法中各點權值的初始值，即Pi0＝cosθi［8－9］。

2.3 距離定權

點到平面的距離反映了點與平面的相關性，即點距離平面越近，其與平面的相關性越強，在擬合過程中所占權重應越大；反之，點距離平面越遠，其與平面的相關性越弱，在擬合過程中所占的權重應越小。因此，可根據點到平面的距離確定該點的擬合權值。

本文以距離定權構造了WTLSD 法，其先利用TLS對觀測數據進行擬合，之后計算各點相對于擬合平面的距離，根據式（10）將點到平面距離的倒數轉換成［0，1］之間的數值，并以此作為迭代算法中各點權值的初始值，即Pi0＝d′i：

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

若點云數據在x、y、z三個方向等精度獲取，則依據以上定權方法，觀測值權陣PZ及系數陣A的行向量權陣PXY定義如下：

同時，由于系數矩陣A第3列為常數，不需要修正，則A的列向量權陣為：

3 平面參數迭代解算及精度評定

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

式中，λ為n×1維拉格朗日乘數，通過對式（14）求導、求極值可計算出λ及ξ的值。

對于WTLSPIA－2法與WTLSD 法，模型參數在迭代解算過程中均可自適應地修正觀測向量權陣與系數矩陣權陣，同時以3倍距離標準差為閾值，剔除異常數據，其參數解算過程如下［8］：

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

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

式中，Q0為P0的廣義逆；

6）利用式（18）計算點至平面的距離di、距離標準差σd：

式中，n為總觀測數，若di＞3σd，則認為該點為異常點，刪除；反之，保留；

對于以強度值為權值的RWTLSI 法與WTLSPIA－1法，除無需步驟4）與6）外，解算方法與WTLSPIA－2法一致。

經迭代運算得到擬合平面的參數解之后，可利用式（19）計算單位權中誤差及平面擬合精度，進行精度評定：

4 算例分析

為對本文4種定權方法與由其構造的4種平面擬合算法的可行性及適用性進行驗證，利用徠卡C10三維激光掃描儀分別掃描距其30m 的標準反射板、普通木板，獲得兩組平面點云數據（圖1）。

圖1 實驗樣本的點云數據Fig.1 Point clouds in experiments

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

表1 標準反射板平面參數及擬合精度Tab.1 Plane parameters and fitting precision of standard reflecting board

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

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

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

5 結 語

1）不論擬合何種反射材質的平面點云數據，LS法與TLS由于具有良好的穩健性，均能得到較合理的參數解。

2）WTLS法以觀測值的點位精度作為其擬合權重，當選擇的定權方法可準確表達各點的點位精度時，本文4 種加權平面擬合算法的擬合效果均優于傳統的LS法及TLS法；若定權方法不合理，即當點的擬合權重與實際權重存在偏差時，擬合效果較差，所得參數解精度較低，如WTLSI法、RWTLSI法及WTLSPIA－1法。

3）WTLSPIA－2法及WTLSD法分別依據其定權準則，在參數解算的過程中，可自適應地修正觀測向量權陣及系數矩陣權陣，使權陣合理化。因此，用以上兩種方法擬合標準反射板、普通模板兩種材質的點云數據，均得到了良好的擬合效果。

4）以距離定權的WTLSD 法利用點與平面之間的相關性，可使參與擬合的各點的權更準確地表達其點位精度，因而擬合效果最好，各項精度評定指標最優，適用性最強。

［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］.武漢大學學報：信息科學版，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］官云蘭，劉紹棠，周世健.基于整體最小二乘的穩健點云數據平面擬合［J］.大地測量與地球動力學，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，Felus Y A.An Iterative Solution of Weighted Total Least－Squares Adjustment［J］.Journal of Geodesy，2010，85（4）：229－238

［6］袁慶，樓立志，陳瑋嫻.基于加權整體最小二乘的點云數據平面擬合法［J］.測繪通報，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］蒼桂華，岳建平.基于加權總體最小二乘法的點云平面擬合［J］.激光技術，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］歐江霞，李明峰，王永明.基于穩健加權總體最小二乘的點云數據平面擬合［J］.大地測量與地球動力學，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］蒼桂華，李明峰，岳建平.以入射角定權的點云數據加權總體最小二乘平面擬合研究［J］.大地測量與地球動力學，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