基于粒子群優(yōu)化算法的圓柱體擬合方法

袁 輝

（1.中鐵第四勘察設(shè)計(jì)院集團(tuán)有限公司，湖北 武漢 430063）

在工程測(cè)量中，經(jīng)常涉及圓柱體擬合問題，如支柱的半徑測(cè)量與垂直度檢測(cè)、管道變形監(jiān)測(cè)與姿態(tài)測(cè)量、隧道變形安全監(jiān)測(cè)等。最小二乘法是解決圓柱體擬合問題的常用方法[1-4]，由于圓柱體方程是非線性的，因此最小二乘法的實(shí)現(xiàn)復(fù)雜，且擬合結(jié)果受初值選取影響很大。袁建剛[5]等利用主成分分析法和線性最小二乘法計(jì)算得到圓柱體參數(shù)的擬合初值，再利用非線性最小二乘法循環(huán)迭代求解圓柱體模型參數(shù)，該方法過程繁瑣且計(jì)算量巨大；申旭[6]等提出了一種圓柱體擬合的參數(shù)初值查找算法，但該算法受柱面點(diǎn)云分布和噪聲點(diǎn)影響較大。目前也有基于智能算法的圓柱體擬合方法研究，如利用遺傳算法擬合得到圓柱體的模型參數(shù)[7-8]，但遺傳算法容易早熟收斂；粒子群優(yōu)化（PSO）算法也是一種智能算法[9-10]，相較于遺傳算法，其原理更簡單、收斂速度更快，且無需交叉變異操作，參數(shù)更少，實(shí)現(xiàn)更方便。因此，本文提出了基于PSO算法的圓柱體擬合方法。

1 圓柱體模型

在三維空間中，圓柱體是到中心軸線距離為某一常數(shù)的點(diǎn)集合，如圖1所示。

圖1 圓柱體模型

圓柱體包括7個(gè)參數(shù)，分別為圓柱體的中軸線約束方向向量s(a,b,c)、該中軸線上某一點(diǎn)的坐標(biāo)p0(x0,y0,z0)和圓柱體半徑r。根據(jù)圓柱體的定義，其方程可表示為：

約束方向向量s為單位向量，即a2+b2+c2=1，因此式（1）可簡化為：

空間中一點(diǎn)pj(xj,yj,zj)，其誤差方程可表示為：

因此，擬合的目標(biāo)函數(shù)為：

式中，N為擬合點(diǎn)的個(gè)數(shù)。

2 PSO算法流程

PSO算法[9]源于對(duì)鳥群捕食的行為研究，通過群體中個(gè)體之間的協(xié)作和信息共享來尋找全局最優(yōu)解。該算法利用粒子來模擬鳥群中的鳥，每個(gè)粒子從隨機(jī)解出發(fā)單獨(dú)尋找最優(yōu)解，將粒子個(gè)體的最優(yōu)解與群體其他粒子共享找到群體當(dāng)前最優(yōu)解；然后所有粒子根據(jù)自身最優(yōu)解與群體最優(yōu)解來調(diào)整自身的速度和位置，多次迭代后搜尋到群體最優(yōu)解。PSO算法流程如圖2所示。

圖2 PSO算法流程圖

本文利用PSO算法對(duì)圓柱體進(jìn)行擬合，具體過程為：定義粒子由圓柱體的7個(gè)參數(shù)組成，粒子i為oi=(xi,yi,zi,ri,ai,bi,ci)，其中i∈[1,n]；n為粒子群中粒子總數(shù)；算法迭代次數(shù)限制為kmax，第k∈[1,kmax]次迭代時(shí)，粒子位置為粒子的 每個(gè)維度值均需在給定值域范圍內(nèi)，即

式中，xmin和xmax、ymin和ymax、zmin和zmax、rmin和rmax、amin和amax、bmin和bmax、cmin和cmax分別為圓柱體7個(gè)參數(shù)值域的上下限。

通常xmin和xmax、ymin和ymax、zmin和zmax的取值可為點(diǎn)云數(shù)據(jù)的空間范圍，即

