程效军，李 杰

（同济大学 测绘与地理信息学院，上海200092）

三维激光扫描作为一种新型的测量技术，可以快速获取对象表面的点云数据，为目标的三维重建提供了一种便捷方法.真实场景中含有大量平面特征，点云平面拟合是三维重建的基础步骤，在此基础上衍生了众多应用.蔡来良等［1］利用点云平面拟合方法进行了建筑物变形监测应用研究；陈磊［2］和叶珉吕［3］提出了基于点云平面拟合的点云滤波算法；程效军［4］利用最小二乘平面拟合方法来检测墙面平整度.因此，探讨点云平面拟合的稳健算法对三维激光扫描的应用具有重要意义.

对于含有粗差的点云平面拟合，目前常用的方法有粗差探测法、选权迭代法、稳健特征值法和随机采样一致性算法等.粗差探测法在假定只有一个观测粗差的前提下其理论上是严密的，但实际获取的点云数据往往不止一个粗差.若逐次采用粗差探测法来剔除粗差，则存在以下问题：由于粗差对每个观测值都有影响，尤其当存在多个粗差时，第一步数据探测中统计量最大的那个观测值不一定含有粗差，如果将其剔除就会造成错误的判断［5］，这就造成了数据探测法应用的局限性.选权迭代法是目前测量数据处理中剔除粗差点的稳健估计方法，可以保证所估计的参数少受模型误差尤其是粗差的影响.一般通过最小二乘估计来确定第一次平差后的残差，由于极少的粗差也可能使最小二乘方法崩溃，因此，在一定程度上削弱了选权迭代法的抗差能力，而且随着粗差含量的增加，选权迭代法得到的拟合平面将明显偏离真实值.稳健特征值法能抵抗一定含量的粗差影响，得到较为准确的平面.随机采样一致性算法的采样次数视具体数据而定，没有固定的参考值，而且恰当的阈值难以确定.

最小截断二乘法（least trimmed squares，LTS）是一种稳健估计方法，常用于统计学中处理含有局外点的线性回归问题［6-9］.LTS对局外点的影响不敏感，可以达到50%的崩溃点［10］.与LTS方法相近的还有最小中位二乘法（least median squares，LMS），Rousseeuw［8］指出，LTS回归算法比 LMS更有优势：LTS的目标函数更平滑，使其对局部影响较小；其次LTS估值近似正态，收敛速度比LMS更快，因而LTS的统计效果更好.Doornik［9］讨论了样本量与计算时间的关系，当样本量为105时，计算时间为6 800S.因此当数据量较大时，LTS估计效率不高.在处理含有粗差的点云平面拟合问题上，本文尝试在选权迭代法的基础上结合LTS，以提高选权迭代法的稳健性，同时保持较高的估计效率.

1 具有LTS稳健初值的点云平面拟合算法

常用的选权迭代法能抵抗点云中一定含量的粗差，当粗差含量达到一定值后，选权迭代法便不能得到正确的结果.为了得到统计意义上更为稳健的估计方法，同时保持较高的估计效率，考虑将LTS估值作为初始值，计算第一次平差后的残差，选择较稳健 的 绝 对 偏 差 中 位 数 MAD（the median of all absolute deviations）作为第一次计算的单位权中误差，然后再利用此残差进行选权迭代，这样既保留了选权迭代法较高的估计效率，又继承了LTS的高崩溃点，因而在一定程度上改善了选权迭代法的估计性能.

1.1 LTS估计

Rousseeuw［11］提 出 LTS 的 定 义 式 为 arg min（本文中带括号的下标表示排 序后的 顺序），将平方后的残差r2按照从小到大的顺序进行排列后得到，即r2（1）≤…≤r2（i）≤…≤r2（h）.其中，r为残差，h=int［（n+p+1）／2］，n为点的个数，p为待估参数个数.LTS的主要思想是使目标函数升序排列的前h个残差平方值之和达到最小［12］.LTS估计当h=int［（n+p+1）／2］时其崩溃点在合理情况下可达到50%，是一种非常稳健的方法.

