陳 娟，何 歆，張起航，孫義環(huán)，陸利正

(浙江工商大學(xué)統(tǒng)計(jì)與數(shù)學(xué)學(xué)院，浙江 杭州 310018)

光順曲線具有簡單的曲率輪廓(例如螺線的曲率是單調(diào)的),它的構(gòu)造方法是計(jì)算機(jī)輔助設(shè)計(jì)和相關(guān)領(lǐng)域的核心問題[1-2].在路徑設(shè)計(jì)和軌道規(guī)劃時(shí),使用光順曲線以滿足幾何造型的需求就顯得尤為重要[3].通常來說,曲線光順由某些內(nèi)在的能量函數(shù)來衡量;除了滿足一定端點(diǎn)條件外,曲線表示的其余自由度往往通過優(yōu)化某些能量函數(shù)來確定[4].在曲線擬合時(shí),常見做法也是采用能量作為目標(biāo)函數(shù)的構(gòu)成部分,以改善擬合曲線的光順程度.

光順曲線的構(gòu)造方法大致上可分為兩大類.第一類方法使用螺線分段擬合自由曲線,因此具有分段單調(diào)的曲率分布.Walton等[5]借助歐拉螺線提出G1插值方法.但更多研究側(cè)重于使用多項(xiàng)式曲線構(gòu)造G2插值的螺線[6-9],以適應(yīng)計(jì)算機(jī)輔助設(shè)計(jì)(CAD)系統(tǒng)的要求.然而,這類方法須滿足一定端點(diǎn)條件才能保證得到的插值曲線是螺線,這意味著在某些端點(diǎn)條件下不存在G2插值的螺線,因此不適用于任意的端點(diǎn)數(shù)據(jù),從而給外形設(shè)計(jì)造成不便.

然而,已有研究工作主要集中在平面曲線情形,很少專門針對空間曲線的系統(tǒng)研究.部分G1插值構(gòu)造方法[10-11,13-14]可直接推廣用于空間曲線的G1插值問題.由于空間曲線必須考慮在插值點(diǎn)處的曲率和Frenet標(biāo)架,已有平面曲線G2插值構(gòu)造方法[6-9]還不能應(yīng)用于空間曲線的G2插值問題.

本文提出基于jerk能量最小化的空間五次G2插值曲線構(gòu)造方法.對于任意給定的G2數(shù)據(jù),提出端點(diǎn)滿足特定曲率和Frenet標(biāo)架且含4個(gè)自由參數(shù)的空間五次Bézier曲線模型;然后,借助能量優(yōu)化確定這4個(gè)自由參數(shù),以生成光順曲線.目前這方面研究工作仍不多,Lu[12]的工作由于采用最小化應(yīng)變能優(yōu)化曲線,導(dǎo)致目標(biāo)曲線的曲率輪廓在許多情況下不夠理想,并且只適用于平面曲線的插值問題.本文主要的目的是在Lu[12]的工作基礎(chǔ)上做出改進(jìn),使其更好地滿足實(shí)際需要.

借鑒文獻(xiàn)[14-17]的做法,本文使用jerk能量近似表示曲線的曲率變化.首先,定義目標(biāo)函數(shù)為jerk能量和曲線長度的加權(quán)組合,將曲線長度作為修正項(xiàng)以抑制曲線長度；然后,化簡后的目標(biāo)函數(shù)是二元四次多項(xiàng)式,并且很容易求得它的梯度和Hessian矩陣；最后,使用投影牛頓法確定參數(shù)的最優(yōu)取值.相比于Lu的工作[12],本文方法生成曲線的曲率變化更加簡單平緩,因此更符合光順要求.另外,本文方法也適用于空間曲線情形,滿足幾何造型的實(shí)際需求.

1 空間五次G2插值曲線模型

假設(shè)兩端點(diǎn)處的G2數(shù)據(jù)為{P0,T0,κ0N0;P1,T1,κ1N1}，其中κi,Ti和Ni分別表示曲線在端點(diǎn)Pi的曲率、單位切向量和單位法向量.那么曲線在端點(diǎn)處的Frenet標(biāo)架分別表示為{T0,N0,B0}和{T1,N1,B1},Bi=Ti×Ni.

空間五次G2插值曲線的構(gòu)造問題是指針對給定的G2數(shù)據(jù),構(gòu)造一條滿足插值條件的空間五次Bézier曲線

(1)

圖1 空間五次G2插值方法Fig.1 Spatial quintic G2 interpolating method

假設(shè)五次Bézier曲線b(t)的長度為L以及弧長參數(shù)表示為b(s),s=s(t)∈[0,L]，那么曲線b(t)的導(dǎo)數(shù)計(jì)算為：

b′(t)=s′(t)b′(s),

