改進差分進化算法求解B樣條曲線曲面擬合問題

李小艷　陳紹平

(武漢理工大學(xué)理學(xué)院　武漢 湖北 430000)

0　引　言

隨著計算機運行速度的提高和逆向工程的快速發(fā)展，從20世紀90年代中期起，B樣條曲線曲面擬合一直是CAGD領(lǐng)域的研究熱點，也是CAGD領(lǐng)域的經(jīng)典問題之一[1]。B樣條具有幾何變換不變性、局部支撐性、凸包性、變差縮減性等諸多優(yōu)良性質(zhì)，這使得它成為形狀數(shù)學(xué)描述的主流方法，并廣泛應(yīng)用于數(shù)據(jù)可視化與逼近、醫(yī)學(xué)圖像輪廓重構(gòu)和幾何建模等領(lǐng)域[2]。

解決B樣條曲線曲面擬合問題的關(guān)鍵在于B樣條節(jié)點矢量和控制頂點的選取，合適的節(jié)點矢量與控制頂點能保證擬合曲線的光滑度和擬合效果。把節(jié)點矢量和控制頂點均視為變量時，擬合問題就變?yōu)橐粋€多維多變量高度非線性的最優(yōu)化問題。解決該問題的方法大致可以被分為三類：傳統(tǒng)的節(jié)點插入和節(jié)點刪除方法、支配點和特征點方法以及人工智能算法。傳統(tǒng)的方法一般是將節(jié)點向量進行固定，得到一個關(guān)于控制頂點的線性系統(tǒng)，但是這些方法并不能處理帶有不連續(xù)點、尖點和拐點的復(fù)雜型曲線。支配點和特征點方法是通過提取數(shù)據(jù)點的形狀信息(包括拐點、折痕點和曲率極值點等)用于節(jié)點的選擇，但是這些方法往往僅限于連續(xù)曲線。近年來，利用人工智能算法求解并優(yōu)化B樣條擬合問題受到諸多學(xué)者的關(guān)注。如Yoshimoto等[3]提出了一種實數(shù)編碼的遺傳算法對平面曲線進行擬合，該方法能自動的確定節(jié)點的數(shù)量和位置，不僅可以處理光滑的數(shù)據(jù)擬合問題，也可以處理帶尖點和不連續(xù)點的數(shù)據(jù)擬合問題；周明華等[4]采用遺傳算法同時對數(shù)據(jù)點的參數(shù)化和節(jié)點矢量進行優(yōu)化，與傳統(tǒng)方法相比取得了較高的擬合效果；穆國旺等[5]提出一種改進的遺傳算法，將誤差精度考慮在內(nèi)，實現(xiàn)了在給定誤差下能使用最少的控制頂點進行B樣條曲線擬合，但是該方法的運算時間較長；Adi等[6]采用粒子群優(yōu)化(PSO)方法進行NURBS曲線擬合，該方法將控制頂點作為待求解的變量，但是誤差比較大；Galvez等[7]利用PSO通過固定節(jié)點個數(shù)自適應(yīng)的確定節(jié)點的位置優(yōu)化擬合曲線；Ulker等[8]提出將人工免疫算法(AIS)應(yīng)用于求解B樣條曲線擬合問題上，得到較好的結(jié)果；Zhao等[9]在B樣條曲線逼近問題中采用了一種新穎的隨機優(yōu)化方法，利用高斯混合模型和聚類技術(shù)來實現(xiàn)自適應(yīng)節(jié)點配置，但該方法僅對封閉曲線進行了擬合實驗。何兵朋等[10]對差分進化算法的變異操作進行改進，設(shè)計出一種新的B樣條曲線擬合方法。近幾年，許多學(xué)者也不斷的對這些智能算法進行改進來改善B樣條擬合效果。如Trejo-Caballero G 等[11]設(shè)計出一種分層編碼方法，將節(jié)點與控制頂點分別進行二進制編碼和實數(shù)編碼，利用遺傳算法進行擬合，得到了較好的擬合效果。

