连续小线段的高精度B样条曲线拟合方法

封雨鑫,钟 菲,余 强

（深圳市大族智能控制科技有限公司，广东 深圳 518101）

1 引言

随着智能制造行业的迅猛发展，现代数控系统正向着高速度、高精度、高可靠性等方面发展。如何在保证加工精度的同时最大可能的提高加工效率，是当前数据控系统研究的热点。传统的数控系统在加工曲线曲面轨迹时，为了满足精度需求，通常将它们离散成一系列首尾相连的微小线段来逼近原始加工轨迹，不可避免地产生了大量加工代码。这样不仅增加了数据存储、传输以及处理的负担，还会因为加工这些小线段速度和加速度的频繁变化，降低了零件表面加工质量和效果。现代数控加工中，通常会使用连续性更好的参数曲线，对小线段表示的路径拟合[1]，使得拟合后的曲线路径变长，同时保证了曲线除首末端点外均为二阶连续，即加速度连续，这样不仅可以提高机床加工速度，还可以减少机床加工过程中的振动，提高整体加工质量。

B样条曲线因其表达式的通用性和应用的广泛性，使得B 样条拟合技术成为了国内外许多专家和学者们研究的热点[2-4]。这些研究大多将输入数据作为参数曲线的控制点，未考虑拟合误差和曲线光顺的问题，难以直接应用于数控系统中大量小线段的加工。为了提高数控系统中直线段的加工效率，改善线段衔接处的连续性，研究人员常采用局部过渡(如：倒圆角[5]，PH曲线[6]，羊角曲线[7]等)方法或全局拟合方法对连续直线段进行平滑。而对于连续小线段路径而言，局部平滑过渡会显著降低加工效率。对此，一些研究学者根据B样条特性对连续小线段的拟合进行了研究[8-10]。Deng[8]等将迭代逼近和最小二乘法相结合，提出了一种渐进迭代逼近最小二乘法的B样条拟合方法，该算法通过建立型值点与B样条曲线对应点距离的平方和为目标函数，利用给定步长最速下降法来求解目标函数。但该方法未考虑弓高误差的约束，同时存在数值计算不稳定的问题，难以应用于实际加工中。为了提高样条拟合精度和数值计算的稳定性，Qing[9]等提出了一种带弓高误差约束迭代逼近B 样条拟合方法。该方法在样条拟合过程中通过渐进迭代逼近和局部细化的方法，使拟合曲线满足弓高误差约束，同时保证了样条曲线的凸包性。Du[10]等提出了一种基于优势点的误差有界b 样条曲线逼近方法。该方法根据连续小线段的几何特征和加工工艺，采用双边弓高误差选择合适的拟合点，然后根据最小二乘法迭代来迭代逼近原始路径。本文在前人的基础之上，提出了一种新的三次B 样条高精度拟合连续微小线段的方法。为了减少迭代和误差调整对加工小线段的效率，该方法首先采用连续小线段的转角限制和双边弓高误差限制，选择满足连续小线段可被拟合的型值点集；接着，采用凹凸性准则对可拟合点集分段后拟合成参数曲线；最后，利用插值拟合和误差分析模块，保证拟合的参数曲线在误差允许的范围以内，并通过实验验证了本文算法的有效性。

2 B样条曲线的基本概念和性质

2.1 B样条曲线定义

B样条曲线是一种描述曲线曲面造型技术，在计算机辅助设计中应用广泛。一条k次B样条曲线表达式[2]：

其中，Pi(i=0，1，…，n)是曲线的控制顶点，控制顶点顺序连成的折线段称为B样条曲线的控制多边形。Ni，k(u)称为k次B样条曲线的基函数，由德布尔递推公式计算如下：

2.2 型值点集选择准则