b″(t)=s″(t)b′(s)+(s′(t))2b″(s).

由于曲線在兩端點(diǎn)處須滿足給定的曲率和Frenet標(biāo)架,再借助Frenet公式[19],所以：

其中自由參數(shù)αi=s′(i),βi=s″(i),i=0,1.再根據(jù)Bézier曲線的性質(zhì)[1],控制頂點(diǎn)計(jì)算為：

(2)

(3)

(4)

(5)

只要確定αi和βi的參數(shù)值,就可代入計(jì)算得到五次Bézier曲線的控制頂點(diǎn).對于任意給定的G2數(shù)據(jù),該五次模型不僅能生成相應(yīng)的插值曲線,而且為后續(xù)形狀優(yōu)化提供4個(gè)可修改參數(shù).同時(shí),該五次模型只取決于相鄰兩個(gè)插值點(diǎn)的端點(diǎn)條件,適合自由曲線的分段構(gòu)造.

2 jerk能量和曲線優(yōu)化構(gòu)造

對于曲線b(t),t∈[0,1]， Farin[13]提出曲率變化e是較好衡量光順的觀點(diǎn),其定義為

(6)

雖然b(t)是五次多項(xiàng)式曲線,但κ(t)是非線性的,[κ′(t)]2是次數(shù)很高的有理分式,這不利于后續(xù)優(yōu)化求解.尤其當(dāng)碰到自由曲線由許多段構(gòu)成時(shí),計(jì)算量就非常大且很難做到實(shí)時(shí)求解.因此,使用jerk能量

(7)

近似表示曲率變化.當(dāng)表示曲線的參數(shù)是弧長參數(shù)時(shí),該曲線的jerk能量和曲率變化剛好相等；對于一般的多項(xiàng)式曲線而言,jerk能量是曲率變化一個(gè)較好的近似.

另一方面,外形設(shè)計(jì)時(shí)曲線長度也是需要考慮的重要因素.對于一系列插值點(diǎn)來說,構(gòu)造的插值曲線應(yīng)盡量接近插值點(diǎn)表示的曲線輪廓,同時(shí)避免過多振蕩.這要求兩個(gè)插值點(diǎn)之間的曲線形狀不能變化太大,因此在構(gòu)造光順曲線時(shí)就要求曲線長度也保持適當(dāng)小.基于上述考慮,構(gòu)造目標(biāo)函數(shù)

(8)

通常取λ為較小的正數(shù),因此稱第2項(xiàng)為修正項(xiàng),用于抑制目標(biāo)曲線的長度.在目標(biāo)函數(shù)中加入修正項(xiàng)有兩個(gè)好處:第一,可進(jìn)一步增強(qiáng)目標(biāo)函數(shù)的凸性,提高優(yōu)化求解時(shí)的穩(wěn)定性;第二,通過抑制‖b′(t)‖的整體變化使曲線參數(shù)更接近弧長參數(shù),從而使jerk能量更接近曲率變化.

利用Bézier曲線的導(dǎo)數(shù)公式[1]

其中Δbi=bi+1-bi和Δ3bi=bi+3-3bi+2+3bi+1-bi是控制頂點(diǎn)的差分算子,以及Bernstein多項(xiàng)式的積分公式[20]

式(8)展開為關(guān)于控制頂點(diǎn)的二次函數(shù):

(9)

其中

[gij]=

[hij]=

將控制頂點(diǎn)的表達(dá)式(2)～(5)代入目標(biāo)函數(shù)(9)后,f本質(zhì)上是α0,α1的四次函數(shù)和β0,β1的二次函數(shù).顯然,曲線b(t)在兩端點(diǎn)處切向量的長度分別等于α0和α1.為使曲線盡量光順,α0和α1的取值不能太大也不能太小.因此,限定(α0,α1)的可行域?yàn)?/p>

D={(α0,α1):d0≤α0≤d1,d2≤α1≤d3},

其中di>0是給定的界.最后通過求解如下約束優(yōu)化問題

(10)

確定4個(gè)未知數(shù)α0,α1,β0,β1， 從而最終得到五次Bézier曲線的控制頂點(diǎn).

3 算法實(shí)現(xiàn)

目標(biāo)函數(shù)(9)的梯度計(jì)算如下.首先,

根據(jù)式(2)～(5),bi的偏導(dǎo)數(shù)中只有6個(gè)不為零,羅列如下:

?b1/?α0=(1/5)T0,

?b2/?α0=(2/5)T0+(κ0α0/10)N0,

?b4/?α1=(-1/5)T1,

?b3/?α1=(-2/5)T1+(κ1α1/10)N1,

?b2/?β0=(1/20)T0,

?b3/?β1=(1/20)T1.

