苍桂华 李明峰 岳建平

1)南京工业大学测绘学院，南京 210009

2)河海大学地球科学与工程学院，南京 210098

以入射角定权的点云数据加权总体最小二乘平面拟合研究*

苍桂华1)李明峰1)岳建平2)

1)南京工业大学测绘学院，南京 210009

2)河海大学地球科学与工程学院，南京 210098

顾及入射角对点云数据点位精度的影响，提出入射角定权的点云数据加权总体最小二乘平面拟合方法。该方法以加权总体最小二乘方法为基础，利用入射角信息确定点云数据中各点的平面拟合权重，并根据系数矩阵结构特点定义其权阵。通过3种不同材质平面样本数据的计算分析表明，与最小二乘法、总体最小二乘法相比，该方法可得到更高的平面拟合精度。

点云数据;平面拟合;加权总体最小二乘;入射角定权;系数矩阵

地面三维激光现实扫描场景中含有大量平面特征，可用于前期数据预处理或后期建模的数据简化，还可用于多视点云数据配准，实现建模。平面特征的提取是利用点云数据进行3D建模的一项重要工作，提取精度直接关系到后续应用成果的精度。以往提取平面特征常用的方法为最小二乘平面拟合法。此方法根据点的三维坐标(xi，yi，zi)建立平面方程z=ax+by+c，假设误差e只存在于观测向量Z 中，建立高斯-马尔科夫(Gauss-Markov，G-M)模型，采用最小二乘(least squares，LS)方法拟合平面，获取平面拟合参数。然而，由于仪器设备、外界环境、地物特性等因素影响，使得实际点云数据的xi、yi、zi三个观测值中均存在误差，从而使包含变量xi、yi的系数矩阵A也含有误差。因此，利用最小二乘法进行平面数据拟合的结果不是最优，而是有偏的［1］。针对这种观测向量和系数矩阵均包含误差的模型，即所谓的变量中的误差(error-in-variables，EIV)模型，Golub等［2］最早提出了总体最小二乘(total least squares，TLS)估计方法。但总体最小二乘估计仅在设计矩阵和残差元素均服从独立等精度分布时才是最优估计［3］，而实际点云数据各点的精度是不等的，因此简单的总体最小二乘方法并非最优估计［4］。为了解决矩阵和残差的不等精度估计问题，Markovsky等［5］提出了加权总体最小二乘(weighted total least squares，WTLS)方法。Schaffrin等［6］则进一步扩展了WTLS方法。国内学者也从理论、应用等方面对 WTLS方法进行了深入研究［7－10］。本文根据入射角对点云数据点位精度的影响，提出利用入射角定权的加权总体最小二乘点云数据平面拟合方法，简称角度加权总体最小二乘(angle weighted total least squares，AWTLS)方法。

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

点云数据经过粗差剔除预处理后得到一组平面点云的三维坐标(xi，yi，zi)(i=1，2，…，n)。设空间平面方程式为:

式中，a、b、c为待求的平面拟合参数。将式(1)写成矩阵形式为:

若同时考虑观测向量Y和系数矩阵A的误差，则建立EIV函数模型［6］:

式中，eY和EA分别表示观测向量Y和系数矩阵A的随机误差矩阵;vec是指将矩阵按列拉直所得到的列向量，排列顺序为从左到右;eA=vec(EA)∈R3n×1表示矩阵EA按列向量化后得到的向量;σ20为未知的单位权方差;QY=P－1Y，QA=P－1A，分别表示观测向量和系数阵列向量化向量的协因素阵。

根据eA及其相关权阵的设置，采用不同的估计方法进行平面拟合参数的求解。

1)eA≡0。此时只考虑观测向量Y的误差，则EIV模型变为经典测量平差的G-M模型:

应采用最小二乘法进行参数估计。

2)eA≠0，但 QY=In，QA=Inm。此时假设点云数据坐标独立且等精度，即观测向量权阵PY和系数矩阵权阵PA均为单位阵，采用总体最小二乘法进行参数估算。

3)eA≠0，QY、QA均不为单位阵。此时，假设点云数据坐标为不等精度，采用加权总体最小二乘法进行参数估计。根据系数矩阵的结构特点设计权阵P0、Px、PY。P0是3×3矩阵，代表系数矩阵 A 的列向量权阵;PX是n×n矩阵，代表系数矩阵A的行向量权阵;PY是n×n矩阵，代表矩阵Y的权阵。P0、PX、PY相应的协因素矩阵为 Q0、QX、QY，即

式中，⊗表示“kronecker积”。

3 种估算方法的估计准则和单位权方差估计如表1所示。

表1 3种方法的估算准则和单位权方差估计Tab.1 Estimation criterions and unit weight estimation variances with three methods

2 AWTLS方法权阵的确定及其步骤

2.1 权阵的确定

点云数据中各点的精度不等，进行平面拟合时应根据各点精度确定其拟合权重。点位精度越高，参与拟合的权重也越大，从而保证结果的精度。相关文献研究了入射角对点云数据精度的影响，指出点云数据入射角越小则点位精度越高［11－12］。而入射角越小，其余弦值越大，因此可将各点入射角余弦值作为其拟合权重。平面点云数据的各点入射角θ定义为入射光线与平面法向量之间的夹角［13］。设拟合的平面方程如式(1)所示，则平面上i点的入射角余弦值或其拟合权值为:

