曲率单调Bézier曲线G1插值算法及应用

何 川，赵 罡，王 伟，王爱增

(1.北京航空航天大学机械工程及自动化学院，北京 100191；2.北京航空航天大学虚拟现实技术与系统国家重点实验室，北京 100191)

在道路设计、路径规划等问题中，需要避免车辆行驶过程中的突然转向，因此设计师致力于构造曲率连续变化的缓和曲线。欧拉螺旋线(Euler spiral)具有曲率随弧长线性变化的特点，被广泛应用于路径规划[1-3]、曲线过渡[4-5]、形状修复等问题中[6-8]，然而，欧拉螺旋线采用超越方程的形式定义，求解困难，因此出现了许多利用多项式曲线或有理多项式曲线来近似欧拉螺旋线的方法[9-11]。

除了欧拉螺旋线，多项式曲线也常常用来设计缓和曲线。文献[12-14]采用了平面三次Bézier 螺旋线构造过渡曲线，并对文献[12,14]内容进行了推广，增加了Bézier 螺旋线的自由度[15-16]。

由于螺旋线的曲率随弧长线性变化，这一约束使得螺旋线的构造十分困难，工程实际中如果能够保证曲率单调递增或递减，则认为得到了符合光顺性要求的美学曲线。2006 年，FARIN[17]利用Class A矩阵表示Bézier 曲线控制边之间的几何关系，提出了一种构造曲率单调变化(monotone curvature variation，MCV) Bézier 曲线的方法，不过该理论被CAO 和WANG[18]证明不够完善，WANG 和ZHAO[19]给出了反例，并基于几何方法证明了一类平面任意次Bézier 曲线曲率单调的一个充分必要准则[20]。针对欧拉螺旋线G1 插值时的效率与精度问题，本文提出了一种可满足G1 插值条件的MCV Bézier 曲线构造方法，并与基于欧拉螺旋线的插值方法进行了对比，相关的数值实例证明了本文方法的有效性与实用性；且相比于文献[21]采用曲率单调的B 样条曲线实现G1 插值，本文构造方法简单直观，曲线弯曲度更好，算法的适用性与通用性更强。本文算法与图例均在Windows 10 系统下由MATLAB 语言实现。

1 欧拉螺旋线与MCV Bézier 曲线

1.1 欧拉螺旋线

欧拉螺旋线，又称为羊角螺旋线(cornu spiral)或回旋线(clothoid)，其曲率κ可表示为弧长L的线性函数，即

其中，0κ为L=0 时的初始曲率值；σ为欧拉螺旋线的曲率变化率。欧拉螺旋线的切线方向或旋转角θ关于弧长的表达式为[6]

对于规范化的欧拉螺旋线(图1)，其参数形式为

图1 规范化的欧拉螺旋线，(x0,y0)=(0,0)，L=4 Fig.1 Normalized Euler spiral,(x0,y0)=(0,0),L=4

其中，C(L)与S(L)为菲涅尔积分(Fresnel integral)[22]，其表达式为

1.2 MCV Bézier 曲线

定理1.平面k次Bézier 曲线定义式为

且当比例因子s≥1 时，曲线P(t)具有单调递减的曲率；当0＜s＜1 时，曲率单调递增。图2 显示了s≥1 时的一条曲率单减的三次Bézier 曲线及其对应的曲率梳。

图2 曲率单调递减的三次Bézier 曲线及其曲率梳，s=1.5，θ＝，V0=(20,0) ((a)三次MCV Bézier 曲线；(b)曲线的曲率梳)Fig.2 Cubic MCV Bézier curve with decreasing curvature and corresponding curvature comb,s=1.5,θ＝,V0=(20,0) ((a) Cubic MCV Bézier curve;(b) Corresponding curvature comb)

2 MCV Bézier 曲线G1 插值

定理 2.给定 2 个插值端点PA(xA,yA)和PB(xB,yB)，及其相应的方向角θA和θB，存在符合边界约束的k次MCV Bézier 曲线的充分条件为

具有满足s…co sθ≥ 1(s≥ 1)或s≤cosθ(0＜s＜1)的实数解，其中，θ=(θB-θA)/(k-1)，方向角θA和θB以实际绕过的圈数计算，可大于2π。和s是未知变量。

证明.对于次数k一定时，θ=(θB-θA)/(k-1)即为构造MCV Bézier 曲线的旋转角度，若方程组存在满足s…cosθ≥ 1(s≥ 1)或s≤cosθ(0＜s＜1)的实数解，即满足定理1 要求的MCV 条件，则由此构造的Bézier 曲线曲率必单调。证毕.

下面分析边界约束条件对MCV Bézier 曲线解的存在性的影响。

由于Bézier 曲线的仿射不变性，故可将其中一个端点PA固定在原点，其方向角θA作为X轴正方向，并以此建立直角坐标系，则PA=(0,0)，θA=0，记2 端点连线与X轴正方向的夹角为θʹ (图3)，为讨论方便，本文仅讨论PB位于第一象限的情形，即xB＞0，yB＞0。