对于任意连续小线段而言，不能保证用一条符合误差约束的B 样条拟合，因此，需要找出拟合中容易产生超过拟合误差的转折部分，并将其设置为不可拟合的断点，其余部分拟合成平滑的样条曲线。为了提高小线段的拟合率，本文采用相邻线段间的转角，双边弓高误差[11]以及凹凸性准则，用于拟合时筛选能够被拟合的型值点集。

在拟合连续小线段时，为了获得较好的精度，通常对相邻线段采用圆弧近似法来计算拟合的弓高误差。如图1所示，任意相邻两段线段Qi-1Qi和QiQi+1，线段长度分别为li和li+1，所形成的转角为θi，通过判断双边弓高误差εmax=max(δi，δi+1)是否超过给定值以及θi方向变化，来选择适合拟合的连续小线段。

图1 相邻小线段的圆弧近似示意图

2.2.1 转角断点判据

当相邻路径段转角偏转较大时，路径段局部容易形成尖角，导致拟合后的曲线误差较大。为了合理控制拟合精度，相邻路径的夹角应满足：

上式中，θcircle为夹角临界值。当θi≤θcircle时，容易导致轮廓超过限定误差εmax。此时，应将∠Qi-1QiQi+1作为尖角来处理，并对尖角对应的点Qi作为不可拟合的断点。

2.2.2 双边弓高误差判据

为了保证样条高精度拟合效果，需要检测拟合点Qi前后的弓高误差是否超过了最大弓高误差εmax。如果超过εmax，需要将Qi点设置为不可拟合的断点，否则将点Qi放入可拟合的型值点集中。双边误差判据公式如下[12]：

2.2.3 凹凸性准则判据

为了保证样条拟合的凸包性，避免曲线拟合出现拐角或振荡现象，在上述两个准则的基础上，本文增加了对连续小线段凹凸性检测。采用的方法是记录转角θi偏转方向，即：

通过对连续小线段转角θi符号变化的判断，将符号变化的点Qi作为拟合分断点，即可获得具有凹凸特性的拟合点集，从而保证各段拟合曲线的凸包性。

2.3 基于插值拟合的B样条曲线

由式(1)可知，B 样条曲线拟合需要确定节点矢量U={u0，u1，…，un+k+1}和Pi(i=0，1，…，n)控制顶点。而节点矢量和控制顶点求解与2．2 筛选出的每一段型值点集Qi(i=0，1，…，m)有关[2]。

2.3.1 节点矢量计算

为了使一条k次B 样条曲线通过一组型值点Qi(i=0，1，…，m)，其反算过程一般要求曲线的首末端点分别与首末型值点一致，内型值点依次与k次B样条曲线定义域内的节点一一对应，即Qi对应节点值ui+k。因此，B 样条将有m段，所求的k次B 样条插值曲线将有n+1 个控制顶点Pi(i=0，1，…，n)，其中，n=m+k-1。节点矢量U在首末端点处具有k+1阶重复度，相应的节点矢量为：

对于开曲线包括非周期闭曲线，通常两端点取重复度k+1，即u[uk，un+1]=[0，1]。于是u0=u1=…=uk=0，un+1=un+2=…=un+k+1=1，只需计算uk+1，uk+2，…，un内节点矢量。

为了确定与型值点Qi(i=0，1，…，m)相对应的参数值ui+k(i=0，1，…，n)，需要对型值点进行参数化处理。通常对型值点参数化处理有均匀参数化法，向心参数化法和累积弦长参数化法。对于均匀参数化法，当数据分布不均匀时，拟合会产生严重变形。向心参数化法主要适用于数据点急转弯变化的场合。累计弦长参数化法是一种常用的节点矢量计算方法，该方法使得每个节点区间长度与对应曲线上两点之间的弦长相对应，保证了拟合曲线良好的光顺性。因此，本文采用累计弦长参数化方法。

当k=3 时，B 样条曲线而言可达到C2连续，而更高次数的B样条将带来更加复杂计算。因此，本文所述B样条均指三次B样条曲线。采用累计弦长法计算相应节点矢量为：

