陈 娟，何 歆，张起航，孙义环，陆利正

(浙江工商大学统计与数学学院，浙江 杭州 310018)

光顺曲线具有简单的曲率轮廓(例如螺线的曲率是单调的),它的构造方法是计算机辅助设计和相关领域的核心问题[1-2].在路径设计和轨道规划时,使用光顺曲线以满足几何造型的需求就显得尤为重要[3].通常来说,曲线光顺由某些内在的能量函数来衡量;除了满足一定端点条件外,曲线表示的其余自由度往往通过优化某些能量函数来确定[4].在曲线拟合时,常见做法也是采用能量作为目标函数的构成部分,以改善拟合曲线的光顺程度.

光顺曲线的构造方法大致上可分为两大类.第一类方法使用螺线分段拟合自由曲线,因此具有分段单调的曲率分布.Walton等[5]借助欧拉螺线提出G1插值方法.但更多研究侧重于使用多项式曲线构造G2插值的螺线[6-9],以适应计算机辅助设计(CAD)系统的要求.然而,这类方法须满足一定端点条件才能保证得到的插值曲线是螺线,这意味着在某些端点条件下不存在G2插值的螺线,因此不适用于任意的端点数据,从而给外形设计造成不便.

然而,已有研究工作主要集中在平面曲线情形,很少专门针对空间曲线的系统研究.部分G1插值构造方法[10-11,13-14]可直接推广用于空间曲线的G1插值问题.由于空间曲线必须考虑在插值点处的曲率和Frenet标架,已有平面曲线G2插值构造方法[6-9]还不能应用于空间曲线的G2插值问题.

本文提出基于jerk能量最小化的空间五次G2插值曲线构造方法.对于任意给定的G2数据,提出端点满足特定曲率和Frenet标架且含4个自由参数的空间五次Bézier曲线模型;然后,借助能量优化确定这4个自由参数,以生成光顺曲线.目前这方面研究工作仍不多,Lu[12]的工作由于采用最小化应变能优化曲线,导致目标曲线的曲率轮廓在许多情况下不够理想,并且只适用于平面曲线的插值问题.本文主要的目的是在Lu[12]的工作基础上做出改进,使其更好地满足实际需要.

借鉴文献[14-17]的做法,本文使用jerk能量近似表示曲线的曲率变化.首先,定义目标函数为jerk能量和曲线长度的加权组合,将曲线长度作为修正项以抑制曲线长度；然后,化简后的目标函数是二元四次多项式,并且很容易求得它的梯度和Hessian矩阵；最后,使用投影牛顿法确定参数的最优取值.相比于Lu的工作[12],本文方法生成曲线的曲率变化更加简单平缓,因此更符合光顺要求.另外,本文方法也适用于空间曲线情形,满足几何造型的实际需求.

1 空间五次G2插值曲线模型

假设两端点处的G2数据为{P0,T0,κ0N0;P1,T1,κ1N1}，其中κi,Ti和Ni分别表示曲线在端点Pi的曲率、单位切向量和单位法向量.那么曲线在端点处的Frenet标架分别表示为{T0,N0,B0}和{T1,N1,B1},Bi=Ti×Ni.

空间五次G2插值曲线的构造问题是指针对给定的G2数据,构造一条满足插值条件的空间五次Bézier曲线

(1)

图1 空间五次G2插值方法Fig.1 Spatial quintic G2 interpolating method

假设五次Bézier曲线b(t)的长度为L以及弧长参数表示为b(s),s=s(t)∈[0,L]，那么曲线b(t)的导数计算为：

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

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

由于曲线在两端点处须满足给定的曲率和Frenet标架,再借助Frenet公式[19],所以：

其中自由参数αi=s′(i),βi=s″(i),i=0,1.再根据Bézier曲线的性质[1],控制顶点计算为：

(2)

(3)

(4)

(5)

只要确定αi和βi的参数值,就可代入计算得到五次Bézier曲线的控制顶点.对于任意给定的G2数据,该五次模型不仅能生成相应的插值曲线,而且为后续形状优化提供4个可修改参数.同时,该五次模型只取决于相邻两个插值点的端点条件,适合自由曲线的分段构造.

2 jerk能量和曲线优化构造

对于曲线b(t),t∈[0,1]， Farin[13]提出曲率变化e是较好衡量光顺的观点,其定义为

(6)

虽然b(t)是五次多项式曲线,但κ(t)是非线性的,[κ′(t)]2是次数很高的有理分式,这不利于后续优化求解.尤其当碰到自由曲线由许多段构成时,计算量就非常大且很难做到实时求解.因此,使用jerk能量

(7)

近似表示曲率变化.当表示曲线的参数是弧长参数时,该曲线的jerk能量和曲率变化刚好相等；对于一般的多项式曲线而言,jerk能量是曲率变化一个较好的近似.

另一方面,外形设计时曲线长度也是需要考虑的重要因素.对于一系列插值点来说,构造的插值曲线应尽量接近插值点表示的曲线轮廓,同时避免过多振荡.这要求两个插值点之间的曲线形状不能变化太大,因此在构造光顺曲线时就要求曲线长度也保持适当小.基于上述考虑,构造目标函数