在文獻[13]中，他們對DE、PSO和EA在處理數(shù)值優(yōu)化問題時的表現(xiàn)進行了評估，與其他方法相比，DE展示出了較高的優(yōu)越性。因此本文在諸多學(xué)者研究的基礎(chǔ)上，將一種改進的自適應(yīng)差分進化算法應(yīng)用于解決B樣條曲線曲面最小二乘擬合問題。此方法采用了DE/current-to-best/的變異策略和帶有隨機游走(Random Walk)的交叉操作，利用隨機游走的隨機機制彌補DE/current-to-best/容易早熟收斂的缺陷。同時為了跳出局部最優(yōu)解，利用混沌系統(tǒng)的隨機性與遍歷性設(shè)計出一個針對當前種群最優(yōu)解的局部搜索操作，有效地平衡了DE的開發(fā)能力和探索能力。新的差分進化算法在進行B樣條曲線曲面擬合時能夠根據(jù)數(shù)據(jù)點分布特征配置內(nèi)節(jié)點的位置，并在需要時產(chǎn)生多重節(jié)點。試驗結(jié)果表示，與標準的差分進化算法相比，用本文的方法得到的B樣條曲線曲面逼近精度更高。

1　B樣條曲線的最小二乘擬合

(1)

式中：j=0,1,…,M-1;xj=[a,b]，f(x)為采樣函數(shù)(對于散亂數(shù)據(jù)來說，f(x)是未知的)，ε為噪聲誤差，N為采樣數(shù)量，n為內(nèi)節(jié)點個數(shù)，[a,b]為擬合區(qū)間。k次B樣條目標曲線為：

(2)

(3)

式中：pi(i=0,1,…,n+k)為控制頂點，Ni,k(x)(i=0,1,…,n+k)為由de-Boor遞推公式定義的k次B樣條基函數(shù)(本文選擇三次B樣條曲線)：

‖·‖是歐氏距離。端節(jié)點設(shè)置為k+1重節(jié)點：u0=u1=…=uk,un+k+1=un+k+2=…=un+2k+1對于一組固定的內(nèi)節(jié)點，采用標準的最小二乘技術(shù)，欲使式(3)中的Q達到最小，應(yīng)使它關(guān)于n+k-1個控制頂點的導(dǎo)數(shù)為0，于是就得到了n+k-1個以pi(i=1,2,…,n+k-1)為未知量的線性方程組：

(NTN)P=NTR

(4)

這里N是(M-2)×(n+k-1)的標量矩陣：

2　差分進化算法原理

差分進化算法DE(Differential Evolution)與其他遺傳類型算法相似，以迭代的方式對種群中個體進行變異、交叉和選擇操作直到根據(jù)適應(yīng)度函數(shù)找出最佳個體。標準的DE算法描述如下：

2.1　差分變異操作

變異就是指從當前種群中隨機選擇幾個個體得到差分矢量乘上變異系數(shù)作用于目標個體對其產(chǎn)生擾動，生成變異個體。五種最常用的差分變異策略如下：

DE/current-to-best/:

(5)

2.2　交叉操作

式中：CR∈[0，1],rand1(0,1)為均勻分布的隨機數(shù)，jrand為介于1到D的隨機整數(shù)下標索引，可以確保試驗個體中至少有一個基因是由變異個體貢獻的。

2.3　選擇操作

DE根據(jù)個體的適應(yīng)度函數(shù)值對目標個體和試驗個體進行貪婪選擇，以最小化問題為例，具體操作如下：

(6)

3　改進差分進化算法的B樣條擬合

(1) 差分變異操作。選擇式(5) DE/current-to-best/的變異策略，該策略具有局部搜索能力強，收斂速度快的優(yōu)點，但全局搜索能力較弱，容易陷入局部最優(yōu)。縮放因子F由下式?jīng)Q定：

