60工程勘察GeotechnicalInvestigation＆Surveying2017年第12期一种圆柱面拟合方法袁建刚1，潘轶*2(1.江苏城乡建设职业学院，江苏常州213147;2.河海大学地球科学与工程学院，南京211100)摘要:针对建(构)筑物点云建模过程中可能碰到的圆柱面拟合问题，本文探讨了一种基于几何特性的圆柱面拟合方法。该方法利用主成分分析法和线性最小二乘法确定拟合初值，再利用非线性最小二乘法循环迭代求解圆柱模型参数，实现圆柱面拟合。它的特点在于误差方程式的建立，本文算法在传统算法的基础上对误差方程稍作改进。实验结果表明，相对于传统方法，采用本文算法计算所得结果的残差与标准差有所降低，精度提高。关键词:三维激光扫描;点云建模;柱面拟合;主成分分析;线性最小二乘;非线性最小二乘中图分类号:TB113文献标识码:AAcylindricalsurfacefittingmethodYuanJiangang1，PanYi*2(1.JiangsuUrbanandＲuralConstructionCollege，Changzhou213147，China;2.CollegeofEarthSciencesandEngineering，HohaiUniversity，Nanjing211100，China)Abstract:Inordertoaddressthecylindersurfacefittingproblemduringthepointcloudmodeling，thispaperproposesanovelsurfacefittingmethodbasedonthegeometriccharacteristics．Thismethodfirstusestheprincipalcomponentanalysisandlinearleastsquaresmethodtodeterminetheinitialvalue，andthen，aniterationsolutionofthecylindricalmodelparameterscanbeachievedbymeansofnonlinearleastsquaresmethod．Itischaracterizedbyimprovingtheerrorequationofthetraditionalalgorithmslightly．Theexperimentalresultsshowthat，comparedtothetraditionalmethod，theresidualsandthestandarddeviationofthefittingarereduced，andthefittingaccuracyisimproved．Keywords:3Dlaserscanning;pointcloudmodeling;cylindricalfitting;principalcomponentanalysis;linearleastsquares;nonlinearleastsquares收稿日期:2016-10-13;修订日期:2016-12-06作者简介:袁建刚(1974－)，男(汉族)，江苏靖江人，硕士，副教授．*通讯作者:潘轶(1994－)，硕士，E-mail:834920245@qq.com．0引言在大型工程的建设过程中，需要对建(构)筑物的形体进行测量，以确保工程施工质量。实际应用中，多数建筑物由规则几何体构成，如长方体、圆柱体等［1］。采用基于三维激光扫描技术的形体测量系统对建筑物进行扫描，对获取点云进行模型重建的过程中，可能会碰到圆柱面拟合问题［2］。圆柱面拟合方法主要包含遗传算法、Gause-Newton优化算法、特征值算法、距离函数参数化算法等［3～5］。这些算法大多原理复杂、计算繁琐，而根据圆柱面几何原理进行拟合，简便易行，易于理解，非常适用于平时的点云数据处理。因此，本文基于圆柱面的几何特性，对误差方程稍作改进，研究了一种圆柱面拟合方法，拟合出圆柱面的空间位置和形态参数，并对使用该方法所得结果的残差、标准差与传统方法进行了对比分析。1圆柱方程由圆柱面的几何特性可得，圆柱面上的点到其轴线的距离恒等于半径r0［6］，如图1所示。其中，P为圆柱面上任意一点，P0(x0，y0，z0)为圆柱轴线上一点，(a，b，c)为圆柱轴线向量，r0为圆柱底圆半径。圆柱面上任意一点到其轴线的距离为半径r0，即:(x－x0)2+(y－y0)2+(z－z0)2－［a(x－x0)+b(y－y0)+c(z－z0)］2=r20(1)

