基于粒子群優化算法的圓柱體擬合方法

袁 輝

（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算法的圓柱體擬合方法，用于解決工程測量中圓柱體模型參數擬合問題。實驗結果表明，該方法的擬合精度與最小二乘法相當，但其實現更簡單，也避免了參數初值選擇問題，具有有效性和實用性。在本文實驗中，慣性權重、學習因子、比例系數等參數均采用已有文獻中的推薦值，在后續實驗中需根據具體問題進行更準確的確定，使算法 收斂更快。