式中，θi为 i点的入射角，(xi，yi，zi)为点 i的三维坐标，(a，b，－1)为拟合平面的法向量。

设点云在x、y、z三个方向等精度获取，对于平面的系数阵列向量和观测向量，有σx=σy=σz。结合系数矩阵A的特点，设置相应权阵:

式中，P0的第三个对角元素为0，表示系数矩阵A中第三列不需要改正;其余对角线元素为1，表示系数矩阵A中的第一、二列中的元素是等精度获取的;PX、PY与入射角有关。

2.2 方法步骤

由式(9)可以看出，权阵PX、PY与平面拟合参数相关，因此需要迭代求解。AWTLS方法的具体步骤如下:

1)将最小二乘平面拟合参数作为起始解，根据式(8)计算各点的入射角余弦值，组成相关矩阵P0、PX、PY。

5)重新计算每点的入射角余弦值，设置相关权阵 P0、PX、PY，计算(i+1):

3 算例分析

为验证AWTLS方法的可行性，采用3种不同材质的平面点云样本数据进行实验。实验的平面样本有标准反射板(反射率为35%)、普通木板及一般的水泥建材模板。3种材质样本中，标准反射板均质性最好，板面最平整;普通木板均质性和平面度次之;建材水泥模板均质性和平面度最差，但其平面偏差目测在1 mm以内。利用徕卡HDS ScanStation2分别对样本进行扫描，获取平面点云数据(图1)。

根据点云数据特点确定平面形式［14］。分别利用最小二乘法、总体最小二乘法和AWTLS法对各样本点云数据进行平面拟合，获得平面拟合参数、、以及单位权中误差(0)。假设确定的平面形式如式(1)所示，由式(15)计算点云数据各点到拟合面的距离(di)以及平面拟合精度(P)，并得到点到拟合面的最大距离(max(di)):

式中，n为拟合平面上点的个数。

将单位权中误差、点到拟合面的最大距离、平面拟合精度作为估算方法优劣的评判指标。用3种方法对3种材质的样本数据进行平面拟合，结果如表2所示。

从表2可以看出，不论何种材质样本，LS方法的单位权中误差都远远大于TLS方法和AWTLS方法，而TLS和AWTLS方法的单位权中误差为同一数量级，这是因为LS方法与TLS方法以及AWTLS方法的单位权中误差计算方法不同所致。另外，对于3种材质的样本数据，利用AWTLS方法获得的3个精度判定指标均要比LS方法和TLS方法小得多。以平面度和均质性最差的水泥模板为例，利用AWTLS方法得到的单位权中误差比TLS方法提高了65%，平面拟合精度比LS方法和TLS方法分别提高了63%和19%，点到拟合面的最大距离也由LS方法和TLS方法的9.1 mm和7.0 mm 降为5.3 mm。

图1 实验样本的平面点云数据Fig.1 Plane point clouds used in the experiments

表2 三种方法下平面点云数据拟合的相关计算结果Tab.2 Fitting results of plane samples

4 结论

1)AWTLS方法是在加权总体最小二乘方法基础上，利用入射角对点云数据精度的影响确定各点的拟合权重。从实验结果看，该定权方法符合点云数据平面拟合的实际，结果也更加合理。

2)由于LS方法、TLS方法以及AWTLS方法的单位权中误差计算方法不同，因此不能简单地以单位权中误差的大小来判定某个估算方法的优劣，应结合其他精度评判指标进行评定。

3)AWTLS方法系数矩阵的行向量权阵PX和观测向量权矩阵PY与拟合参数有关，因此在计算过程中需要迭代运算。

利用WTLS方法进行点云数据平面拟合的关键是如何确定各点的拟合权值。本文只进行了利用入射角定权的WTLS方法在点云数据平面拟合中的适用性研究，后续将进一步研究其他权值确定方法，并与本文的AWTLS方法进行比较，研究各方法的优劣。

1 邱卫宁，等.测量数据处理理论与方法［M］.武汉:武汉大学出版社，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 周拥军，邓才华.加权和不加权TLS方法及其在不等精度坐标变换中的应用［J］.武汉大学学报:信息科学版，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 王乐洋，许才军，鲁铁定.病态加权总体最小二乘平差的岭估计解法［J］.武汉大学学报:信息科学版，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 赵辉.基于加权总体最小二乘的GPS数据高程拟合［J］.大地测量与地球动力学，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 周拥军，朱建军，邓才华.附参数的条件平差与按行独立的加权总体最小二乘法估计的一致性研究［J］.测绘学报，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 陈玮娴，等.加权总体最小二乘在三维激光标靶拟合中的应用［J］.大地测量学与地球动力学，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 王解先，季凯敏.工业测量拟合［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

国家自然科学基金项目(41174002);江苏省测绘科研项目(JSCHKY201303)。

苍桂华，女，1971年生，讲师，主要研究方向为摄影测量与遥感、大地测量与测绘工程数据处理。E－mail:ghuac@163.com。