根据B样条曲线的定义可知，其定义区间不包括节点矢量构成的前三个和后三个参数区间。因此，在端点处存在四重节点U=[u0=u1=u2=u3=u0，un+3=uk(k=1，2，…，(n-1)，un+3=un+4=un+5=un+6=un=1]。从而保证了选取的节点矢量与弦长有良好的对应关系，同时也使得对每个型值点所产生影响的控制点不会完全相同。

2.3.2 控制点的反算

给定型值点Qi(i=0，1，…，m)，节点矢量U=[0，0，0，0，uk+1，uk+2，un，1，1，1，1]，采用切点矢量的边界条件，反算Pi(i=0，1，…，n)为所求控制顶点，属于式(1)的逆过程。

根据三次B样条插值条件：

三次样条的首末节点重复度为4，所以B 样条曲线的首末控制顶点与首末型值点相同，即Q0=P0，Qm=Pn。将上述线性方法改写为矩阵形式：

通过高斯消元法即可求出全部未知控制点。

2.4 误差控制

误差控制是连续小线段拟合后实现高精度加工的关键步骤。定义连续小线段允许最大拟合拟合误差εallowmax和平均误差εavg，C(u)表示连续小线段拟合的参数曲线。如果直接求解连续小线段到C(u)的最大距离εmax和平均误差εavg，一般没有解析解，只能通过迭代算法求解，难以满足数控系统对实时性的要求；另一种简化方法是先计算曲线C(u)与拟合型值点间组成的每一个线段之间的最大弦差εi，max和平均误差εi，avg，后计算所有线段之间最大误差值εmax=max(εi，max)和平均值εavg=(1/n)εi，avg分别为拟合曲线C(u)最大近似误差参考值和平均误差参考值。其中，每一个线段之间的最大弦差εi，max和平均误差εi，avg的计算过程如图2所示。在曲线段C(ui+k)和C(ui+k+1)之间将u[ui+k，ui+k+1]等分成若干等分，得到一系列C(u)的上的点集，计算这些点到对应直线段Qi-1Qi的距离εi，将其最大距离作为εi，max和平均距离作为εi，avg，以此来控制拟合曲线精度。

图2 误差分析示意图

3 仿真及实验

3.1 B样条拟合流程

依据以上方法，针对连续小线段建立完成的拟合流程图如图3所示。

图3 连续小线段拟合流程图

该算法首先提取连续小线段端点作为待拟合型值点集，并采用转角，双边弓高误差以及凸包特性判别拟合分段点，见图4(b)；接着，统计分段点间的型值点数目m，判断拟合点数是否满足三次B样条拟合的最小数目3。对满足拟合条件的型值点，设定最大拟合点数M来控制拟合曲线计算规模，避免计算超时的情况；然后，采用一条三次B 样条拟合这些型值点。而对于不能拟合的型值点，采用直线直接连接。最后，采用误差控制模块，对拟合曲线的误差进行分析。若拟合曲线与原连续小线段的误差没有超过精度限制，则输出拟合误差结果和拟合曲线参数。否则，将超过误差最大的点设置为断点，针对新的断点重新拟合和误差分析，至到所有型值点拟合的样条曲线在误差允许的范围内。

3.2 拟合效果验证

以“蝴蝶”图形为例，对本文所提算法的拟合效果进行验证，其验证结果如图4所示。图4(a)为数控编程NCG01绘制的“蝴蝶”原始轮廓图，实验设定最大拟合误差εallowmax=20um，拟合转角θcircle=40°时，采用本文算法拟合后的“蝴蝶”图形如图4(b)所示，其中，(b)图中的“o”表示分段点。对比图4(a)和图4(b)，可以发现拟合后的“蝴蝶”非常光顺平滑，说明本文算法的有效性。

图4 “蝴蝶”轮廓图及其分段点和拟合曲线

为了进一步验证本文算法的有效性，实验对比了“蝴蝶”拟合前后的数据，见表1。拟合前，“蝴蝶”轮廓由1000 段长短不一的小线段首尾连接而成的，采用本文算法拟合后，“蝴蝶”轮廓仅需38条三次B样条即可表示，小线段数量减少了96．20%，并且拟合后的曲线长度分别较拟合前小线段的最大路径长度整体提升96．17%，最短长度提升98．89%，平均路径长度提升96．19%。说明本文算法对减少小线段数量，使得拟合的路径变长的效果显著。并且，整体拟合时间只花费了0．2919秒，实时性好，可以实现小线段的在线拟合。

表1 “蝴蝶”轮廓拟合数据

lmax，lmin，l平均分别为最大，最短线段长度和平均线段长度。

3.3 误差分析

按照2．4的误差控制模块，将图4(b)“蝴蝶”拟合每一段B 样条曲线，按△u=0．001 下进行等分密化扫描，计算“蝴蝶”拟合图形的最大误差和平均误差，其结果如图5 所示。图中“*”表示每一段B 样条拟合的误差值。由图5(a)可知，坐标(1，16．7561)表示最大误差在第一段拟合曲线上，其值为εmax=16．7561um＜εallowmax=20um，小于拟合允许误差。由图5(b)可知，坐标(38，0．55336)表示最大平均误差在最后一段，其值为εavg=0．55336um，其值远远小于拟合允许误差，说明拟合的曲线非常贴合原始路径。

图5 “蝴蝶”拟合最大误差和平均拟合误差

3.4 实用性分析

为了进一步验证本文算法的实用性，试验在HAN’S801激光切割机床上进行，试验台与工控机相连，进给轴最大速度30m/min，最大加速度为1．5g，对“蝴蝶”轮廓分别采用小线段加工方法和本文所提三次B 样条拟合加工方法进行切割。为了更加清晰地描述实际切割的细节部分，我们使用“黑色”矩形标签框任意挑选切割工件边缘一个区域，并放大显示在对应图像的右侧。图6(a)为小线段直接加工的“蝴蝶”图形，6(b)为6(a)“蝴蝶”加工细节部分。图6(c)为采用本文方法拟合后加工的“蝴蝶”图形，6(d)为6(c)“蝴蝶”加工细节部分。通过对比，不难发现，图6(b)切割边缘毛刺较明显，边缘切割质量差，图6(d)切割边缘较平滑，表面切割质量明显优于图6(a)的方法。

图6 不同算法下的“蝴蝶”加工效果图

表2给出了不同方法拟合“蝴蝶”轮廓后，加工时间和速度。直接采用小线段加工“蝴蝶”轮廓，实际加工时间约为9s，而采用本文所提方法拟合后，由于平均加工速度的提高，实际加工时间减少到了2s，加工效率提高了78%，从实际加工速度与指令加工速度的均方差可以看出，小线段加工过程的速度波动较大，而采用B样条曲线加工的速度波动则小的多，而且小线段加工过程中的最大加速度也比B样条曲线加工过程中的要大。在实际加工过程中，样条曲线拟合的目的是在控制加工轮廓在允许范围内，尽量拉伸单段加工轨迹的长度，提高单段轨迹加工速度，从而提高整个轮廓的加工效率。

表2 不同算法加工“蝴蝶”时间和速度

4 结束语

本文提出的连续小线段的高精度B样条拟合方法，通过将连续小线段拟合成参数曲线来实现数控加工路径变长，同时保证参数拟合曲线精度可控。采用误差对加工路径进行转角和双边弓高检验，快速确定拟合型值点集。运用凹凸性检测准则，保证拟合样条的凸包性。运用误差控制模块，保证拟合的参数曲线控制在精度允许范围以内。实验结果表明，所提出的算法可以将路径数量减少96．20%，同时，在误差控制模块下，保证每一段参数曲线的拟合误差不超过精度限制，而且整个求解过程无迭代，实时性好。通过在数控机床上的切割试验，也进一步验证了本文方法的实用性。