图3 θB＜θʹ，无解 Fig.3 θB＜θʹ,no solution exists

2.1 次数k=2

图4 θ B＞，无解 Fig.4 θ B＞,no solution exists

图5 二次MCV Bézier 曲线控制边Fig.5 The control edge of quadratic MCV Bézier curve

图6 约束条件：PB=(20,20)，θA=0 Fig.6 Constraints:PB=(20,20),θA=0

图7 二次MCV Bézier 曲线，s=2.732 051，||V0||=8.452 995，θ＝θB＝ ((a)二次MCV Bézier 曲线；(b)曲线曲率图)Fig.7 Quadratic MCV Bézier curve,s=2.732 051,||V0||=8.452 995,θ＝θB＝ ((a) Quadratic MCV Bézier curve;(b) Corresponding curvature plot)

图8 二次MCV Bézier 曲线，s=0.565 685，=50，θ＝θB＝ ((a)二次MCV Bézier 曲线；(b)曲线曲率图)

Fig.8Quadratic MCV Bézier curve,s=0.565 685,=50,θ＝θB＝((a) Quadratic MCV Bézier curve; (b) Corresponding curvature plot)

2.2 次数k 可变时

命题2.对于任意次数k，若θB＜θʹ (图9)，则不存在满足约束的MCV Bézier 曲线。

图9 θB＜θʹ，无MCV Bézier 曲线解 Fig.9 θB＜θʹ,without solution of MCV Bézier curve

证明.θB＜θʹ的情形如图9 所示，由k次MCV Bézier 曲线的构造过程可知，控制边Δbi(i=1,2,…,k-1)的方向角θi=(i-1)…θ是单调递增的序列，且倒数第二个控制点bk-1位于虚线QPB上第一象限部分，其方向角θk-1=θB。

假设下标为i的控制点bi位于虚线PAPB下侧(即三角形ΔPAPPB内)，下标为i+1 的控制点bi+1位于虚线PAPB上侧，则控制边Δbi=bi+1-bi将由下至上穿过PAPB，其方向角满足θi＞θʹ＞θB=θk-1，与已知控制边的方向角的递增序列矛盾，故θB＜θʹ时不存在满足约束的MCV Bézier 曲线证毕.

命题2 说明MCV Bézier 曲线不存在S 型曲线。

命题3.当θB＞θʹ时，总能找到合适的次数k(k≥2)，使得k次MCV Bézier 曲线满足约束。

证明.单调递增的序列，当次数k为自由变量时，MCV Bézier 插值曲线的方程组为

推论.当θB＞θʹ时，若k=m≥2 时存在MCV Bézier 曲线满足约束，则当k≥m+1 时，也存在MCV Bézier 曲线满足约束，即MCV Bézier 曲线存在多解，且解的极限为0θ→，1s→ 。

2.3 算法

基于G1 边界约束的k次MCV Bézier 曲线插值算法的步骤如下：

3 数值实例

图10 G1 插值MCV Bézier 曲线，k=6 ((a)本文算法求解的MCV Bézier 曲线;(b)曲线曲率图) Fig.10 G1 interpolation MCV Bézier curve,k=6 ((a) MCV Bézier curve obtained by our algorithm; (b) Corresponding curvature plot)

4 与基于欧拉螺旋线的G1 插值算法比较

图11 G1 插值MCV Bézier 曲线，k=18 ((a)本文算法 求解的MCV Bézier 曲线；(b)曲线曲率图) Fig.11 G1 interpolation MCV Bézier curve,k=18 ((a) MCV Bézier curve obtained by our algorithm; (b) Corresponding curvature plot)

图12 G1 插值MCV Bézier 曲线，k=26 ((a)本文算法求解的MCV Bézier 曲线；(b)曲线曲率图) Fig.12 G1 interpolation MCV Bézier curve,k=26 ((a) MCV Bézier curve obtained by our algorithm; (b) Corresponding curvature plot)

欧拉螺旋线的曲率随弧长线性变化，是光顺 性非常好的一类曲线，然而其复杂的数学形式只能用数值方法进行近似求解，因此带来了求解效率与计算精度等问题。而本文方法采用的Bézier曲线，属于简单的多项式函数，计算结果精确，效率高。与传统基于欧拉螺旋线的G1 插值算法相比，优点如下：

(1) 构造简单，计算高效；

(2) 计算精确，插值结果无误差；

(3) 与现有的NURBS方法兼容。

图13～16 展示了相同约束条件下采用欧拉螺旋线插值与本文MCV Bézier 曲线插值的结果。

给定约束PA=(0,0)，θA=0，PB=(50,50)，欧拉螺旋线与MCV Bézier 曲线求解结果分别如图13和图14 所示。欧拉螺旋线初始曲率0κ=0.007406，弧长L=85.919605，曲率变化率σ=0.000395。本文算法MCV Bézier 曲线插值结果如图14 所示，k=5，s=0.797364，表1 为2 种方法的误差与时间对比结果。

表1 误差与时间对比1Table 1 First comparison of position error and time