遍历n个观测值中任意h个观测值的所有可能组合，对每种组合进行最小二乘估计，取其中最小残差平方和对应的解作为LTS估计的准确解［6］.从n个观测值中任取h个观测值一共有Chn种可能，当n和h较大时，计算量非常大，因而也导致了LTS的估计效率不高.对此，为提高计算效率，王彤等［7］提出了重复抽样算法以得到LTS的近似解.Rousseeuw等［8］利用选择迭代和数据分块的策略提出了快速LTS方法来提高LTS的效率.由于点云数据量较大，因此采用王彤等［7］提出的方法得到LTS近似解以提高计算效率.步骤如下：

（1）确定重复抽样次数SN；

（2）从n个观测值中任取p+1个观测值（p为待估参数个数，如平面方程为ax+by+cz=1，其中a，b，c为平面方程系数，则p=3）进行平面拟合，得到相应的平面方程系数；

（3）利用步骤（2）得到的平面方程计算所有点的残差平方值r2i，i=1，…，n，并排序为；

（5）重复步骤（2）～（4）直至抽样次数SN，在所得到的SN个残差平方和中取最小值，此时对应的平面方程系数即作为拟合平面的LTS估计.

1.2 权函数的确定

权函数的选择对选权迭代法的抗差性能十分关键.李德仁等［13］指出，关于权函数的选择应满足三点要求：①通过迭代，含粗差观测值的权应逐步趋近于零；②不含粗差观测值的权，在迭代终止时应等于该组观测值的权；③权函数的选择应保证迭代过程能以较快的速度收敛.常见的权函数有Huber权函数、Hampel权函数、丹麦权函数、IGG（institute of geodesy and geophysics）权函数等.其中 Huber权函数和丹麦权函数没有淘汰段，因而抗差能力减弱.Hampel权函数将阈值区间分为四段，计算较为复杂.王任享等［14］指出残差并不能使改正数在同一概率水平，标准化残差属于同方差母体的子样，用标准化残差代替残差作为权函数的参数在统计判断上更为严密.标准化残差的计算式为

式中：ωi为标准化残差；ri为残差；gii为标准化因子，其值为残差协因数阵Qvv的主对角线元素；σ0为单位权中误差.因此采用改进的IGG权函数.

式中：P（ω）为权函数；ω为标准化残差；k为常数.

第一次平差计算时单位权中误差可以采用Rousseeuw［15］推荐的尺度估计S：

式中：1.482 6为调制因子；n为观测值个数；p为待估参数个数；r为残差；median为r2i（i=1，2，…，n）的中值.也可以采用较稳健的绝对偏差中位数MAD作为初始的单位权中误差M［16］.

其中，1.483为改正因子，使MAD在正态分布情况下无偏，r为残差，median表示取中值.本文采用MAD作为初始单位权中误差.综上所述，具有LTS稳健初值的点云平面拟合算法的计算步骤如下：

（1）计算点云平面方程系数的LTS初值，迭代次数i=1，得到残差r（1），并计算初始单位权中误差M；

（2）计算标准化残差ω，根据IGG权函数确定权阵P；

（3）迭代次数i=i+1，计算平面方程系数的加权最小二乘估计p（i）p，并计算残差r（i）；

（6）输出平面方程系数的最佳估值pp.

2 算法验证

为了验证本文算法的可靠性和稳健性，利用地面三维激光扫描仪扫描一平面，截取其中的一部分作为实验数据，所有算法采用Matlab编程实现.为检验拟合平面的准确性，在点云处理软件Geomagic Studio中进行去噪，得到净化的点云数据（图1），再根据最小二乘方法计算得到标准平面方程为：

对粗差的模拟，利用Matlab向净化的点云数据中随机添加粗差，粗差含量从5%逐次添加到40%（如图2所示为随机添加5%粗差后的点云）.首先采用选权迭代法和特征值迭代法对不同粗差含量的点云数据进行拟合，结果见表1和表2，表中，a，b，c为平面方程系数，σ为中误差.

图1 不含粗差的点云平面图Fig.1 Point clouds without outliers

图2 含5%粗差的点云平面图Fig.2 Point clouds with 5%outliers

表1 选权迭代法的拟合结果Tab.1 Fitting result of selecting weight iteration

根据图表可以得出：当粗差含量在15%以内时，选权迭代法和特征值迭代法均稳健，当粗差含量达到20%时，选权迭代法和特征值迭代法的计算结果与标准系数的偏差分别为21.99%和7.04%，迭代崩溃.为定出这两种方法的迭代崩溃点，对粗差含量为16%～19%的点云再次进行测试.试验结果见表3和表4.若以偏差2%为限，则选权迭代法的崩溃点可定为18%，特征值迭代法的崩溃点可定为17%.将表1和表2中的数据绘制成图，如图3所示.

