基于粒子群优化算法的圆柱体拟合方法

袁 辉

（1.中铁第四勘察设计院集团有限公司，湖北 武汉 430063）

在工程测量中，经常涉及圆柱体拟合问题，如支柱的半径测量与垂直度检测、管道变形监测与姿态测量、隧道变形安全监测等。最小二乘法是解决圆柱体拟合问题的常用方法[1-4]，由于圆柱体方程是非线性的，因此最小二乘法的实现复杂，且拟合结果受初值选取影响很大。袁建刚[5]等利用主成分分析法和线性最小二乘法计算得到圆柱体参数的拟合初值，再利用非线性最小二乘法循环迭代求解圆柱体模型参数，该方法过程繁琐且计算量巨大；申旭[6]等提出了一种圆柱体拟合的参数初值查找算法，但该算法受柱面点云分布和噪声点影响较大。目前也有基于智能算法的圆柱体拟合方法研究，如利用遗传算法拟合得到圆柱体的模型参数[7-8]，但遗传算法容易早熟收敛；粒子群优化（PSO）算法也是一种智能算法[9-10]，相较于遗传算法，其原理更简单、收敛速度更快，且无需交叉变异操作，参数更少，实现更方便。因此，本文提出了基于PSO算法的圆柱体拟合方法。

1 圆柱体模型

在三维空间中，圆柱体是到中心轴线距离为某一常数的点集合，如图1所示。

图1 圆柱体模型

圆柱体包括7个参数，分别为圆柱体的中轴线约束方向向量s(a,b,c)、该中轴线上某一点的坐标p0(x0,y0,z0)和圆柱体半径r。根据圆柱体的定义，其方程可表示为：

约束方向向量s为单位向量，即a2+b2+c2=1，因此式（1）可简化为：

空间中一点pj(xj,yj,zj)，其误差方程可表示为：

因此，拟合的目标函数为：

式中，N为拟合点的个数。

2 PSO算法流程

PSO算法[9]源于对鸟群捕食的行为研究，通过群体中个体之间的协作和信息共享来寻找全局最优解。该算法利用粒子来模拟鸟群中的鸟，每个粒子从随机解出发单独寻找最优解，将粒子个体的最优解与群体其他粒子共享找到群体当前最优解；然后所有粒子根据自身最优解与群体最优解来调整自身的速度和位置，多次迭代后搜寻到群体最优解。PSO算法流程如图2所示。

图2 PSO算法流程图

本文利用PSO算法对圆柱体进行拟合，具体过程为：定义粒子由圆柱体的7个参数组成，粒子i为oi=(xi,yi,zi,ri,ai,bi,ci)，其中i∈[1,n]；n为粒子群中粒子总数；算法迭代次数限制为kmax，第k∈[1,kmax]次迭代时，粒子位置为粒子的 每个维度值均需在给定值域范围内，即

式中，xmin和xmax、ymin和ymax、zmin和zmax、rmin和rmax、amin和amax、bmin和bmax、cmin和cmax分别为圆柱体7个参数值域的上下限。

通常xmin和xmax、ymin和ymax、zmin和zmax的取值可为点云数据的空间范围，即

半径参数rmin=0，rmax根据点云在x、y、z方向上的长度确定，rmax=max{xmax-xmin,ymax-ymin,zmax-zmin}。由于约束方向向量s(a,b,c)为单位向量，进一步约束s(a,b,c)方向朝上，即c≥0，因此amin=-1、amax=1、bmin=-1、bmax=1、cmin=0、cmax=1。

粒子i在第k次迭代时的速度为同样地，每个维度的速度值也要求在给定的值域范围内，用来限制粒子探索解空间的步幅大小。粒子各维度参数的速度限制为其中分别为粒子各维度参数的速度最大值。若速度最大值过大，则粒子可能飞过最优位置；若速度最大值过小，则粒子探索范围较小，可能陷入局部最优。通常每个维度的最大速度为每维变化范围的10%～20%[10]，因此定义该比例系数为β，则有

粒子i在第k次迭代的适应度计算公式为：

本文中粒子适应度越低，粒子位置越优。粒子i在迭代过程中的最优位置为obesti，对应的适应度为Fmini；粒子群的全局最优位置为gbest，对应的适应度为Fmgin。粒子i在第k次迭代过程中的运动速度为：

式中，wk为惯性权重，体现了粒子上次速度对当前速度的影响；λ1、λ2为学习因子，表示粒子的加速度，用于调节学习步长；γ1∈[0,1]、γ2∈[0,1]为两个随机数，用于增加搜索随机性。

