本发明涉及高质量加工中参数样条曲线的拟合及椭圆弧段的识别与合并，属于数控加工技术领域。
背景技术：
椭圆是圆锥曲线中的一种，在工件加工程序中较为常见。由于椭圆弧长的计算涉及到椭圆积分问题，且其原函数不能用有限形式表示，所以无法直接对椭圆进行精确插补。又因其具有较强的工程实用性，国内外的学者对此进行了广泛的研究。
如映射法，将椭圆映射到某一平面，使其投影为圆，通过对圆进行插补，然后利用坐标变换得到椭圆上的插补点。离心角增量法，在弧长增量相对椭圆周长非常小的情况下，通过下一插补点与当前插补点和离心角增量的关系进行插补。上述方法虽简单可行，但都存在轮廓误差和速度波动无法控制的问题。于是提出了等误差方法，虽然提高了椭圆的加工精度，但仍无法有效地控制速度波动。因此大量学者对数值方法进行了研究，通过较为精确地推导出插补点与弧长的关系，来有效地降低速度波动，但这些方法因涉及到多次迭代与复杂数值计算，实时性较差，所以均无法应用于数控系统。
随着设计和制造技术的发展，越来越多的人使用计算机辅助设计(cad)系统来进行复杂零部件的设计。但由于绝大多数的数控系统不支持参数样条数据的传输，通常利用计算机辅助制造(cam)系统将cad设计的自由曲线或曲面，在特定的容差范围内用一系列折线去覆盖，从而生成由大量指令点构成的数控加工程序。
通过对微小直线段进行数控加工，可以不用考虑复杂的椭圆弧长计算问题。但由于加速度和加加速度的频繁变化，会引起机床的震动，降低加工的效率和质量。而现在对于微小线段高速加工的研究主要分为两种方法。一种是在相邻微小线段的拐角处插入过渡样条曲线。例如，通过在拐角处插入三次b样条曲线或四次bézier曲线的方式，提高拐角处的速度，从而提高加工效率，但由于插补点在直线段和过渡样条曲线段上的循环出现，使得插补步长不一致，从而导致加工速度的波动，如果指令点越密集，速度的波动就会越频繁。另一种则是通过插值或者逼近的方式将离散的指令点拟合成平滑的加工路径。例如，将由连续微小线段指定的加工路径转化为由二次bézier曲线表示的加工路径，通过对bézier曲线插补来实现对自由曲线的高速高精加工。虽然这种方式能较好地逼近原设计曲线，但由于拟合曲线较复杂，无法精确计算插补步长对应的插补参数，导致加工速度波动，降低了加工精度。
技术实现要素：
为了克服直接插补算法和现有小线段插补算法不能高精加工和存在速度波动的不足，本发明的目的是提供一种参数样条拟合方法，通过参数样条曲线特征，识别加工路径中的椭圆弧段，通过对椭圆弧段进行插补，来保证加工的精度，同时降低加工速度的波动。
本发明解决其技术问题所采用的技术方案是：一种适用于高速高精加工的椭圆弧平滑压缩插补算法，包括以下步骤：
识别连续微小线段加工区域；
在连续微小线段加工区域中选取型值点并对其进行拟合得到二次有理bézier曲线；
根据曲线特征识别椭圆弧及其几何形式转换；
将属于同一椭圆的相邻椭圆弧合并，得到插补曲线；
在插补曲线上进行椭圆弧插补。
所述进行拟合包括以下步骤：
2-1)求取二次有理bézier曲线的权重；
2-2)根据权重通过求取平均值得到二次有理bézier的拟合曲线。
所述权重通过下式得到：
其中，w1为权重，p0为首型值点、p2为末型值点、p1为控制点、p为指令点，、u为p0q与qp2的比值，q为以p1为投影中心，直线段[p0p2]投影到二次有理bézier拟合曲线的投影点。
所述根据权重通过求取平均值得到二次有理bézier的拟合曲线包括以下步骤：
4-1)根据权重的求取公式得到求取二次有理bézier曲线对应的权值wk及相应的肩点坐标sk＝wk/(1+wk),k＝i+1,...,j-1；
4-2)通过对sk求平均值得到平均肩点坐标s和平均权值w＝s/(1-s)；
4-3)根据p0，p1，p2和w确定型值点qi与qj之间二次有理bézier的拟合曲线。
所述根据曲线特征识别椭圆弧及其几何形式转换包括以下步骤：
5-1)当|p0p1|！＝|p1p2|且0＜w1＜1时，则该曲线为椭圆弧；
5-2)通过以下步骤得到椭圆弧的几何信息：
ci(u)为加工路径中某段的标准型二次有理bézier曲线，且其对应于一个椭圆弧段，根据ci(u)可得到该椭圆弧段的起点坐标p_start和终点坐标p_end；该椭圆弧段的中心坐标为，
p1+ε(s+t)
该椭圆弧段的长、短半径分别为，
假设λ2＞λ1＞0，且是如下二次方程的根，
2δλ2-(kη+4β)λ+2(k-1)＝0
δ＝|s×t|2,η＝|s-t|2，β＝s·t
该椭圆弧段长轴上的两个点为，
q1＝p1+(ε+r1x0)s+(ε+r1y0)t
q2＝p1+(ε-r1x0)s+(ε-r1y0)t
根据q1、q2，可以得到长轴的斜率kl，和长轴与x轴正半轴的夹角d_kl；以长轴为x′轴，短轴为y′轴建立局部坐标系，利用q1、q2、中心坐标、椭圆弧段起点坐标和椭圆弧段终点坐标，求得椭圆弧段起点与中心连线、椭圆弧段终点与中心连线分别相对于局部坐标系x′正半轴的起始角度d_start、终点角度d_end，以及椭圆弧的方向turn；
当根据曲线特征识别椭圆弧及其几何形式转换结束后，得到一个椭圆弧段数组ellipse_arcs[]，数组中每个数据的结构如下：
本发明具有以下有益效果及优点：
1.本发明方法控制简单，可以有效地降低椭圆弧加工速度的频繁波动，实现椭圆弧的高质量加工。
2.压缩量大，平滑度高。本发明方法能够识别出数控加工程序中的椭圆弧，并能用椭圆弧的几何形式进行表示，极大地减少程序段的数量，同时椭圆弧的几何表示形式提高了加工路径的光滑度。
3.加工精度和加工效率高。本发明方法用几何形式的椭圆弧进行插补，能够较为精确地计算出弧长对应的插补参数，降低了计算的复杂度和加工中速度的频繁波动，提高了加工质量和加工效率。
附图说明
图1为本发明的方法流程图
图2为连续微小线段加工区域的识别示意图；
图3为局部曲率最大值的识别示意图；
图4为二次有理bézier曲线示意图；
图5为型值点切矢的计算原理图；
图6为椭圆弧几何信息的计算原理图；
图7为区域及变量的划分示意图。
具体实施方式
下面结合实施例对本发明做进一步的详细说明。
本发明是一种适用于高速高精加工的椭圆弧平滑压缩插补算法，根据双弓高误差限制将加工路径划分为非连续微小线段加工区域和连续微小线段加工区域。对于非连续微小线段加工区域，直接在相邻指令点构成的直线段上进行插补计算，以保证加工精度。对于连续微小线段加工区域，根据离散指令点的曲率值，对曲率极值点和拐点进行拟合，将折线加工路径转化为平滑的二次有理bézier曲线加工路径；然后，利用二次有理bézier曲线特征识别出椭圆弧，并转换为几何形式；最后，将相邻椭圆弧段合并后，进行插补计算。
本发明提供了一种适用于高速高精加工的椭圆弧平滑压缩插补算法，包括以下步骤：
1.加工区域的识别，根据双弓高误差判断条件，将相邻两个不满足条件的点及它们之间的指令点称为连续微小线段加工区域。
2.型值点的选取，对于连续微小线段加工区域中的指令点，通过离散点曲率计算公式,计算其曲率值，并根据相邻指令点的曲率值和判断条件,找出加工路径中局部曲率最大值点和拐点,将连续微小线段加工区域的两端点、局部曲率最大值点和拐点标记为型值点。
3.型值点的拟合，在保证加工精度的条件下，根据型值点的坐标值和单位切矢，将指令点指定的折线加工路径转化为平滑的二次有理bézier曲线加工路径。
4.椭圆弧的识别，根据二次有理bézier曲线的特征识别出椭圆弧段,并将椭圆弧的二次有理bézier形式转换为几何形式。
5.椭圆弧段的合并，对相邻的椭圆弧段，根据其几何信息，识别他们是否属于同一椭圆，对于属于同一椭圆的相邻椭圆弧段进行合并。
6.椭圆弧的插补，通过在几何形式的椭圆弧上进行插补计算，实现椭圆弧的高速高精加工。
如图1所示，本发明提出了一种适用于高速高精加工的椭圆弧平滑压缩插补算法，解决了小线段插补的问题，方法由加工区域的识别、型值点的选取、型值点的拟合、椭圆弧的识别与几何形式的转换、椭圆弧段的合并和椭圆弧的插补6个部分组成，提高了加工的质量和效率。
根据双弓高误差判断条件，将相邻两个不满足条件的点及它们之间的指令点称为连续微小线段加工区域。
对于连续微小线段加工区域中的指令点，通过离散点曲率计算公式,计算其曲率值，并根据相邻指令点的曲率值和判断条件,找出加工路径中局部曲率最大值点和拐点,将连续微小线段加工区域的两端点、局部曲率最大值点和拐点标记为型值点。
保证加工精度的条件下，根据型值点的坐标值和单位切矢，将指令点指定的折线加工路径转化为平滑的二次有理bézier曲线加工路径。
根据二次有理bézier曲线的特征识别出椭圆弧段，并将椭圆弧的二次有理bézier形式转换为几何形式。
对相邻的椭圆弧段，根据其几何信息，识别他们是否属于同一椭圆，对于属于同一椭圆的相邻椭圆弧段进行合并。
椭圆弧的插补通过在几何形式的椭圆弧上进行插补计算，实现椭圆弧的高速高精加工。
本发明的具体步骤如下：
1.连续微小线段加工区域的识别
如图2所示，pi-1、pi和pi+1为顺序相邻的三个指令点，l1、l2为小线段的段长，θ为小线段间的拐角，双弓高误差判断条件如下所示，
其中，δ1、δ2分别为小线段pi-1pi和pipi+1的弓高误差，φ1为opi-1和opi夹角的一半。φ2为opi和opi+1夹角的一半。
如果δ1或δ2大于设定的最大弓高误差值δmax，则pi为断点；若两个相邻的断点间存在指令点，那么两个断点连同它们之间的指令点就构成了一个连续微小线段加工区域。
2.型值点的选取
为了减少对连续微小线段加工区域的拟合次数、增加程序段的压缩量，通过如下三个步骤来选取型值点。
(1)将连续加工区域的开始点和结束点，即断点，标记为型值点。
(2)将局部曲率最大值点标记为型值点。
如图2所示，三个相邻指令点的坐标为pi-1(xi-1,yi-1)、pi(xi,yi)和pi+1(xi+1,yi+1)，离散指令点pi的曲率值由下式确定,
其中δpi-1pipi+1为带符号的三角型面积由下式确定,
假设kl为pi左边的局部曲率最小值，kr为pi右边的局部曲率最小值，如果满足下面两个条件，则pi为局部曲率最大值点，如图3所示。
1)|ki|＞|kl|并且|ki|＞|kr|
2)|ki|-|kl|≥δf或者|ki|-|kr|≥δf，δf为设定的最大曲率差值。
(3)将加工路径弯曲方向改变的点，即拐点，标记为型值点。
利用第(2)步计算出的指令点曲率值进行判断，如果ki-1ki＞0并且kiki+1＜0，ki-1、ki+1分别为指令点pi-1和pi+1的曲率值，则pi为拐点。
3.型值点的拟合
对于连续微小线段加工区域中的n个型值点，可用一条n-1段的二次有理bézier曲线进行拟合，以达到平滑压缩加工路径的目的。
标准型的二次有理bézier曲线如下所示,
其中，p0、p1和p2为控制点，w1为权值，u为变量，范围为[0，1]。
由图4可知，当给定首末端点p0和p2，以及这两点处的切线方向t0和t2时，可以很容易地通过直线[p0t0]和[p2t2]的交点得到p1，再给定一点p就能唯一确定该曲线，因而确定了w1。
将所求的曲线看成是由点p0，p1和p2确定的抛物线的投影，p1为投影中心。如图4所示，将直线段[p0p2]投影到要求的曲线上，则点p和q为相对应的投影点。令w1＝0，得到直线段l(u)＝[p0p2]，即
l(u)是点p0和p2的凸组合，因此|p0q|和|qp2|的比值为u2:(1-u)2，从而推出
将u和p带入二次有理bézier曲线，得到w1，从而得到所求曲线。
假定型值点qi与qj之间由指令点qi+1,qi+2,…,qj-2与qj-1指定的折线加工路径组成，且型值点切矢已知，则型值点qi与qj之间的拟合详细步骤如下所述。
(1)构造以p0＝qi，p1和p2＝qj为控制点的拟合曲线，使其分别插值于qk(k＝i+1,...,j-1)。这将产生中间的权值wk及相应的肩点坐标sk＝wk/(1+wk),k＝i+1,...,j-1；
(2)通过对sk求平均值得到逼近曲线的肩点坐标，即
则中间的权值为w＝s/(1-s)；
(3)根据p0，p1，p2和w即可确定型值点qi与qj之间的拟合曲线。
4.型值点切矢的计算
由于数控加工程序中并不提供指令点处的切矢量，但可以通过型值点以及它周围的四个指令点来计算型值点处的切矢量。如图5所示，qi为型值点，qi-2、qi-1、qi+1、qi+2为它周围的四个指令点，qi的单位切矢表达式如下
由于计算公式中未使用参数u，所求的单位切矢只能被看作切矢的方向。利用下式，即可得到t0，t1和tn-1，tn。
q0＝2q1-q2,q-1＝2q0-q1
qn+2＝2qn+1-qn,qn+1＝2qn-qn-1
5.拟合精度的控制
虽然曲线段经过型值点qi与qj，但不能保证它到qi与qj之间的所有指令点的距离都满足最大拟合误差。因此，将qi与qj之间的所有指令点投影到拟合的曲线上，检查它们的偏差是否都在允许的误差范围内，即小于等于系统设定的最大轮廓误差δc，如果满足误差要求，进行下一段的曲线拟合；否则，将偏差最大的指令点设置为新的型值点，用新的型值点重新进行曲线拟合，再次检验拟合精度，重复这个过程，直至满足误差要求。
6.椭圆弧的识别与几何形式的转换
由于标准型的二次有理bézier曲线只有一个权因子，其表达能力比非有理形式的bézier曲线强，可以表达多种曲线，当|p0p1|！＝|p1p2|且0＜w1＜1时，则所表达的曲线为椭圆弧。
如果ci(u)为加工路径中某段的标准型二次有理bézier曲线，且其对应于一个椭圆弧段，如图6所示，根据ci(u)可得到该椭圆弧段的起点坐标p_start和终点坐标p_end。由椭圆的标准型二次有理bézier形式与其几何形式的关系可知，该椭圆弧段的中心坐标为，
p1+ε(s+t)
该椭圆弧段的长、短半径分别为，
假设λ2＞λ1＞0，且是如下二次方程的根，
2δλ2-(kη+4β)λ+2(k-1)＝0
δ＝|s×t|2,η＝|s-t|2，β＝s·t
该椭圆弧段长轴上的两个点为，
q1＝p1+(ε+r1x0)s+(ε+r1y0)t
q2＝p1+(ε-r1x0)s+(ε-r1y0)t
根据q1、q2，可以得到长轴的斜率kl，和长轴与x轴正半轴的夹角d_kl；以长轴为x′轴，短轴为y′轴建立局部坐标系，利用q1、q2、中心坐标、椭圆弧段起点坐标和椭圆弧段终点坐标，求得椭圆弧段起点与中心连线、椭圆弧段终点与中心连线分别相对于局部坐标系x′正半轴的起始角度d_start、终点角度d_end，以及椭圆弧的方向turn；x0、y0为椭圆弧段起点坐标。
当连续微小线段加工区域完成椭圆弧的识别和几何形式转换后，会得到一个椭圆弧段数组ellipse_arcs[]，其中椭圆弧段的数据结构如下：
7.椭圆弧的合并
为了提高椭圆弧加工的效率和精度，将连续微小线段加工区域内的椭圆弧段合并，如果相邻的椭圆弧段具有相同的中心center、长半径r_l、短半径r_s、方向turn、长轴与x正半轴的夹角d_kl，并且上一椭圆弧段的终点坐标p_endi-1与下一椭圆弧段的起点坐标p_starti相等，则它们属于同一椭圆弧上连续的椭圆弧段，通过修改上一椭圆弧段的终点角度d_endi-1和终点坐标p_endi-1将它们合并为一个椭圆弧段。重复上述步骤，直至椭圆弧段不能再合并。
8.椭圆弧的插补
8.1椭圆弧长的计算
假设椭圆方程为
以第一象限y为变量为例，
根据高斯勒让德求积公式，可将上式转换为，
其它象限x或y为变量的情况，其弧长的推导方法相同。
8.2弧长与变量的关系
以第一象限y为变量，从点e0(a,0)到点
(1)在变量y的区间
(2)根据8.1的弧长公式分别计算yi+1＝0、y1、y2和
(3)将(s0,0)、(s1,y1)、(s2,y2)和
y＝a0+a1s+a2s2+a3s3，s∈[s0,s3]
其中，a0～a3为系数。其它象限x或y为变量的情况，其弧长与变量关系的推导方法相同。
8.3椭圆弧的插补
通过上述过程便可将数控加工路径中的椭圆弧识别出来，并用几何的形式表示。对于识别出来的椭圆弧段，按照四个区域对其进行划分，并分别计算出各个区域内弧长与变量的关系，然后利用这些弧长与变量的关系来计算第i个插补周期的插补点坐标(xi′,yi′)，具体步骤如下所示。
假定插补周期为t，p0′(x′i-1,y′i-1)和si-1分别为第i-1个插补周期的插补点坐标与弧长，第i个插补周期的进给速度为vi-1，则弧长增量为δs＝vi-1t，第i个插补周期的弧长为si＝si-1+δs，将si代入弧长与变量的关系中，得到对应变量xi′或yi′的值，再将其代入椭圆方程中，从而求得第i个插补周期的插补点坐标p′(xi′,yi′)。