CN104597847A - 基于Akima样条曲线拟合的前瞻插补方法 - Google Patents

Abstract

本发明是基于Akima样条曲线拟合的前瞻插补方法。Akima样条拟合首先利用向心参数化方法计算每个插补数据点对应的参数，然后采用阿克玛五点插值法将数据点拟合成光滑的Akima样条曲线。前瞻插补根据加工精度的要求以及数控***的编程速度和最大加速度信息，通过前瞻插补计算，向实时插补返回不同的数值。实时插补则首先根据前瞻插补返回的值调整实时进给速度，然后根据已知的Akima样条曲线方程和插补周期等信息，创建构造方程，预测下一个插补参数的迭代初值后，再利用牛顿迭代法精确计算插补参数。本发明方法采用Akima样条曲线拟合待加工的数据点，执行效率高，光滑性好，在保证了插补的实时性的同时，可以提高拟合精度，具有很高的准确性和灵活性，可以在较大程度上还原工件的实际状态。

Description

基于Akima样条曲线拟合的前瞻插补方法

技术领域

本发明属于数控加工技术领域，具体涉及基于Akima样条曲线拟合的前瞻插补方法。

背景技术

目前，复杂曲面零部件的生产和制造在航空航天、汽车、轮船、刀具和模具等行业具有特别重要的意义。对于复杂曲面、曲线的加工，传统数控***采用离散化的微小直线段或圆弧段逼近加工曲线，这种方法容易造成进给速度轮廓的不连续和波动，破坏工件表面的光滑性，而且产生的大量程序增加了CAD/CAM和计算机数控***CNC之间的通信负担，影响了插补的实时性。样条曲线虽然在CAD***中得到广泛应用，但是在CNC***中的发展相对滞后，只有一些高档数控***如SIEMENS、FANUC等实现了样条曲线插补功能。FANUC***采用G01指令点生成平滑的样条曲线，并通过对所生成的样条曲线进行插补来提高自由曲面的加工质量。SIEMENS的840D高级编程手册中则提出了压缩器的概念，通过连接一系列G01指令点，并将其压缩以形成平滑的样条曲线进行插补，具体有A(Akima)样条插补、B样条插补(NURBS插补)和C样条插补，但是这些技术和方法一直是对外保密的。

由于NURBS曲线是国际标准化组织规定的CAD/CAM的数据交换标准，且NURBS曲线的曲率变化平稳，因此针对NURBS样条曲线插补方法的研究比较广泛。但是NURBS曲线的拟合过程和数学计算比较复杂，难以满足插补的实时性。另外一种方法是将连续微小线段指定的加工路径拟合成Bézier曲线，并通过对Bézier样条曲线插补来提高工件表面的加工质量，但是这种方法需要计算额外的Bézier曲线控制点，计算也比较复杂。还有一种方法是用Cubic样条曲线拟合刀具轨迹进行插补，但是Cubic样条曲线本身具有产生摆动的倾向，改变一个数据点可能会对许多部分产生影响。

发明内容

针对现有方法中存在的上述不足，本发明要解决的技术问题是提供了一种基于Akima样条拟合的样条曲线插补方法。

本发明为实现上述目的所采用的技术方案是：一种基于Akima样条曲线拟合的前瞻插补方法，包括以下步骤：

Akima样条曲线拟合，首先利用向心参数化方法计算每个插补数据点对应的参数，然后采用阿克玛五点插值法将插补数据点拟合成光滑的Akima样条曲线；

前瞻插补，根据加工精度要求的最大弦高误差和数控***的编程进给速度计算当前插补参考进给速度，使刀具按照所述插补参考进给速度前瞻插补，之后的每一步插补都以最大的逆向加速度进行减速，计算每次插补的参考速度，并根据不同的情况向实时插补返回数值或者插补到额定前瞻数；

实时插补，根据前瞻插补返回的值调整实时进给速度，根据已知的Akima样条曲线方程和插补周期创建构造方程，并用一阶泰勒展开式，预测下一个插补参数的迭代初值，然后利用该初值，应用牛顿迭代法计算插补参数，直到满足迭代终止条件。

所述Akima样条曲线，在参数[u_i,u_i+1]之间的形式为：