(2) 交叉操作。算法早熟收斂是因為隨著DE的進化，種群中個體漸漸聚集，差異減小。為了避免局部最優(yōu)，本文借鑒文獻[19]中的帶有隨機游走的交叉操作，隨機游走是一個可以向任何方向前進的、到達任何位置的無規(guī)則運動，這種隨機性機制在用于DE搜索過程時可以增加種群的多樣性，有效避免早熟現(xiàn)象。交叉操作的表達式如下：

(7)

交叉概率CR和參數(shù)RW的表達式如下：

在DE中差分矢量縮放因子F和交叉概率CR對優(yōu)化性能有著重要的影響。Storn等[15]建議過參數(shù)的選取范圍應(yīng)為:F∈[0.5,1]、CR∈[0.8,1]，本文對參數(shù)的設(shè)置符合這個條件。

(3) 混沌局部搜索操作。設(shè)計出一種混沌局部搜索操作，在當前最優(yōu)個體的鄰域范圍內(nèi)搜索更優(yōu)解，以跳出局部最優(yōu)。通過一維Logistic映射生成混沌序列，迭代公式如下：

Zi+1=μZi(1-Zi)i=0,1,2,…

(8)

(9)

式中：λ為局部搜索系數(shù)，針對不同的實例λ取不同的值。若在搜索結(jié)果中找到了優(yōu)于當前代最優(yōu)個體的內(nèi)節(jié)點向量，則對其進行取代。

(4) 適應(yīng)度函數(shù)。使用統(tǒng)計學(xué)中貝葉斯信息準則BIC[23]作為適應(yīng)度函數(shù)，表達式為：

BIC=M×(lnQ)+(lnM)×(2n+k+1)

(10)

式中:M為數(shù)據(jù)點數(shù)量，Q為由式(3)計算所得的殘差平方和，n+k+1代表控制頂點個數(shù)，n代表內(nèi)節(jié)點個數(shù)。適應(yīng)度函數(shù)值越小，代表找到的節(jié)點矢量越優(yōu)，擬合誤差越小。

在B樣條曲線擬合問題中，節(jié)點數(shù)量是變化的，控制頂點的個數(shù)也會隨之改變，如果控制頂點太少，無法還原出原數(shù)據(jù)點分布特征，擬合誤差就會增大；控制頂點太多，擬合曲線的光滑度會降低。因此如果只以殘差平方和Q為目標函數(shù)，就會忽略模型復(fù)雜度和控制頂點個數(shù)對曲線幾何形狀的影響。選取BIC為適應(yīng)度函數(shù)，它的第一項能夠保證模型函數(shù)的精度，第二項添加了對模型函數(shù)自由參數(shù)個數(shù)的懲罰，通過精度與計算復(fù)雜度之間的平衡，獲得最優(yōu)的逼近模型。BIC的另一個優(yōu)勢在于它避免了使用主觀參數(shù)(如誤差界限和平滑因子)來判斷擬合模型。

基于該算法的三次B樣條曲線擬合流程如下：

輸出：最佳內(nèi)節(jié)點。

Step2利用式(8)求出種群中每個內(nèi)節(jié)點個體的適應(yīng)度值，獲得最優(yōu)個體及其對應(yīng)的適應(yīng)度值。

Step3判斷迭代是否達到Tmax，若是，退出，否則執(zhí)行下一步。

Step4對種群中個體根據(jù)式(5)和式(7)進行變異和交叉操作。

Step5按照式(6)進行選擇操作，生成新一代個體，令t=t+1。

Step6判斷種群是否陷入局部最優(yōu)，若是，則根據(jù)式(8)和式(9)進行混沌局部搜索，更新種群中最優(yōu)個體。若否，返回Step2。

4　仿真實驗與分析

為了驗證本文算法在解決B樣條曲線曲面最小二乘擬合問題上的有效性，選擇5個測試函數(shù)進行試驗分析。

4.1　曲線擬合