表2 特征值迭代法的拟合结果Tab.2 Fitting result of eigenvalue iteration

图3 与标准系数的偏差及迭代次数的变化Fig.3 Deviation from standard plane and iteration times variation

表3 选权迭代法的拟合结果（16%～19%粗差含量）Tab.3 Fitting result of selecting weight iteration（With 16%～19%outliers）

表4 特征值迭代法的拟合结果（16%～19%粗差含量）Tab.4 Fitting result of eigenvalue iteration（With 16%～19%outliers）

试验结果表明，随着粗差含量的增加，迭代计算并不能取得满意的结果，当粗差含量达到崩溃点后，由于初始迭代模型偏差过大，导致迭代失败；两种方法解算出的系数与标准系数的偏差总体趋势是增加的，但由于粗差的随机性，在某些点的偏差有所减小；选权迭代法的迭代次数随着粗差含量的增加先增加后减小，而特征值迭代法的迭代次数相对较稳定.因此，当粗差含量在崩溃点以内时，选权迭代法和特征值迭代法都可以得到准确的结果.在实际应用中，由于点云的粗差含量事先无法预知，为了验证本文算法的有效性，对算法的测试从粗差含量为5%开始.

在计算LTS初值时，首先需要确定重复抽样次数.对平面拟合来说，每次应从所有观测值中随机抽取4个观测值进行计算.设粗差含量为ν，则抽样一次抽到粗差的概率为ν，随机抽取u次，则抽到不含粗差的样本的概率为P=1-νu，当样本的粗差含量为40%时，取u=40，计算可得P≈1，因此在测试数据中确定重复抽样次数为40.LTS初值结果见表5.

表5 LTS初值Tab.5 LTS initial value

由于LTS初值与标准平面系数仍有偏差，因此在LTS初值基础上再进行选权迭代，对模型进行精炼.最终获得平面方程系数见表6.

根据表6可以得出：对于粗差含量较低的点云，具有LTS初值的选权迭代法可以得到准确的平面方程系数，其精度与选权迭代法和特征值迭代法相当.当粗差含量为20%～40%时，本文算法依然稳健，且点云平面拟合系数与标准系数的偏差均小于2%，满足精度要求.在整个计算过程中，最大迭代次数为34次.因此，具有LTS初值的选权迭代法较上述两种迭代法具有更高的崩溃点，更为稳健，且保证了计算效率.

表6 具有LTS初值的选权迭代法的拟合结果Tab.6 Fitting result of selecting weight iteration with LTS initial value

为探讨20%和35%处出现的偏差波动现象是由粗差的随机性引起，还是由系统误差造成，选择另一平面按照具有LTS初值的选权迭代法进行试验.标准平面方程为：0.023 910 4x+0.133 495 1y-0.003 445z=1，试验结果见表7.

表7 具有LTS初值的选权迭代法的拟合结果Tab.7 Fitting result of selecting weight iteration with LTS initial value

根据表7可知，点云平面拟合系数与标准系数的偏差均小于2%，满足精度要求.随着粗差含量的增加，偏差值并未呈现相应的增大，而是呈现不规律的波动，故推断由于粗差的随机性而导致偏差变化的不规律.

3 结语

对于粗差含量不高的点云平面拟合问题，采用选权迭代法、特征值迭代法和具有LTS初值的选权迭代法都可以获得理想的结果.当粗差含量增加到一定值后，由于初始迭代模型偏差过大，导致选权迭代法和特征值迭代法失效.具有LTS初值的选权迭代法可以获得稳健的平面初值，在此基础上再进行选权迭代，这样既采纳了LTS的高崩溃点，又利用了选权迭代法的高估计效率，从而获得较为精确的平面方程系数.在进行LTS初值计算时采用了随机抽样的方法求得近似解，以提高计算效率，这对于数据量庞大的点云来说十分重要.实验表明，当粗差含量达到很高的比例（如40%）时，具有LTS初值的选权迭代法依然稳健.由于本文算法采用的平面方程为通式ax+by+cz=1，因此该算法可用于水平或铅垂等特殊点云平面拟合.如果将平面模型改为曲线曲面等模型，则具有LTS稳健初值的选权迭代法还可以应用到其他工程数据拟合方面.