應(yīng)用鏈?zhǔn)椒▌t,得

令?f/?β0=?f/?β1=0， 得到關(guān)于β0,β1的線性方程組(其中cij是依賴于兩端點(diǎn)處G2數(shù)據(jù)的常數(shù)):

從該方程組可知,β0,β1是α0,α1的二次函數(shù),即β0=β0(α0,α1),β1=β1(α0,α1).這意味著當(dāng)f取到極值時(shí),β0,β1總是二次依賴于α0,α1.

最后,將β0=β0(α0,α1)和β1=β1(α0,α1)代入目標(biāo)函數(shù)(9)后,函數(shù)f是α0,α1的二元四次多項(xiàng)式.于是,問題(10)最終化簡為

(11)

使用投影牛頓法[21]對上述問題進(jìn)行數(shù)值求解,通常取初值為(α0,α1)=(1,1).函數(shù)f關(guān)于α0,α1的梯度和Hessian矩陣根據(jù)鏈?zhǔn)椒▌t重新計(jì)算;雖然具體公式表達(dá)式有些復(fù)雜,但借助數(shù)學(xué)軟件容易計(jì)算.

4 實(shí) 例

4.1 平面曲線

圖中圓點(diǎn)為控制頂點(diǎn)，折線表示控制多邊形，下同.圖2 本文方法(虛線)與Lu[12]方法(實(shí)線)的比較Fig.2 Comparison of the method of this article(dashed) and Lu’s[12](solid)

下面給出一些平面曲線例子來驗(yàn)證本文方法的可行性,并與文獻(xiàn)[12]進(jìn)行比較.在所有例子中,取權(quán)重λ=0.01和可行域D=[0.1,5]2.

考慮6個(gè)代表性的G2插值問題,不失一般性,令插值點(diǎn)為P0=(0,0),P1=(1,0)， 以及相應(yīng)的單位切向量為T0=(cosφ0,-sinφ0),T1=(cosφ1,sinφ1).圖2所示為本文方法與Lu[12]的方法的對比結(jié)果:曲線圖列在左邊,而對應(yīng)的曲率圖列在右邊.對比表明,本文方法得到曲線具有曲率變化更加平緩的特點(diǎn),因此更符合光順要求.表1列出這兩種方法得到曲線的曲率變化和jerk能量.對比數(shù)據(jù)可以看出,本文方法得到曲線的曲率變化值更小.

表1 圖2中兩種方法得到曲線的曲率變化和jerk能量Tab.1 Curvature variation and jerk energy of thecurve obtained by the two methods in figure 2

歐拉螺線具有線性曲率，是一類廣泛使用的光順曲線,常用于路徑設(shè)計(jì).考慮一段歐拉螺線(圖3(a)中實(shí)線):

首先將該螺線分成3段,然后分別用文獻(xiàn)[6,12]中的方法和本文方法得到五次G2插值曲線.從圖3可知,本文方法得到的曲線不僅更接近歐拉螺線,而且具有更完美的曲率分布.

圖3 歐拉螺線的五次G2分段插值Fig.3 Quintic G2 piecewise interpolation of Euler spiral

4.2 空間曲線

本文方法適用于空間曲線的分段插值構(gòu)造.例如,對于Viviani曲線和Eightknot曲線[19],其定義分別為：

圖4 Viviani曲線的五次G2分段插值Fig.4 Quintic G2 piecewise interpolation of Viviani curve

圖5 Eightknot曲線的五次G2分段插值Fig.5 Quintic G2 piecewise interpolation of Eightknot curve

5 結(jié) 論

本文提出基于jerk能量最小化的空間五次G2插值曲線構(gòu)造方法.衡量曲線光順的標(biāo)準(zhǔn)有幾個(gè),而曲率變化是較好的標(biāo)準(zhǔn).用jerk能量近似表示曲率變化以降低約束優(yōu)化求解時(shí)的難度,因此可滿足大規(guī)模插值問題快速求解的需要.許多對比實(shí)例表明本文方法能夠生成更光順的曲線,也可用于解決一系列空間數(shù)據(jù)點(diǎn)的插值曲線構(gòu)造問題.下一步研究重點(diǎn)是開展在路徑設(shè)計(jì)等方面的應(yīng)用[3].

由于目標(biāo)函數(shù)是jerk能量和曲線長度的加權(quán)組合,所以權(quán)重λ的取值對最終結(jié)果有一定影響.如果λ越小,那么生成曲線的jerk能量就越小.對于大多數(shù)例子建議取默認(rèn)值λ=0.01，而少部分例子需要交互修改.

此外,本文關(guān)于空間曲線高質(zhì)量插值構(gòu)造方法可用于曲面的插值與擬合問題.未來將進(jìn)一步開展曲面問題的研究.