所有曲線擬合的試驗中，采樣數(shù)據(jù)點均由測試函數(shù)根據(jù)式(1)得到，噪聲誤差εj取均值為0，方差為1的標準正太分布隨機數(shù)。差分進化算法的參數(shù)均設(shè)置為種群規(guī)模NP=100，最大迭代次數(shù)Tmax=200，由于本文實驗擬合區(qū)間均為[0,1],局部搜索系數(shù)我們?nèi)≥^小的值λ=0.025。對于每一個測試函數(shù)我們都在取值范圍內(nèi)選取M=201個數(shù)據(jù)點，試驗在同等條件下重復(fù)30次。

該函數(shù)是多峰值的分段函數(shù)，x=0.6是它的不連續(xù)點。

圖1(a)、圖2(a)、圖3(a)分別為本文算法對取自函數(shù)f1(x)、f2(x)、f3(x)的數(shù)據(jù)點進行試驗得到的BIC值與內(nèi)節(jié)點數(shù)之間的關(guān)系，從圖中可以得到f1(x)、f2(x)、f3(x)對應(yīng)的最佳擬合內(nèi)節(jié)點個數(shù)分別為4、8、5。圖1-圖3中的(b)顯示了BIC值與殘差平方和Q的變化情況，實線表示適應(yīng)度BIC，虛線表示殘差平方和Q；圖1-圖3中的(c)分別繪出了例1-例3由最佳內(nèi)節(jié)點得到的最佳擬合效果圖，三角形標出了最佳內(nèi)節(jié)點的位置，×是擬合數(shù)據(jù)點，實線是B樣條擬合曲線。

圖1　例1函數(shù)擬合結(jié)果

圖2　例2函數(shù)擬合結(jié)果

圖3　例3函數(shù)擬合結(jié)果

圖1(b)、圖2(b)、圖3(b)為本文算法對例1-例3進行擬合過程中隨著迭代次數(shù)的增加，最佳個體的適應(yīng)度值和殘差平方和的變化曲線。從圖中可以看出，對于每一個測試函數(shù)，適應(yīng)度值基本都在在100代左右收斂，說明我們的算法收斂速度快且結(jié)果穩(wěn)定。圖1-圖3中(c)對應(yīng)例1-例3的擬合效果圖與文獻[3,10,12]中的擬合結(jié)果一致，并且擬合效果較好。必須說明的是對于三次B樣條曲線，如果某個內(nèi)節(jié)點的重復(fù)度為4，就說明該三次B樣條曲線在該點處是不連續(xù)的，函數(shù)f2(x)的擬合結(jié)果就說明了這種情況。這表示本文的算法能夠根據(jù)數(shù)據(jù)點的分布特征來自動分配節(jié)點位置，使擬合結(jié)果與數(shù)據(jù)點的分布特征一致。三個擬合效果圖可以說明本文算法不僅能夠較好地擬合如例1這種在某點處發(fā)生突變的光滑曲線，也能夠處理例2、例3這種帶有間斷和尖點的情況。例2的函數(shù)圖像在x=0.6處發(fā)生間斷，我們的方法相應(yīng)的在該點處產(chǎn)生了四重節(jié)點，根據(jù)B樣條曲線的性質(zhì)，擬合曲線在x=0.6處也發(fā)生了間斷。

將傳統(tǒng)的帶五種不同變異操作的差分進化算法分別與本文方法進行比較。進行30次試驗，表1、表2給出的是在相同節(jié)點數(shù)量的情況下對應(yīng)于例1-例3的試驗結(jié)果，總結(jié)的是30次試驗的最小、最大以及平均的適應(yīng)度值(BIC)和殘差平方和(Q)。從表1、表2中可以看出，由本文方法找到的最佳節(jié)點矢量對應(yīng)的適應(yīng)度值(BIC)和殘差平方和(Q)都比傳統(tǒng)的差分進化算法更小。這說明我們的方法在處理B樣條曲線最小二乘擬合時，擬合效果更好、擬合精度更高。

表1　傳統(tǒng)差分進化算法與本文算法的比較