半徑參數(shù)rmin=0，rmax根據(jù)點(diǎn)云在x、y、z方向上的長度確定，rmax=max{xmax-xmin,ymax-ymin,zmax-zmin}。由于約束方向向量s(a,b,c)為單位向量，進(jìn)一步約束s(a,b,c)方向朝上，即c≥0，因此amin=-1、amax=1、bmin=-1、bmax=1、cmin=0、cmax=1。

粒子i在第k次迭代時(shí)的速度為同樣地，每個(gè)維度的速度值也要求在給定的值域范圍內(nèi)，用來限制粒子探索解空間的步幅大小。粒子各維度參數(shù)的速度限制為其中分別為粒子各維度參數(shù)的速度最大值。若速度最大值過大，則粒子可能飛過最優(yōu)位置；若速度最大值過小，則粒子探索范圍較小，可能陷入局部最優(yōu)。通常每個(gè)維度的最大速度為每維變化范圍的10%～20%[10]，因此定義該比例系數(shù)為β，則有

粒子i在第k次迭代的適應(yīng)度計(jì)算公式為：

本文中粒子適應(yīng)度越低，粒子位置越優(yōu)。粒子i在迭代過程中的最優(yōu)位置為obesti，對(duì)應(yīng)的適應(yīng)度為Fmini；粒子群的全局最優(yōu)位置為gbest，對(duì)應(yīng)的適應(yīng)度為Fmgin。粒子i在第k次迭代過程中的運(yùn)動(dòng)速度為：

式中，wk為慣性權(quán)重，體現(xiàn)了粒子上次速度對(duì)當(dāng)前速度的影響；λ1、λ2為學(xué)習(xí)因子，表示粒子的加速度，用于調(diào)節(jié)學(xué)習(xí)步長；γ1∈[0,1]、γ2∈[0,1]為兩個(gè)隨機(jī)數(shù)，用于增加搜索隨機(jī)性。

式（8）中第一部分是粒子的慣性速度；第二部分是粒子的“認(rèn)知”，表示粒子本身的思考，可理解為粒子當(dāng)前位置與自身最優(yōu)位置之間的距離；第三部分是粒子的“社會(huì)”，表示粒子間信息共享與合作，可理解為粒子當(dāng)前位置與全局最優(yōu)位置之間的距離[9,11]。

在全局搜索過程中，通常前期需要粒子具有較強(qiáng)的探索能力，以快速探索到全局最優(yōu)位置附近，而后期則需要粒子具有較強(qiáng)的局部搜索能力，以加快收斂速度，因此wk可設(shè)置為隨迭代次數(shù)線性減小[12]，即

式中，wmin、wmax分別為慣性權(quán)重最小、最大值。

對(duì)于更新的粒子速度，按照各維度最大速度值進(jìn)行約束改正，再基于粒子運(yùn)動(dòng)速度更新粒子的當(dāng)前位置，則有：

根據(jù)相應(yīng)的閾值范圍，對(duì)更新的粒子位置的每個(gè)維度值進(jìn)行約束改正。基于PSO算法的圓柱體擬合方法的詳細(xì)步驟為：

1）定義粒子群規(guī)模n，隨機(jī)生成粒子?i∈[1,n]的初始位置和初始速度=粒子i的最優(yōu)位置=粒子的最優(yōu)適應(yīng)度。粒子群全局最優(yōu)適應(yīng)度，粒子群全局最優(yōu)位置gbest為對(duì)應(yīng)的粒子位置。設(shè)置學(xué)習(xí)因子λ1、λ2，慣性權(quán)重wmin、wmax，比例系數(shù)β等參數(shù)值，給定最大迭代次數(shù)kmax，迭代次數(shù)k=1。

2）根據(jù)式（7）計(jì)算每個(gè)粒子的適應(yīng)度。

4）根據(jù)式（8）更新粒子速度，根據(jù)式（10）更新粒子位置。

5）迭代次數(shù)k=k+1，若全局最優(yōu)適應(yīng)度Fming足夠好或k≥kmax則停止迭代，否則轉(zhuǎn)到步驟2）。

3 實(shí)驗(yàn)分析