(8)

通常取λ为较小的正数,因此称第2项为修正项,用于抑制目标曲线的长度.在目标函数中加入修正项有两个好处:第一,可进一步增强目标函数的凸性,提高优化求解时的稳定性;第二,通过抑制‖b′(t)‖的整体变化使曲线参数更接近弧长参数,从而使jerk能量更接近曲率变化.

利用Bézier曲线的导数公式[1]

其中Δbi=bi+1-bi和Δ3bi=bi+3-3bi+2+3bi+1-bi是控制顶点的差分算子,以及Bernstein多项式的积分公式[20]

式(8)展开为关于控制顶点的二次函数:

(9)

其中

[gij]=

[hij]=

将控制顶点的表达式(2)～(5)代入目标函数(9)后,f本质上是α0,α1的四次函数和β0,β1的二次函数.显然,曲线b(t)在两端点处切向量的长度分别等于α0和α1.为使曲线尽量光顺,α0和α1的取值不能太大也不能太小.因此,限定(α0,α1)的可行域为

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

其中di>0是给定的界.最后通过求解如下约束优化问题

(10)

确定4个未知数α0,α1,β0,β1， 从而最终得到五次Bézier曲线的控制顶点.

3 算法实现

目标函数(9)的梯度计算如下.首先,

根据式(2)～(5),bi的偏导数中只有6个不为零,罗列如下:

∂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.

应用链式法则,得

令∂f/∂β0=∂f/∂β1=0， 得到关于β0,β1的线性方程组(其中cij是依赖于两端点处G2数据的常数):

从该方程组可知,β0,β1是α0,α1的二次函数,即β0=β0(α0,α1),β1=β1(α0,α1).这意味着当f取到极值时,β0,β1总是二次依赖于α0,α1.

最后,将β0=β0(α0,α1)和β1=β1(α0,α1)代入目标函数(9)后,函数f是α0,α1的二元四次多项式.于是,问题(10)最终化简为

(11)

使用投影牛顿法[21]对上述问题进行数值求解,通常取初值为(α0,α1)=(1,1).函数f关于α0,α1的梯度和Hessian矩阵根据链式法则重新计算;虽然具体公式表达式有些复杂,但借助数学软件容易计算.

4 实 例

4.1 平面曲线

图中圆点为控制顶点，折线表示控制多边形，下同.图2 本文方法(虚线)与Lu[12]方法(实线)的比较Fig.2 Comparison of the method of this article(dashed) and Lu’s[12](solid)

下面给出一些平面曲线例子来验证本文方法的可行性,并与文献[12]进行比较.在所有例子中,取权重λ=0.01和可行域D=[0.1,5]2.

考虑6个代表性的G2插值问题,不失一般性,令插值点为P0=(0,0),P1=(1,0)， 以及相应的单位切向量为T0=(cosφ0,-sinφ0),T1=(cosφ1,sinφ1).图2所示为本文方法与Lu[12]的方法的对比结果:曲线图列在左边,而对应的曲率图列在右边.对比表明,本文方法得到曲线具有曲率变化更加平缓的特点,因此更符合光顺要求.表1列出这两种方法得到曲线的曲率变化和jerk能量.对比数据可以看出,本文方法得到曲线的曲率变化值更小.

表1 图2中两种方法得到曲线的曲率变化和jerk能量Tab.1 Curvature variation and jerk energy of thecurve obtained by the two methods in figure 2

欧拉螺线具有线性曲率，是一类广泛使用的光顺曲线,常用于路径设计.考虑一段欧拉螺线(图3(a)中实线):

首先将该螺线分成3段,然后分别用文献[6,12]中的方法和本文方法得到五次G2插值曲线.从图3可知,本文方法得到的曲线不仅更接近欧拉螺线,而且具有更完美的曲率分布.

图3 欧拉螺线的五次G2分段插值Fig.3 Quintic G2 piecewise interpolation of Euler spiral

4.2 空间曲线

本文方法适用于空间曲线的分段插值构造.例如,对于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 结 论

本文提出基于jerk能量最小化的空间五次G2插值曲线构造方法.衡量曲线光顺的标准有几个,而曲率变化是较好的标准.用jerk能量近似表示曲率变化以降低约束优化求解时的难度,因此可满足大规模插值问题快速求解的需要.许多对比实例表明本文方法能够生成更光顺的曲线,也可用于解决一系列空间数据点的插值曲线构造问题.下一步研究重点是开展在路径设计等方面的应用[3].

由于目标函数是jerk能量和曲线长度的加权组合,所以权重λ的取值对最终结果有一定影响.如果λ越小,那么生成曲线的jerk能量就越小.对于大多数例子建议取默认值λ=0.01，而少部分例子需要交互修改.

此外,本文关于空间曲线高质量插值构造方法可用于曲面的插值与拟合问题.未来将进一步开展曲面问题的研究.