表2　傳統(tǒng)差分進化算法與本文算法的比較

4.2　曲面擬合

將本文的算法延伸到曲面擬合問題上，以例4、例5為測試函數(shù)選取雙三次B樣條曲面進行擬合數(shù)值試驗。試驗均將擬合區(qū)域0≤x≤1、0≤y≤1分別32等分，進行數(shù)據(jù)點的直接采樣(無噪聲)；算法參數(shù)均設(shè)置為NP=50、Tmax=200。取平均平方誤差MSE為適應(yīng)度函數(shù)，定義如下：

式中:S(x,y)為擬合曲面，f(x,y)為采樣函數(shù)。

式中:x,y∈[0,1]。該函數(shù)為帶尖點的曲面。

圖4-圖7分別為根據(jù)例4和例5函數(shù)解析式直接用Matlab繪制的圖形和用本文算法繪制擬合曲面圖形的對比圖。由圖可知，使用本文算法得到的雙三次B樣條擬合曲面幾乎能準確構(gòu)造出采樣數(shù)據(jù)的原始圖像。多次試驗的數(shù)據(jù)顯示，對于測試函數(shù)4，大約迭代到第10～20次，平均平方誤差就已達到10-6，最終收斂到10-8，最佳內(nèi)節(jié)點個數(shù)為8×8。對于測試函數(shù)5，迭代到第10代左右，平均平方誤差達到10-9，最終收斂到10-12，最佳內(nèi)節(jié)點個數(shù)為5×5。這說明本文的算法收斂速度快，擬合精度高。

圖4　例4的根據(jù)函數(shù)解析式繪制的原圖

圖5　例4的本文算法的擬合效果圖

圖6　例5的根據(jù)函數(shù)解析式繪制的原圖

圖7　例5的本文算法的擬合效果圖

5　結(jié)　語

本文基于差分進化算法在處理數(shù)值優(yōu)化問題時所表現(xiàn)出的收斂速度快、求解精度高的優(yōu)勢，將差分進化算法進行改進，應(yīng)用于B樣條曲線曲面的最小二乘擬合問題。該方法選擇帶有隨機游走的交叉操作來克服DE/current-to-best易早熟收斂的缺點；并且設(shè)計出一種帶有混沌系統(tǒng)的局部搜索操作來跳出局部最優(yōu)解。將該方法與傳統(tǒng)的差分進化算法進行比較，試驗數(shù)據(jù)以及效果圖顯示，該方法在處理帶間斷和尖點的數(shù)據(jù)時能夠自適應(yīng)地確定節(jié)點的位置，并且與傳統(tǒng)的差分進化算法相比在擬合誤差精度上也有顯著的提高。其不足之處有兩點：首先，本文方法并沒有完全自適應(yīng)地同時確定節(jié)點的數(shù)量與位置，這主要是受限于差分進化算法的變異操作不允許我們將節(jié)點矢量設(shè)置為變長度的個體。其次，由于B樣條曲線和曲面不能精確地表示除拋物線外的二次曲線圓弧、和除拋物線面外的二次曲面，出于工業(yè)實際(特別是飛機外形的設(shè)計)精確表示二次曲線弧和二次曲面的需要，因此需要我們進一步研究有理B樣條的曲線曲面擬合問題。鑒于此，作者將進一步研究差分進化算法的改進，解決B樣條曲線曲面最小二乘擬合問題，實現(xiàn)節(jié)點數(shù)量和位置對不同實例的自適應(yīng)調(diào)整，并研究有理B樣條的曲線曲面擬合問題。

[1] 施法中. 計算機輔助幾何設(shè)計與非均勻有理B樣條[M].北京:高等教育出版社, 2013:301-336.

[2] Piegl L, Tiller W. The NURBS Book[M]. Berlin: Springer-Verlag, 1995.

[3] Yoshimoto F, Harada T, Yoshimoto Y. Data fitting with a spline using a real-coded genetic algorithm[J]. Computer-Aided Design, 2003, 35(8):751-760.