為了驗(yàn)證基于PSO算法的圓柱體擬合方法的有效性，本文將其與最小二乘法進(jìn)行對(duì)比試驗(yàn)。實(shí)驗(yàn)數(shù)據(jù)為徠卡P40三維激光掃描儀采集的武漢地鐵5號(hào)線某區(qū)間隧道的點(diǎn)云數(shù)據(jù)。點(diǎn)云數(shù)據(jù)經(jīng)過重采樣、去噪處理，從中截取得到某一環(huán)片上的點(diǎn)云如圖3所示，該環(huán)片點(diǎn)云約有35 000個(gè)點(diǎn)、點(diǎn)間距為3 cm，環(huán)片長度為1.6 m，隧道設(shè)計(jì)直徑為5.5 m。

圖3 隧道環(huán)片點(diǎn)云數(shù)據(jù)

基于Visual Studio 2015，采用C#編程語言實(shí)現(xiàn)了本文算法和最小二乘法。實(shí)驗(yàn)平臺(tái)配置為Intel Core i7 3.4 GHz CPU、8 GB內(nèi)存、Windows 7（64位）操作 系統(tǒng)。在本文算法中，粒子群規(guī)模n=50，學(xué)習(xí)因子λ1=2、λ2=2[13]，慣性權(quán)重wmin=0.4、wmax=0.9[13]，比例系數(shù)β=0.2，最大迭代次數(shù)kmax=100。

圓柱體模型參數(shù)擬合結(jié)果如表1所示，可以看出，對(duì)于圓柱體中軸線上點(diǎn)p0(x0,y0,z0)的坐標(biāo)，本文算法與最小二乘法選取的位置不同，為了方便比較，將最小二乘法擬合的p0沿中軸線平移到z0=10.411 9，平移后的x0、y0分別為529 806.760 8、383 858.974 8，因此x0、y0與本文算法擬合結(jié)果的差值分別為0.5 mm和0.4 mm；兩種算法擬合得到的圓柱體半徑分別為2.749 8 m和2.750 1 m，差值為0.3 mm，差值很小；本文算法擬合的中軸線與最小二乘法擬合的中軸線夾角為0.01e，差值也非常小。因此，本文算法的擬合結(jié)果雖為近似解，但收斂后的結(jié)果能滿足實(shí)際工程測(cè)量的精度要求。

表1 圓柱體模型參數(shù)擬合結(jié)果

由于PSO算法具有隨機(jī)搜索過程，算法每次收斂迭代次數(shù)不同，因此本文進(jìn)行了10次重復(fù)試驗(yàn)，平均收斂迭代次數(shù)為74.6，平均運(yùn)行耗時(shí)為55.19 s。以其中一次典型的實(shí)驗(yàn)過程為例，其全局最優(yōu)適應(yīng)度變化曲線如圖4所示，可以看出，由于算法在粒子群初始化階段都是隨機(jī)位置，因此初始適應(yīng)度很大，在前 20次迭代過程中，全局最優(yōu)適應(yīng)度迅速減小，說明算法收斂速度較快；經(jīng)過約40次迭代后，粒子已探索到全局最優(yōu)位置附近；后續(xù)迭代進(jìn)行局部搜索優(yōu)化，不斷提高全局最優(yōu)解的精度。由此可見，本文算法的收斂速度較快，結(jié)果精度較高，能有效解決圓柱體擬合問題。

圖4 全局最優(yōu)適應(yīng)度變化曲線

4 結(jié) 語

本文提出了基于PSO算法的圓柱體擬合方法，用于解決工程測(cè)量中圓柱體模型參數(shù)擬合問題。實(shí)驗(yàn)結(jié)果表明，該方法的擬合精度與最小二乘法相當(dāng)，但其實(shí)現(xiàn)更簡單，也避免了參數(shù)初值選擇問題，具有有效性和實(shí)用性。在本文實(shí)驗(yàn)中，慣性權(quán)重、學(xué)習(xí)因子、比例系數(shù)等參數(shù)均采用已有文獻(xiàn)中的推薦值，在后續(xù)實(shí)驗(yàn)中需根據(jù)具體問題進(jìn)行更準(zhǔn)確的確定，使算法 收斂更快。