［1］ 蔡来良，吴侃，张舒.点云平面拟合在三维激光扫描仪变形监测中的应用［J］.测绘科学，2010，35（5）：231.CAI Lailiang，WU Kan，ZHANG Shu.Application of point cloud plane fitting to deformation monitoring using 3D laser scanner［J］.Science of Surveying and Mapping，2010，35（5）：231.

［2］ 陈磊，赵书河.一种改进的基于平面拟合的机载LiDAR点云滤波方法［J］.遥感技术与应用，2011，26（1）：117.CHEN Lei，ZHAO Shuhe.An improved plane fitting based filtering algorithm for airborne LiDAR data［J］.Remote Sensing Technology and Application，2011，26（1）：117.

［3］ 叶珉吕，花向红，陈西江，等.基于正交整体最小二乘平面拟合的点云数据去噪方法研究［J］.测绘通报，2013（11）：37.YE Minlv，HUA Xianghong，CHEN Xijiangetal.Research on method for the denoising of point cloud based on orthogonal total least squares fitting［J］.Bulletin of Surveying and Mapping，2013（11）：37.

［4］ 程效军，唐剑波.基于最小二乘拟合的墙面平整度检测方法［J］.测绘信息与工程，2007，32（4）：19.CHENG Xiaojun，TANG Jianbo.Method for estimating metope smoothing grade based on least squares fitting［J］.Journal of Geomatics，2007，32（4）：19.

［5］ 张勤，张菊清，岳东杰.近代测量数据处理与应用［M］.北京：测绘出版社，2011.ZHANG Qin，ZHANG Juqing，YUE Dongjie.Advanced theory and application of surveying data［M］.Beijing：Surveying and Mapping Press，2011.

［6］ Giloni A，Padberg M.Least trimmed squares regression，least median squares regression，and mathematical programming［J］.Mathematical and Computer Modeling，2002，35（9）：1043.

［7］ 王彤，何大卫.LTS回归与M估计稳健性的比较［J］.中国卫生统计，1999，16（2）：79.WANG Tong，HE Dawei.Comparison of M-estimator and LTS estimator in regression ［J］.Chinese Journal of Health Statistics，1999，16（2）：79.

［8］ Rousseeuw P J，Van Driessen K.Computing LTS regression for large data sets［J］.Data Mining and Knowledge Discovery，2006，12（1）：29.

［9］ Doornik J A.Robust estimation using least trimmed squares［EB／OL］.（2011-03-14）［2014-08-08］.http：／／econ.au.dk／fileadmin／site＿files／filer＿oekonomi／subsites／creates／Seminar＿Papers／2011／ELTS.pdf.

［10］ Mount D M，Netanyahu N S，Piatko C D，etal.On the least trimmed squares estimator［J］.Algorithmica，2014，69（1）：148.

［11］ Rousseeuw P J.Least median of squares regression［J］.Journal of the American Statistical Association，1984，79（388）：871.

［12］ 杨飚，张曾科，孙政顺.基于LTS稳健初值的选权迭代法［J］.科学技术与工程，2005，5（22）：1720.YANG Biao，ZHANG Zengke，SUN Zhengshun.Selecting weight iteration method with initial value by LTS［J］.Science Technology and Engineering，2005，5（22）：1720.

［13］ 李德仁，袁修孝.误差处理与可靠性理论［M］.武汉：武汉大学出版社，2002.LI Deren，YUAN Xiuxiao.Error processing and reliability theory［M］.Wuhan：Wuhan University Press，2002.

［14］ 王任享.选权迭代定位粗差时权函数参数之功能［J］.武汉测绘科技大学学报，1988，13（4）：42.WANG Renxiang.Effects of parameters of weight function for blunder detection by the recursive weighted least squares method［J］.Journal of Wuhan Technical University of Surveying and Mapping，1988，13（4）：42.

［15］ Rousseeuw P J，Leroy A M.Robust regression and outlier detection［M］.New Jersey：John Wiley &Sons，2005.

［16］ Rousseeuw P J，Hubert M.Robust statistics for outlier detection［J］.Wiley Interdisciplinary Reviews：Data Mining and Knowledge Discovery，2011，1（1）：73.