R(u)=A+B(u-u_i)+C(u-u_i)²+D(u-u_i)³

其中R(u)为Akima样条曲线，A，B，C，D为Akima曲线的系数矩阵，u为变量参数，u_i为数据点P_i对应的参数。

所述向心参数化方法的公式为：

$\left\{\begin{array}{c}{u}_0=0\\ u_i=u_{i-1}+\sqrt{|P_i-P_{i-1}|}i=\mathrm{1,2}Ln-1\end{array}\right.$

其中，u₀为初始参数，P_i为第i个数据点，u_i为数据点P_i对应的参数，P_i-1为第i-1个数据点，u_i-1为数据点P_i-1对应的参数。

所述阿克玛五点插值法用于计算Akima样条曲线的系数，其公式为：

$\left\{\begin{array}{c}{A}_{\mathrm{Cd}}=C{d}_i\\ B_{\mathrm{Cd}}=g_{\mathrm{Cd},i}\\ C_{\mathrm{Cd}}=\frac{3u_{\mathrm{Cd},i}-2g_{\mathrm{Cd},i}-g_{\mathrm{Cd},i+1}}{u_{i+1}-u_i}\\ D_{\mathrm{Cd}}=\frac{-2u_{\mathrm{Cd},i}+g_{\mathrm{Cd},i}+g_{\mathrm{Cd},i+1}}{{(u_{i+1}-u_i)}^2}\end{array}\right.$

其中，Cd表示x,y或者z的坐标符号，若Cd表示x，则A_Cd表示系数矩阵A中的A_x；Cd_i表示点P_i的横坐标x的值，g_Cd,i表示在u_i处R(u)中的Cd量对于u的一阶导数值，其计算公式为：

$g_{\mathrm{Cd},i}=\frac{|u_{\mathrm{Cd},i+1}-u_{\mathrm{Cd},i}|u_{\mathrm{Cd},i-1}+|u_{\mathrm{Cd},i-1}-u_{\mathrm{Cd},i-2}|u_{\mathrm{Cd},i}}{|u_{\mathrm{Cd},i+1}-u_{\mathrm{Cd},i}|+|u_{\mathrm{Cd},i-1}-u_{\mathrm{Cd},i-2}|}$

所述插补参考进给速度V_r(u_i)为：

$V_r\left(u_i\right)=\left\{\begin{array}{cc}{V}_e\left(u_i\right),& V_e\left(u_i\right)<F\\ F,& V_e\left(u_i\right)\&GreaterEqual;F\end{array}\right.$

其中，u_i为当前插补点P_i对应的插补参数，F为数控***的编程进给速度；V_e(u_i)为精度要求下的速度，由以下公式决定：

$V_e\left(u_i\right)=\frac{T}{2}\sqrt{{\mathrm{\&rho;}}_i^2-{({\mathrm{\&rho;}}_i-\mathrm{ER})}^2}$

其中ρ_i为曲率半径，T为数控***的插补周期，ER为加工要求的最大弦高误差。

所述前瞻插补向实时插补返回值由以下情况决定：

若减速插补到额定前瞻数之前能减速到0，则返回0；

若在减速插补的过程中，进给速度大于参考速度，则返回1；

若减速插补到额定前瞻数的整个期间，没有进给速度大于参考速度，则返回0。

所述额定前瞻数为：

$N=\frac{F}{A_{\max }}$

其中，F为数控***的编程进给速度，A_max为数控***最大的加速度。

所述实时插补根据前瞻插补返回的数值调整实时进给速度的规则如下：

若返回值为0，则实时进给速度为参考速度V(u_i)；

若返回值为1，则实时插补从当前插补点以最大的逆向加速度减速，即V(u_i)=V(u_i-1)-A_maxT。

所述构造方程用实际弦长逼近理论弦长的方法构建，其表达式为：

F(ξ)=||R(ξ)-R(u_i)||-V(u_i)T

其中，u表示当前的插补参数，ξ表示要计算的下一个插补参数，F(ξ)即为构造方程，R(ξ)为参数ξ处的Akima方程的值，R(u_i)为参数u处的Akima方程的值。

所述应用牛顿迭代法计算插补参数的迭代公式为：

${\mathrm{\&xi;}}_{k+1}={\mathrm{\&xi;}}_k-\frac{F\left({\mathrm{\&xi;}}_k\right)}{F^{\text{'}}\left({\mathrm{\&xi;}}_k\right)}$

其中，ξ_k表示经过第k次迭代后ξ的值，ξ_k+1表示第k次迭代后ξ的值，F’(ξ_k)表示参数ξ_k处F(ξ)的一阶导数值。

所述迭代终止条件为：

δ_i+1<Δ

k=K_max

其中，δ_i+1为当前的u_i+1对应的速度波动率的值，Δ为速度波动率上限，k为当前迭代的次数，K_max为最大迭代次数。

本发明具有以下优点及有益效果：

1.拟合精度高。本发明方法采用Akima样条曲线拟合待加工的数据点，能够精确穿过每个数据点，拟合精度高，光滑性好，具有相当的准确性和灵活性，可以在较大程度上还原零件的实际状态。

2.保证了加工精度。本发明方法通过前瞻插补，使得刀具在尖锐拐点之前能够充分减速，从而保证了加工的轮廓精度，并且保证了进给速度和进给加速度不会超过***要求的范围，减小了机床震动，能够满足加工精度的要求。

3.平滑度高，加工质量高。本发明方法采用向心参数化方法计算各个数据点对应的参数，考虑了数据点相邻弦线的折拐情况，使得拟合的Akima曲线更加广顺平滑；而且实时插补时采用了牛顿迭代法计算插补参数，保证了插补轨迹的平滑性，减小了速度波动，在一定程度上提高了加工质量。

4.执行效率高。本发明方法采用Akima样条拟合插补数据点，而Akima样条曲线的方程为三次多项式，其数学运算较NURBS曲线要简单的多，保证了插补的实时性。

附图说明

图1为待加工的数据点及其折线图；

图2为本发明方法的算法流程图；

图3为本发明方法加工的弦高误差曲线；

图4为本发明方法加工的进给速度曲线；

图5为本发明方法加工的加速度曲线；

图6为本发明方法同NURBS拟合插补的加工对比图；

图7为本发明方法同NURBS拟合插补的在折点处的加工对比放大图。

具体实施方式

下面结合附图及实施例对本发明做进一步的详细说明。

本发明方法包括Akima样条曲线拟合、前瞻插补和实时插补。

Akima样条拟合首先利用向心参数化方法计算每个插补数据点对应的参数，然后采用阿克玛五点插值法将数据点拟合成光滑的Akima样条曲线。

前瞻插补根据加工精度要求的最大弦高误差和数控***的编程进给速度计算当前插补参考进给速度，并让刀具按照该速度前瞻插补，之后的每一步插补都以最大的逆向加速度进行减速，计算每次插补的参考速度，并根据不同的情况向实时插补返回数值或者插补到额定前瞻数。

实时插补则首先根据前瞻插补返回的值调整实时进给速度，然后根据已知的Akima样条曲线方程和插补周期等信息，创建构造方程，并用一阶泰勒展开式，预测下一个插补参数的迭代初值，然后利用该初值，应用牛顿迭代法，迭代计算插补参数，直到满足迭代终止条件。

所述的Akima样条曲线，在参数[u_i,u_i+1]之间的形式为：

R(u)=A+B(u-u_i)+C(u-u_i)²+D(u-u_i)³

其中R(u)为Akima样条曲线，A，B，C，D为Akima曲线的系数矩阵，u为变量参数，u_i为数据点P_i对应的参数。

所述的向心参数化方法的公式为：

$\left\{\begin{array}{c}{u}_0=0\\ u_i=u_{i-1}+\sqrt{|P_i-P_{i-1}|}i=\mathrm{1,2}Ln-1\end{array}\right.$

所述阿克玛五点插值法计算Akima曲线的系数，其公式：

$\left\{\begin{array}{c}{A}_{\mathrm{Cd}}=C{d}_i\\ B_{\mathrm{Cd}}=g_{\mathrm{Cd},i}\\ C_{\mathrm{Cd}}=\frac{3u_{\mathrm{Cd},i}-2g_{\mathrm{Cd},i}-g_{\mathrm{Cd},i+1}}{u_{i+1}-u_i}\\ D_{\mathrm{Cd}}=\frac{-2u_{\mathrm{Cd},i}+g_{\mathrm{Cd},i}+g_{\mathrm{Cd},i+1}}{{(u_{i+1}-u_i)}^2}\end{array}\right.$

其中，Cd表示x,y或者z的坐标符号，若Cd表示x，则A_Cd表示系数矩阵A中的A_x，Cd_i表示点P_i的横坐标x的值，g_Cd,i表示在u_i处，R(u)中的Cd量对于u的一阶导数值，其计算公式为：

$g_{\mathrm{Cd},i}=\frac{|u_{\mathrm{Cd},i+1}-u_{\mathrm{Cd},i}|u_{\mathrm{Cd},i-1}+|u_{\mathrm{Cd},i-1}-u_{\mathrm{Cd},i-2}|u_{\mathrm{Cd},i}}{|u_{\mathrm{Cd},i+1}-u_{\mathrm{Cd},i}|+|u_{\mathrm{Cd},i-1}-u_{\mathrm{Cd},i-2}|}$

所述的插补参考进给速度V(u_i)由以下公式决定：

$V_r\left(u_i\right)=\left\{\begin{array}{cc}{V}_e\left(u_i\right),& V_e\left(u_i\right)<F\\ F,& V_e\left(u_i\right)\&GreaterEqual;F\end{array}\right.$

其中，u_i为当前插补点P_i对应的插补参数，F为数控***的编程进给速度，V_e(u_i)为精度要求下的速度，由以下公式决定：

$V_e\left(u_i\right)=\frac{T}{2}\sqrt{{\mathrm{\&rho;}}_i^2-{({\mathrm{\&rho;}}_i-\mathrm{ER})}^2}$

其中ρ_i为曲率半径，T为数控***的插补周期，ER为加工要求的最大弦高误差。

所述的前瞻插补向实时插补返回值由以下情况决定：

(1)若减速插补到额定前瞻数之前能减速到0，则返回0；

(2)若在减速插补的过程中，进给速度大于参考速度，则返回1；

(3)若减速插补到额定前瞻数的整个期间，没有进给速度大于参考速度，则返回0。

所述的额定前瞻数N由以下公式决定：

$N=\frac{F}{A_{\max }}$

其中，A_max为数控***最大的加速度。

所述的实时插补根据前瞻插补返回的数值调整实时进给速度规则如下：

(1)若返回值为0，则实时进给速度为参考速度V(u_i)；

(2)若返回值为1，则实时插补从当前插补点以最大的逆向加速度减速，即V(u_i)=V(u_i-1)-A_maxT；

所述的构造方程用实际弦长逼近理论弦长的方法构建，其表达式为：

F(ξ)=||R(ξ)-R(u_i)||-V(u_i)T

所述的应用牛顿迭代法计算插补参数，迭代公式为：

${\mathrm{\&xi;}}_{k+1}={\mathrm{\&xi;}}_k-\frac{F\left({\mathrm{\&xi;}}_k\right)}{F^{\text{'}}\left({\mathrm{\&xi;}}_k\right)}$

其中，ξ_k表示经过第k次迭代后ξ的值，ξ_k+₁表示第k次迭代后ξ的值，F’(ξ_k)表示参数ξ_k处F(ξ)的一阶导数值。

所述的迭代终止条件为：

(a)δ_i+1<Δ

(b)k=K_max

在迭代计算时，只要满足上述条件的(a)或(b)，则迭代终止。

这里δ_i+1为当前的u_i+1对应的速度波动率的值，Δ为速度波动率上限，k为当前迭代的次数，K_max为最大迭代次数。

本实施例将本发明方法在PC上进行仿真验证，所用的编程软件为MicrosoftVisual C++6.0，使用C语言编写程序，这里仿真插补选用的数据点如图1所示。

测试环境的主要技术插补参数如下：

操作***：Microsoft Windows XP

CPU：Pentium(R)Dual-Core

主频：2.93GHz

内存：2G

数控***插补参数如下：

编程进给速度F=100mm/s；

最大加速度A=5000mm/s2；

最大弦高误差ER=0.001mm；

最大速度波动率Δ=1.0×10-4%；

插补周期T=4ms。

本发明方法是基于Akima样条曲线拟合的前瞻插补方法，其流程图如图2所示。

本发明方法包括Akima样条曲线拟合、前瞻插补和实时插补三部分。

Akima样条拟合首先利用向心参数化方法计算每个插补数据点对应的参数，然后采用阿克玛五点插值法将数据点拟合成光滑的Akima样条曲线。

前瞻插补根据加工精度要求的最大弦高误差和数控***的编程进给速度计算本次插补参考进给速度，并假设刀具按照该速度进行前瞻插补，之后的每一步都以最大的逆向加速度进行减速插补，并计算每一步的精度要求下的参考速度，并根据以下情况向实时插补返回数值。

(1)若减速到额定前瞻数之前能减速到0，则向实时插补返回0；

(2)若在减速插补的过程中，出现进给速度大于参考速度的情况，则向实时插补返回1；

(3)若在减速插补到额定前瞻段的整个期间，没有速度大于参考速度，则向实时插补返回0。

实时插补则首先根据前瞻插补返回的值调整实时进给速度，若返回值为0，则实时进给速度为精度要求下的参考速度；若返回值为1，则从当前插补点以最大的逆向加速度减速。在得到进给速度后，根据已知的Akima样条曲线方程和插补周期等信息，根据理论进给弦长逼近实际进给弦长的策略创建构造方程，然后应用一阶泰勒展开式，预测下一个插补参数的初值，然后利用该初值，应用牛顿迭代法，迭代计算插补参数，直到速度波动率能达到要求的范围内或者迭代的次数超过了最大限定次数。

采用本发明方法，对图1所示的数据点进行仿真加工，得到的整个加工的弦高误差曲线如图3所示，速度曲线如图4所示，加速度曲线如图5所示，。为了验证本发明方法拟合刀具轨迹的效果，仿真实验对比了应用Akima样条曲线和NURBS曲线加工图1所示的离散数据点所代表的曲线，实验结果如图6所示，图7是在二者在拐点处的放大图，能够进一步体现二者的拟合效果。

通过分析，可以得到如下结论：

1.本发明方法可以把加工的弦高误差限制在加工精度要求的范围内。本发明方法通过前瞻插补，使得刀具在尖锐拐点之前能够充分减速，从而保证了加工的轮廓精度。从图3可以看出，加工过程中的所有插补点的弦高误差均未超过***要求的最大弦高误差，所以本发明方法能够满足加工精度的要求。

2.并保证插补时的进给速度、进给加速度均不超过数控***的限制。从图4可以看出，所有加工过程的进给速度一直限制在编程速度内，从图5可以看出，加工过程的加速度也限制在数控***的最大加速度内，所以本发明方法可以使机床运行平稳。

3.本发明方法的拟合精度高，平滑程度高，加工质量高。从图6和图7可以看出，用本发明算法加工离散数据点控制的曲线时，加工轨迹能够精确穿过每个数据点，拟合精度高，并且拟合轨迹的光滑性好。通过对比本发明方法和NURBS曲线拟合加工的效果，可以看出，用Akima样条拟合加工轨迹的误差更小，拟合轨迹更加接近数据点的折线图，因此本发明方法具有相当的准确性和灵活性，可以在较大程度上还原零件的实际状态。

4.本发明方法控制简单，执行效率高。本发明方法采用Akima样条拟合插补数据点，而Akima样条曲线的方程为三次多项式，其数学运算较NURBS曲线要简单的多，统计本发明算法的插补时间较NURBS曲线拟合插补的时间少了一个数量级，因此本发明方法保证了插补的实时性。

Claims (11)

1.一种基于Akima样条曲线拟合的前瞻插补方法，其特征在于，包括以下步骤： 

Akima样条曲线拟合，首先利用向心参数化方法计算每个插补数据点对应的参数，然后采用阿克玛五点插值法将插补数据点拟合成光滑的Akima样条曲线； 

前瞻插补，根据加工精度要求的最大弦高误差和数控***的编程进给速度计算当前插补参考进给速度，使刀具按照所述插补参考进给速度前瞻插补，之后的每一步插补都以最大的逆向加速度进行减速，计算每次插补的参考速度，并根据不同的情况向实时插补返回数值或者插补到额定前瞻数； 

实时插补，根据前瞻插补返回的值调整实时进给速度，根据已知的Akima样条曲线方程和插补周期创建构造方程，并用一阶泰勒展开式，预测下一个插补参数的迭代初值，然后利用该初值，应用牛顿迭代法计算插补参数，直到满足迭代终止条件。 

2.根据权利要求1所述的基于Akima样条曲线拟合的前瞻插补方法，其特征在于，所述Akima样条曲线，在参数[u_i,u_i+1]之间的形式为： 

R(u)=A+B(u-u_i)+C(u-u_i)²+D(u-u_i)³

其中R(u)为Akima样条曲线，A，B，C，D为Akima曲线的系数矩阵，u为变量参数，u_i为数据点P_i对应的参数。 

3.根据权利要求1所述的基于Akima样条曲线拟合的前瞻插补方法，其特征在于，所述向心参数化方法的公式为： 

其中，u₀为初始参数，P_i为第i个数据点，u_i为数据点P_i对应的参数，P_i-1为第i-1个数据点，u_i-1为数据点P_i-1对应的参数。 

4.根据权利要求1所述的基于Akima样条曲线拟合的前瞻插补方法，其特征在于，所述阿克玛五点插值法用于计算Akima样条曲线的系数，其公式为： 

其中，Cd表示x,y或者z的坐标符号，若Cd表示x，则A_Cd表示系数矩阵A中的A_x；Cd_i表示点P_i的横坐标x的值，g_Cd,i表示在u_i处R(u)中的Cd量对于u的一阶导数值，其计算公式为： 

。

5.根据权利要求1所述的基于Akima样条曲线拟合的前瞻插补方法，其特征在于，所述插补参考进给速度V_r(u_i)为： 

其中，u_i为当前插补点P_i对应的插补参数，F为数控***的编程进给速度；V_e(u_i)为精度要求下的速度，由以下公式决定： 

其中ρ_i为曲率半径，T为数控***的插补周期，ER为加工要求的最大弦高误差。 

6.根据权利要求1所述的基于Akima样条曲线拟合的前瞻插补方法，其特征在于，所述前瞻插补向实时插补返回值由以下情况决定： 

若减速插补到额定前瞻数之前能减速到0，则返回0； 

若在减速插补的过程中，进给速度大于参考速度，则返回1； 

若减速插补到额定前瞻数的整个期间，没有进给速度大于参考速度，则返回0。 

7.根据权利要求1所述的基于Akima样条曲线拟合的前瞻插补方法，其特征在于，所述额定前瞻数为： 

其中，F为数控***的编程进给速度，A_max为数控***最大的加速度。 

8.根据权利要求1所述的基于Akima样条曲线拟合的前瞻插补方法，其特征在于，所述实时插补根据前瞻插补返回的数值调整实时进给速度的规则如下： 

若返回值为0，则实时进给速度为参考速度V(u_i)； 

若返回值为1，则实时插补从当前插补点以最大的逆向加速度减速，即V(u_i)=V(u_i-1)-A_maxT。 

9.根据权利要求1所述的基于Akima样条曲线拟合的前瞻插补方法，其特征在于，所述构造方程用实际弦长逼近理论弦长的方法构建，其表达式为： 

F(ξ)=||R(ξ)-R(u_i)||-V(u_i)T 

其中，u表示当前的插补参数，ξ表示要计算的下一个插补参数，F(ξ)即为构造方程，R(ξ)为参数ξ处的Akima方程的值，R(u_i)为参数u处的Akima方程的值。 

10.根据权利要求1所述的基于Akima样条曲线拟合的前瞻插补方法，其特征在于，所述应用牛顿迭代法计算插补参数的迭代公式为： 

其中，ξ_k表示经过第k次迭代后ξ的值，ξ_k+1表示第k次迭代后ξ的值，F’(ξ_k)表示参数ξ_k处F(ξ)的一阶导数值。 

11.根据权利要求1所述的基于Akima样条曲线拟合的前瞻插补方法，其特征在于，所述迭代终止条件为： 

δ_i+1<Δ 

k=K_max

其中，δ_i+1为当前的u_i+1对应的速度波动率的值，Δ为速度波动率上限，k为当前迭代的次数，K_max为最大迭代次数。