[4] 周明華，汪國昭. 基于遺傳算法的B樣條曲線和Bezier曲線的最小二乘擬合[J]. 計算機研究與發(fā)展, 2005, 42(1):135-143.

[5] 穆國旺, 臧婷, 趙罡. 用改進遺傳算法確定B樣條曲線的節(jié)點矢量 [J]. 計算機工程與應(yīng)用, 2006,42(1):88-90.

[6] Adi D I S, Shamsuddin S M, Ali A. Particle Swarm Optimization for NURBS Curve Fitting[C]// Sixth International Conference on Computer Graphics, Imaging and Visualization. IEEE Computer Society, 2009:259-263.

[7] Gálvez A, Iglesias A. Efficient particle swarm optimization approach for datafitting with free knot B-splines[J]. Computer-Aided Design, 2011, 43(12):1683-1692.

[8] ülker E, Arslan A. Automatic knot adjustment using an artificial immune system for B-spline curve approximation [J]. Information Sciences, 2009, 179:1483-1494.

[9] Zhao X, Zhang C, Yang B, et al. Adaptive knot placement using a GMM-based continuous optimization algorithm in B-spline curve approximation [J]. Computer-Aided Design, 2011, 43:598-604.

[10] 何兵朋, 馮仁忠, 余勝蛟. 基于差分進化算法的B樣條曲線曲面擬合[J]. 圖學(xué)學(xué)報, 2016, 37(2):178-183.

[11] Trejo-Caballero G, Rostro H, Garcia-Capulin C H, et al. Automatic Curve Fitting Based on Radial Basis Functions and a Hierarchical Genetic Algorithm[J]. Mathematical Problems in Engineering, 2015, 2015(731207).

[12] Han X, Quan L, Xiong X, et al. Diversity enhanced and local search accelerated gravitational search algorithm for data fitting with B-splines[J]. Engineering with Computers, 2015, 31(2):215-236.

[13] Vesterstrom J, Thomsen R. A comparative study of differential evolution, particle swarm optimization, and evolutionary algorithms on numerical benchmark problems[C]// Evolutionary Computation, 2004. CEC2004. Congress on. IEEE, 2004:1980-1987.

[14] Storn R, Price K. Differential evolution: a simple and efficient heuristic for global optimization over continuous spaces[J]. Kluwer Academic Publishers, 1997, 11(4):341-359.

[15] Storn R, Price K. Differential evolution: a fast and simple numerical optimizer[C]//1996 Biennial Conference of the North American Fuzzy Information Processing Society. New York, 1996:524-527.

[16] Tanabe R, Fukunaga A. Success-history based parameter adaptation for differential evolution[C]// IEEE Congress on Evolutionary Computation, 2013:71-78.

[17] Das S, Abraham A, Chakraborty U K, et al. Differential evolution using a neighborhood-based mutation operator[J]. IEEE Transactions on Evolutionary Computation, 2009, 13(3):526-553.

[18] Cui L, Li G, Lin Q, et al. Adaptive Differential Evolution Algorithm with Novel Mutation Strategies in Multiple Sub-populations[J]. Computers & Operations Research, 2015, 67:155-173.

[19] Zhan Z H, Zhang J. Enhance Differential Evolution with Random Walk [C]// Conference Companion on Genetic & Evolutionary Computation, 2012:1513-1514.

[20] Li G H, Lin Q Z, Cui L Z, et al. A novel hybrid differential evolution algorithm with modified CoDE and JADE[J]. Applied Soft Computing, 2016, 47:577-599.

[21] Liu G, Xiong C Q, Guo Z L. Enhanced differential evolution using random-based sampling and neighborhood mutation[J]. Soft Computing, 2015, 19:2173-2192.

[22] Cai Y, Zhao M, Liao J. Neighborhood guided differential evolution[J]. Soft Computing, 2016:1-44.

[23] Schwarz G E. Estimating the dimension of a model[J]. Annals of Statistics, 1978, 6(2):461-464.