式（8）中第一部分是粒子的惯性速度；第二部分是粒子的“认知”，表示粒子本身的思考，可理解为粒子当前位置与自身最优位置之间的距离；第三部分是粒子的“社会”，表示粒子间信息共享与合作，可理解为粒子当前位置与全局最优位置之间的距离[9,11]。

在全局搜索过程中，通常前期需要粒子具有较强的探索能力，以快速探索到全局最优位置附近，而后期则需要粒子具有较强的局部搜索能力，以加快收敛速度，因此wk可设置为随迭代次数线性减小[12]，即

式中，wmin、wmax分别为惯性权重最小、最大值。

对于更新的粒子速度，按照各维度最大速度值进行约束改正，再基于粒子运动速度更新粒子的当前位置，则有：

根据相应的阈值范围，对更新的粒子位置的每个维度值进行约束改正。基于PSO算法的圆柱体拟合方法的详细步骤为：

1）定义粒子群规模n，随机生成粒子∀i∈[1,n]的初始位置和初始速度=粒子i的最优位置=粒子的最优适应度。粒子群全局最优适应度，粒子群全局最优位置gbest为对应的粒子位置。设置学习因子λ1、λ2，惯性权重wmin、wmax，比例系数β等参数值，给定最大迭代次数kmax，迭代次数k=1。

2）根据式（7）计算每个粒子的适应度。

4）根据式（8）更新粒子速度，根据式（10）更新粒子位置。

5）迭代次数k=k+1，若全局最优适应度Fming足够好或k≥kmax则停止迭代，否则转到步骤2）。

3 实验分析

为了验证基于PSO算法的圆柱体拟合方法的有效性，本文将其与最小二乘法进行对比试验。实验数据为徕卡P40三维激光扫描仪采集的武汉地铁5号线某区间隧道的点云数据。点云数据经过重采样、去噪处理，从中截取得到某一环片上的点云如图3所示，该环片点云约有35 000个点、点间距为3 cm，环片长度为1.6 m，隧道设计直径为5.5 m。

图3 隧道环片点云数据

基于Visual Studio 2015，采用C#编程语言实现了本文算法和最小二乘法。实验平台配置为Intel Core i7 3.4 GHz CPU、8 GB内存、Windows 7（64位）操作 系统。在本文算法中，粒子群规模n=50，学习因子λ1=2、λ2=2[13]，惯性权重wmin=0.4、wmax=0.9[13]，比例系数β=0.2，最大迭代次数kmax=100。

圆柱体模型参数拟合结果如表1所示，可以看出，对于圆柱体中轴线上点p0(x0,y0,z0)的坐标，本文算法与最小二乘法选取的位置不同，为了方便比较，将最小二乘法拟合的p0沿中轴线平移到z0=10.411 9，平移后的x0、y0分别为529 806.760 8、383 858.974 8，因此x0、y0与本文算法拟合结果的差值分别为0.5 mm和0.4 mm；两种算法拟合得到的圆柱体半径分别为2.749 8 m和2.750 1 m，差值为0.3 mm，差值很小；本文算法拟合的中轴线与最小二乘法拟合的中轴线夹角为0.01e，差值也非常小。因此，本文算法的拟合结果虽为近似解，但收敛后的结果能满足实际工程测量的精度要求。

表1 圆柱体模型参数拟合结果

由于PSO算法具有随机搜索过程，算法每次收敛迭代次数不同，因此本文进行了10次重复试验，平均收敛迭代次数为74.6，平均运行耗时为55.19 s。以其中一次典型的实验过程为例，其全局最优适应度变化曲线如图4所示，可以看出，由于算法在粒子群初始化阶段都是随机位置，因此初始适应度很大，在前 20次迭代过程中，全局最优适应度迅速减小，说明算法收敛速度较快；经过约40次迭代后，粒子已探索到全局最优位置附近；后续迭代进行局部搜索优化，不断提高全局最优解的精度。由此可见，本文算法的收敛速度较快，结果精度较高，能有效解决圆柱体拟合问题。

图4 全局最优适应度变化曲线

4 结 语

本文提出了基于PSO算法的圆柱体拟合方法，用于解决工程测量中圆柱体模型参数拟合问题。实验结果表明，该方法的拟合精度与最小二乘法相当，但其实现更简单，也避免了参数初值选择问题，具有有效性和实用性。在本文实验中，惯性权重、学习因子、比例系数等参数均采用已有文献中的推荐值，在后续实验中需根据具体问题进行更准确的确定，使算法 收敛更快。