图13 欧拉螺旋线插值第一例((a)欧拉螺旋线；(b)欧拉螺旋线曲率图(斜率放大20 倍))Fig.13 First interpolation example of Euler spiral ((a) Euler spiral;(b) Corresponding curvature plot (The slope is scaled by 20 times))

图14 本文算法MCV Bézier 曲线插值第一例((a) MCV Bézier 曲线；(b)曲线曲率图) Fig.14 First interpolation example of MCV Bézier curve based on our algorithm ((a) MCV Bézier curve; (b) Corresponding curvature plot)

给定约束PA=(0,0)，θA=0，PB=(40,30)，θB=π-arctan欧拉螺旋线与MCV Bézier 曲线求解结果分别如图15 和图16 所示。欧拉螺旋线初始曲率κ0=-0.010878，弧长L=67.290185，曲率变化率σ=0.001427。本文算法MCV Bézier 曲线插值结果如图16 所示，k=4，s=0.5876780，对比结果见表2。

图15 欧拉螺旋线插值第二例((a)欧拉螺旋线；(b)欧拉螺旋线曲率图(斜率放大20 倍)) Fig.15 Second interpolation example of Euler spiral ((a) Euler spiral;(b) Corresponding curvature plot (The slope is scaled by 20 times))

图16 本文算法MCV Bézier 曲线插值第二例 ((a) MCV Bézier 曲线；(b)曲线曲率图) Fig.16 Second interpolation example of MCV Bézier curve based on our algorithm ((a) MCV Bézier curve; (b) Corresponding curvature plot)

表2 误差与时间对比2Table 2 Second comparison of position error and time

5 分析与讨论

5.1 算法推广

本文给出了基于第一象限2 端点的G1 约束MCV Bézier 样条曲线插值算法，该算法能够推广至更加复杂的端点约束情形，如非第一象限的端点插值及多数据点插值。

针对尾点位于非第一象限的情况，可以通过补充中间插值点的方式进行分段处理，使得每一段的尾点位于局部坐标系的第一象限，将原问题转化为多个第一象限的G1 插值问题。

对于给定一个平面数据点序列的插值问题，同样需要采取分段处理的方式，除首尾2 端点外，中间数据点没有切线方向的限制。由于减少了每一段插值曲线尾点的切向约束(最后一段除外)，因此相邻2 点的插值自由度增加，曲线设计更加灵活。

分段处理引入的一个问题是，分段点处曲率可能会发生突变，且多段拼接而成的整条曲线，只能保证分段曲线的曲率单调性。要解决以上问题，需要引入G2 插值，这是将来的研究方向。

给定平面数据点序列PA=P0=(0,0)，P1=(14,3)，P2=(18,16)，P3=(12,22)，PB=P4=(-2,26) (图17，粉色数据点)，以及首尾切线方向θA=0，θB=π。一种分段插值结果如图17 蓝色曲线所示，该曲线是由4段二次MCV Bézier 曲线组合而成，每个中间数据点即为分段连接点。组合曲线的曲率图分段单调(图18)，在中间3 个数据点处曲率发生突变。

图17 给定平面数据点列的MCV Bézier 曲线插值 Fig.17 MCV Bézier curve interpolation of given planar data point

图18 组合曲线的曲率图Fig.18 Curvature plot of the whole composite curve

5.2 与基于平面三次B 样条曲线G1 插值算法比较

文献[21]提出了一种曲率单调的平面三次B 样条曲线构造算法，且将其应用于形状修复，属于一类G1 插值问题。该算法是通过前2 条控制边矢量的夹角及模长，确定后一条控制边矢量的范围以保证曲率的单调性。但未具体给出如何选取控制边矢量的方向及大小，算法随机性较强；且相邻边矢量的夹角逐渐减小，随着控制点的增多，曲线趋于直线，曲线的角度变化范围有限。

对于5.1 节所示的复杂平面点列，文献[21]算法受限于首尾切线方向的角度变化，无法生成通过各点的插值曲线，相关理论有待进一步研究。而本文方法基于变换矩阵，构造更加快速、直观，曲线整体的弯曲度更大，能够处理的切向约束范围更广，因此算法的适用性与通用性更强。

6 结束语

本文提出了一种基于G1 约束的MCV Bézier 样条曲线插值算法。首先阐述了给定边界约束下MCV Bézier 曲线解的存在性条件并进行证明；接着，理论与实践表明，该算法可以灵活构造满足G1 插值条件的曲率单调Bézier 曲线；对比基于欧拉螺旋线的插值算法，本文算法更加精确且高效；作为处理第一象限2 端点G1 约束问题的基础，该算法可以推广至更加复杂的约束问题求解，相比于曲率单调的平面三次B 样条曲线，MCV Bézier 曲线的构造更加快速、直观，算法的适用性与通用性更强。

未来工作是寻求当PB位置与方向角θB任意取值时，构造满足约束的MCV Bézier 样条曲线的有效解法，以解决复杂情形下的路径规划等工程问题；分段处理将引入曲率的突变，G2 插值问题相比G1 更加重要且复杂，这是未来的研究重点。