苍桂华，岳建平

引 言

点云数据平面拟合是地面3维激光扫描技术应用中常见问题，如建筑物点云数据的特征线提取等［1］。最小二乘(least squares，LS)法是常用的平面拟合方法。设利用点云数据建立的平面方程的矩阵形式为Y=AX，最小二乘法是假设误差e只存在于观测向量Y中，建立经典的高斯-马尔科夫(Gaussian-Markov，G-M)模型对误差方程式进行求解，获得平面参量估值。然而由于模型误差、人为误差、仪器误差等因素影响点云数据在观测数据x，y和z这3个方向上均存在误差，使得包含观测数据的系数矩阵A也含有误差。因此，利用最小二乘法进行平面数据拟合的结果并不是最优，而是有偏的［2］。针对这种观测向量和系数矩阵均包含误差的模型，即所谓的变量中的误差(error-in-variables，EIV)模型，GOLUB等人提出了总体最小二乘(total least squares，TLS)估计方法［3］，但总体最小二乘估计仅在设计矩阵和残差元素均服从独立等精度分布时才是最优估计［4］。实际点云数据中各点坐标精度是不等的，因此简单的总体最小二乘方法并非最优估计［5］。为了解决矩阵和残差的不等精度估计问题，MARKOVSKY等人提出了加权总体最小二乘(weighted totalleastsquares，WTLS)方 法［6］。SCHAFFRIN等人则进一步扩展了WTLS方法，详细介绍了相关权阵的设计方法以及算法的步骤［7］。本文中在SCHAFFRIN等人提出的WTLS方法基础上，根据点云数据中各点反射强度值与其点位精度关系，确定各点的平面拟合权值，得到强度加权总体最小二乘(intensity weighted total least squares，IWTLS)的平面拟合方法，并通过均质性不同的3种样本点云数据，对该方法的适用性进行研究。

1 IWTLS方法及其函数模型

1.1 EIV 模型

设点云数据所建立3维空间平面方程式形式为:

式中，a，b和c为待求的平面拟合参量。

将(1)式写成矩阵形式为:

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

根据系数矩阵特点引入权阵P0，PX和PY。P0是3×3矩阵，代表系数矩阵A的列向量权阵;PX是n×n矩阵，代表系数矩阵A的行向量权阵;PY是n×n矩阵，代表向量 Y的权阵。P0，PX，PY相对应的协因素矩阵为 Q0，QX，QY，即:

由 Q0，QX可以得到 QA，PA:

式中，⊗表示“kronecker积”。

IWTLS估算准则为:

1.2 IWTLS权阵的设置

激光反射强度值与入射角关系为Ii=I0cosαi(i=1，2，…，n)，Ii代表入射角为 αi时点的强度值，I0为垂直入射时(αi=0)点的反射强度值。入射角越小，点位精度越高，点的反射强度值越大［8-10］。因此强度值越大，参与拟合的权重应越大。本次实验数据为.PTS格式，以12bit记录强度值，其强度值范围为［－2047，2048］。设记录的原始强度值为Ii'，按照下式将其值变为［0，1］之间，构成各点的拟合权值:

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

式中，P0的第3个对角元素为0，表示系数矩阵A中的第3列不需要改正，其余对角线元素为1，表示系数矩阵A中的第1列和第2列的数据列中的元素是等精度获取的;PX，PY与强度值有关。

1.3 平面参量求解及精度评定［7，9-10］

目前解决WTLS问题主要采用基于拉格朗日乘数法的迭代解算方法［7］，计算步骤如下。

(1)根据(8)式计算出各点强度值，并根据(9)式设置相关矩阵 P0，PX，PY。

2 算例及分析

为了验证IWTLS方法的适用性，采用了均质性不同的3种平面样本进行实验。实验中的平面样本分别为标准反射板(反射率为90%)、普通木板及一般的水泥建材模板。3种平面样本中，标准反射板均质性最好;普通木板次之;建材水泥模板均质性最差。利用徕卡ScanStation2 C10分别对样本进行扫描，获取点云数据(如图1所示)。