2017年第12期工程勘察GeotechnicalInvestigation＆Surveying61图1圆柱几何特性Fig.1Thecylindricalgeometrycharacteristics式(1)即为圆柱面方程。根据圆柱面上点的坐标(x，y，z)，求出圆柱轴线上一点(x0，y0，z0)、圆柱轴线向量(a，b，c)、圆柱底圆半径r0这七个参数，就可以唯一确定一个圆柱。2本文柱面拟合方法原理圆柱拟合步骤主要包括两步:一是确定柱面模型参数初始值;二是建立误差方程式求解参数值。本文算法结合主成分分析法与线性最小二乘法，确定圆柱轴线向量(a，b，c)、圆柱轴线上一点(x1，y1，z1)、圆柱底圆半径r这七个柱面模型参数初始值，再建立改进误差方程式，求解参数。采用本文算法实现柱面拟合的流程如图2所示。图2圆柱面拟合步骤图Fig.2Stepsofthecylindricalsurfacefittingmethod本文柱面拟合方法详细步骤如下:(1)确定柱面模型参数初始值。圆柱拟合初值的估计原理［7］如下:搜寻圆柱面上任意一点的若干邻近点，将这些点拟合成平面，得到的平面法向量单位化即该点的单位法向量。对圆柱面上每个点处理，得到每个点的单位法向量。将每个点的单位法向量看成点，将这些点拟合成平面，得到平面法向量，即圆柱轴线向量初始值a0、b0、c0，该步骤使用了主成分分析法［8］进行求解。先对圆柱进行坐标转换，使圆柱轴线向量(a0，b0，c0)变换为平行Z轴的向量即为(0，0，(a0)2+(b0)2+(c0)槡2)。设旋转矩阵为Ｒ3×3，要使:Ｒ3×3a0b0c0=00(a0)2+(b0)2+(c0)槡2即Ｒ3×3=b0(a0)2+(b0)槡2－a0(a0)2+(b0)槡20a0*c0(a0)2+(b0)槡2b0*c0(a0)2+(b0)槡2－(a0)2+(b0)槡2a0b0c0将旋转矩阵与圆柱面上的点的坐标(x，y，z)相乘，得到的新坐标(x'，y'，z')，将新坐标中的x'、y'取出得(x'，y')。(x'，y')为新坐标系下圆柱底面圆的坐标，根据(x'，y')可拟合一个圆。圆的方程为(x－x01')2+(y－y01')2=(r0)2。若知道圆心x01'、y01'、半径r0这三个参数，就可以唯一确定并拟合一个圆面。首先，圆坐标(x'，y')准备完毕;其次，找误差方程，将圆的方程展开得到x2+(x01')2－2x(x01')+y2+(y01')2－2y(y01')－(r0)2=0，即2x(x01')+2y(y01')+(r0)2－(x01')2－(y01')2－(x2+y2)=0。可以设x01'、y01'、(r0)2－(x01')2－(y01')2为待求参数。于是，误差方程可列为V=2x(x01')+2y(y01')+(r0)2－(x01')2－(y01')2－(x2+y2)，写成矩阵形式即:Vn×1=Bn×3X3×1－Ln×1(1)其中，B=2x12y112x22y212xn2yn1(2)

62工程勘察GeotechnicalInvestigation＆Surveying2017年第12期X=x01'y01'(r0)2－(x01')2－(y01')2(3)L=x21+y21x22+y22x2n+y2n(4)假设获取的圆的点等精度，即等权。根据线性最小二乘法原理，即VTV=min，求得未知参数X:X=(BTB)－1BTL(5)求出新坐标系下圆柱轴线上一点的初始值x01'，y01'以及底圆半径的初始值r0，再根据公式:Ｒx01y010=x0'1y0'10(6)将新坐标系下的x0'1、y0'1转换成旧坐标系下的坐标x01、y01。无论是在新旧坐标系下，底圆半径的初始值r0都不变。为了尽量减小误差，将圆柱轴线上一点的初始值z01的值设为原坐标系下圆柱点z的最大值与z的最小值和的一半。于是求出圆柱轴线上一点、圆柱底半径这四个参数的初始值x01、y01、z01、r0。(2)建立改进误差方程式求解参数值。圆柱面方程为:(x－x1)2+(y－y1)2+(z－z1)2－［a(x－x1)+b(y－y1)+c(z－z1)］2=r2(7)传统方法，令:f=(x－x1)2+(y－y1)2+(z－z1)2－［a(x－x1)+b(y－y1)+c(z－z1)］2－r2(8)对f线性化，得:f=f0+fx1x1+fy1y1+fz1z1+faa+fbb+fcc+frr(9)这样，误差方程可列为:V=f0+fx1x1+fy1y1+fz1z1+faa+fbb+fcc+frr(10)其中，f0=(x－x01)2+(y－y01)2+(z－z01)2－［a0(x－x01)+b0(y－y01)+c0(z－z01)］2－(r0)2(11)fx1=2*a0*［a0*(x－x01)+b0*(y－y01)+c0*(z－z01)］－2*(x－x01)(12)fa=－2*［a0*(x－x01)+b0*(y－y01)+c0*(z－z01)］*(x－x01)(13)fr=－2*r0(14)同理可得fy1、fz1、fb、fc。本文算法对误差方程式进行改进，令:f=(x－x1)2+(y－y1)2+(z－z1)2－［a(x－x1)+b(y－y1)+c(z－z1)］2－r2(15)这次令:f={(x－x1)2+(y－y1)2+(z－z1)2－［a(x－x1)+b(y－y1)+c(z－z1)］2－r2}/2r(16)对f线性化:f=f0+fx1x1+fy1y1+fz1z1+faa+fbb+fcc+frr(17)这样，误差方程可列为:V=f0+fx1x1+fy1y1+fz1z1+faa+fbb+fcc+frr(18)其中，f0={(x－x01)2+(y－y01)2+(z－z01)2－［a0(x－x01)+b0(y－y01)+c0(z－z01)］2－(r0)2}/2*r0(19)fx1={a0*［a0*(x－x01)+b0*(y－y01)+c0*(z－z01)］－(x－x01)}/r0(20)fa={－［a0*(x－x01)+b0*(y－y01)+c0*(z－z01)］*(x－x01)}/r0(21)fr=－{(x－x01)2+(y－y01)2+(z－z01)2－［a0*(x－x01)+b0*(y－y01)+c0*(z－z01)］2}/2*(r0)2－12(22)同理可得fy1、fz1、fb、fc。将误差方程写成矩阵形式，即:

2017年第12期工程勘察GeotechnicalInvestigation＆Surveying63Vn×1=Bn×7X7×1－Ln×1(23)假设获取的圆柱面上的点等精度，即权相等。根据最小二乘原理VTV=min，最后求得未知参数X［8，9］:X=(BTB)－1BTL(24)这是一个循环的迭代过程，每一次迭代过程中代入的初值都等于上一次的初值加上求出的X的改正值，当X的数值小到满足要求的精度时退出迭代。3算例分析本文采用的原始数据为某个桥墩柱面点云。图3为由点云数据组成的真实圆柱面;图4为利用传统方法依据柱面点云三维坐标拟合出的圆柱面;图5为本文算法依据柱面点云三维坐标拟合出的圆柱面。图3实验点云Fig.3Theexemplarypointclouds图4传统算法圆柱面拟合结果Fig.4Ｒesultsofthetraditionalcylindricalsurfacefittingmethod利用传统算法与上述设计算法，得出该圆柱面圆柱轴线上一点(x0，y0，z0)、圆柱轴线向量(a，b，c)、圆柱底圆半径r0这七个参数。如表1所示。传统算法拟合出的圆柱面方程为:(x+32.5293)2+(y－35.0874)2+(z+0.0007)2－［0.0010(x+32.5293)+0.0012(y－35.0874)+0.9987(z+0.0007)］2=0.66582图5本文算法圆柱面拟合结果Fig.5Ｒesultsofthecylindricalsurfacefittingmethodinthispaper表1柱面拟合参数结果Table1Ｒesultsofcylindricalfittingparameters参数传统算法(m)本文算法(m)x0－32.5293－32.5297y035.087435.0876z0－0.0007－0.0056a0.00100.0011b0.00120.0011c0.99870.9978r00.66580.6666本文算法拟合出的圆柱面方程为:(x+32.5297)2+(y－35.0876)2+(z+0.0056)2－［0.0011(x+32.5297)+0.0011(y－35.0876)+0.9978(z+0.0056)］2=0.66662误差方程式中的V代表残差，它可以作为评价圆柱面拟合质量好坏的标准。通过传统算法、本文算法拟合所得残差结果统计直方图如图6、图7所示。由于残差V数据量过大，所以截取少量V进行对比分析。采用传统算法与上述设计算法所得部分点云残差结果对比如表2所示。图6传统算法残差结果统计直方图Fig.6Statisticalhistogramoftheresidualerrorsusingtraditionalmethod对比图6、图7可得，采用传统方法所得结果残差值范围在±0.08m之间，而本文方法在±0.06m之间。如表2所示，相对于传统方法，采

64工程勘察GeotechnicalInvestigation＆Surveying2017年第12期图7本文算法残差结果统计直方图Fig.7Statisticalhistogramoftheresidualerrorsusingmethodinthispaper表2部分柱面拟合残差结果Table2Partresultsofcylindricalfittingresidualerror点号传统算法(m)本文算法(m)1－0.0250－0.01872－0.0248－0.01863－0.0204－0.01544－0.0245－0.01865－0.0203－0.01556－0.0248－0.01907－0.0257－0.01978－0.0251－0.01929－0.0267－0.020510－0.0249－0.0193用本文方法所得残差值更小。由运行结果可得，传统算法单位权方差σ2为3.4022×102mm2，残差V服从均值u为0，标准差σ为18.4450mm的正态分布;本文算法单位权方差σ2为1.9165×102mm2，残差V服从均值u为0，标准差σ为13.8437mm的正态分布。综上所述，与传统方法相比，本文算法的残差与标准差都小于传统方法，精度有所提高。4结语根据圆柱面的几何原理，本文探讨了一种改进的柱面拟合方法。它适用于空间中任意角度的圆柱(包括直圆柱和斜圆柱)，而且原理简单，易于实现。该方法利用主成分分析法和线性最小二乘法确定拟合初值，再利用非线性最小二乘法循环迭代求解圆柱模型参数，实现圆柱面拟合，与传统方法相比，本文算法的特点在于误差方程式的不同。就本文算法与传统方法对比可得，本文算法的残差明显小于传统方法，标准差σ从18.4450mm降为13.8437mm，精度大大提高。实验结果表明，在基于几何特性的圆柱面拟合过程中，对误差方程式稍作改进，可大大提高最终求解精度。该思想可以推广到下一步的研究工作中，有望进一步提高圆柱面拟合精度。参考文献［1］李广云，李宗春．工业测量系统原理与应用［M］．北京:测绘出版社，2011．LiGuangyun，LiZongchun．Principleandapplicationofindustrialmeasurementsystem［M］．Beijing:SurveyingandMappingPress，2011．(inChinese)［2］郑德华．三维激光扫描数据处理的理论和方法［D］．上海:同济大学，2005．ZhengDehua．Theoryandmethodof3Dlaserscanningdataprocessing［D］．Shanghai:TongjiUniversity，2005．(inChinese)［3］秦世伟，潘国荣，谷川等．基于遗传算法的三维空间柱面拟合［J］．同济大学学报(自然科学版)，2010，38(4):604～607．QinShiwei，PanGuorong，GuChuanetal．3Dcylindricalfittingbasedongeneticalgorithm［J］．JournalofTongjiUniversity(NaturalScience)，2010，38(4):604～607．(inChinese)［4］刘元朋，张定华．逆向工程中圆柱体几何特征参数评估方法的研究［J］．机械科学与技术，2005，24(3):310～311．LiuYuanpeng，ZhangDinghua．Ｒesearchonthemethodofevaluatingthecylinderalgeometriccharacteristicparametersinreverseengineering［J］．MechanicalScience＆Technology，2005，24(3):310～311．(inChinese)［5］王解先．工业测量中一种二次曲面的拟合方法［J］．武汉大学学报(信息科学版)，2007，32(1):47～50．WangJiexian．Aquadricsurfacefittingmethodinindustrialmeasurement［J］．GeomaticsandInformationScienceofWuhanUniversity，2007，32(1):47～50．(inChinese)［6］王崇倡，王晓婉，徐晓昶．圆柱面拟合方法研究［J］．测绘工程，2014，23(3):5～9．WangChongchang，WangXiaowan，XuXiaochang．Studyoncylindricalsurfacefittingmethod［J］．EngineeringofSurveyingandMapping，2014，23(3):5～9．(inChinese)［7］武汉大学测绘学院测量平差组．误差理论与测量平差基础［M］．武汉:武汉大学出版社，2013．SurveyAdjustmentGroupofSurveyingandMappingInstituteofWuhanUniversity．Errortheoryandsurveyadjustmentfoundation［M］．Wuhan:WuhanUniversityPress，2013．(inChinese)［8］陈佩．主成分分析法研究及其在特征提取中的应用［D］．西安:陕西师范大学，2014．ChenPei．Principalcomponentanalysisanditsapplicationinfeatureextraction［D］．Xi’an:ShanxiNormalUniversity，2014．(inChinese)［9］唐利民，朱建军．非线性最小二乘问题数值迭代法的统一模型及其不适定性［J］．长沙交通学院学报，2008，24(2):39～43．TangLimin，ZhuJianjun．Unifiedmodelanditsillfitofnumericaliterativemethodfornonlinearleastsquareproblems［J］．JournalofChangshaJiaotongCollege，2008，24(2):39～43．(inChinese)