李明峰 欧江霞 檀 丁 杨辉辉

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