根据点云数据特点确定平面方程式形式［11］。分别利用LS法、TLS法和IWTLS法对各个样本点云数据进行平面拟合，获得平面拟合参量，，以及单位权中误差。设拟合平面上点的个数为n，计算出各点i(i=1，2，…，n)到拟合面的距离di，获得点到拟合面的最大距离di，max，根据下式计算出平面拟合精度:

将单位权中误差、点到拟合面的最大距离以及平面拟合精度作为估算方法优劣的评判指标。标准反射板、普通木板、水泥模板3种样本数据采用不同平面拟合方法得到的相关结果如表1所示。

Fig.1 Plane fitting of point clouds in experiments

Table 1 Results of plane fitting samples

从表1可以看出:对于均质性较好的标准反射率板和普通木板，利用IWTLS方法获得的3个精度判定指标(σ^0，σ^p，di，max)值均要比 LS 方法和 TLS 方法的相应结果小得多。以普通木头为例，利用IWTLS方法得到的单位权中误差比LS方法和TLS方法分别提高了99%和25%，平面拟合精度比LS方法和TLS方法分别提高了63%和40%，点到拟合面的最大距离也由LS方法和TLS方法的8.3mm和5.8mm降为5.1mm。然而对于一般的建材水泥模板，由于均质性较差，IWTLS 方法计算出的σ^0，σ^p和di，max3 个精度评判指标分别为 0.000323m，3.2mm和6.5mm，虽好于LS方法的相应结果0.037647m，4.2mm和 8.5mm，却比 TLS方法的相关结果0.000304m，1.8mm 和5.8mm 差。这是由于此时各点强度值的差异更多由于材质不同造成，其强度值已经不能代表其点位精度，因此，利用IWTLS方法获得的相关结果差于总体最小二乘方法。

3 小结

点云数据平面拟合中WTLS法虽然从理论上较LS法和TLS法合理，但在拟合时应注意各点拟合权值的设置。确定的各点拟合权值应与其点位精度一致，即点位精度高，拟合权值应越大。如果拟合权重与其实际点位精度情况不一致，会直接影响WTLS的效果。本文中的强度加权总体最小二乘法对于均质性较好的点云平面效果明显，拟合精度较高，而对于均质性较差的点云平面效果不佳，此时应采用TLS方法进行平面拟合。

［1］ YU H X，WU K，AO J F，et al.Extraction of building’s feature lines based on 3-D laser scanning technology［J］.Laser Technology，2012，36(4):553-556(in Chinese).

［2］ QIU W N，TAO B Z，YAO Y B，et al.The theory and method of surveying data processing［M］.Wuhan:Wuhan University Press，2008:161-175(in Chinese).

［3］ GOLUB G H，van LOAN C F.An analysis of the total least squares problem［J］.SIAM Journal on Numerical Analysis，1980，17(6):883-893.

［4］ van HUFFEL S，VANDEWALLE J.The total least squares problem:computational aspects and analysis［M］.Philadelphia，USA:Society for Industrial and Applied Mathematics，1991:263-283.

［5］ ZHOU Y J，DENG C H.Weighted and unweighted total least square methods and applications to heteroscedastic 3-D coordinate transformation［J］.Geomatics and Information Science of Wuhan University，2012，37(8):976-979(in Chinese).

［6］ MARKOVSKY I，RASTELLO M，PREMOLI A，et al.The element-wise weighted total least-squares problem［J］.Computational Statistics and Data Analysis，2006，50(1):181-209.

［7］ SCHAFFRIN B，WIESER A.On weighted total least-square adjustment for linear regression ［J］.Journal of Geodesy，2008，82(7):415-421.

［8］ ZHANG Y.Research on point cloud processing of terrestrial laser scanning［D］.Wuhan:Wuhan University，2008:43-47(in Chinese).

［9］ CHEN W X，CHEN Y，YUAN Q，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(in Chinese).

［10］ YUAN Q，LOU L Z，CHEN W X.Appling weight total leastsquares to the plane point cloud fitting of terrestrial laser scanning［J］.Bulletin of Surveying and Mapping，2011(3):1-3(in Chinese).

［11］ WANG J X，JI K M.Industrial surveying fitting［M］.Beijing:Surveying and Mapping Press，2007:43-45(in Chinese).