适用于高速高精加工的椭圆弧平滑压缩插补算法

李 浩，吴文江，陈渌萍，韩文业，郭 安

1(中国科学院大学，北京 100049) 2(中国科学院沈阳计算技术研究所 高档数控国家工程研究中心， 沈阳 110168) 3(沈阳高精数控智能技术股份有限公司， 沈阳 110168) 4(工业和信息化部计算机与微电子发展研究中心 中国软件评测中心，北京 100048)

1 引 言

椭圆是圆锥曲线中的一种，其拥有几何形式，能够被精确地表示.由于椭圆弧长的计算涉及到椭圆积分问题，且其原函数不能用有限形式表示，所以无法直接对椭圆进行精确插补.又因其具有较强的工程实用性，国内外的学者对此进行了广泛的研究.

苏[1]等提出了映射法，将椭圆映射到某一平面，使其投影为圆，通过对圆进行插补，然后利用坐标变换得到椭圆上的插补点.胡[2]等推导出，在弧长增量相对椭圆周长非常小的情况下，下一插补点与当前插补点和离心角增量的关系.上述方法虽简单可行，但都存在轮廓误差和速度波动无法控制的问题.于[3]等提出了等误差方法，虽然提高了椭圆的加工精度，但仍无法有效地控制速度波动.因此大量学者对数值方法进行了研究，通过较为精确地推导出插补点与弧长的关系，来有效地降低速度波动，如胡[4]提出的幂级数展开法、Toshio[5]提出的半参数和双参数转换法、John[6]提出的经验参数迭代法和Toshio[7]提出的逆椭圆积分法，这些方法因涉及到多次迭代与复杂数值计算，实时性较差，所以均无法应用于数控系统.

随着设计和制造技术的发展，越来越多的人使用计算机辅助设计(CAD)系统来进行复杂零部件的设计，然后利用计算机辅助制造(CAM)系统将CAD设计的自由曲线或曲面，在设定的最大轮廓误差范围内，用微小直线段去逼近，从而生成由大量指令点构成的数控加工程序[8].

通过对微小直线段进行数控加工[9,10]，可以不用考虑复杂的椭圆弧长计算问题.但由于加速度和加加速度的频繁变化，会引起机床的震动，降低加工的效率和质量.而现在对于微小线段高速加工的研究主要分为两种方法.一种是在相邻微小线段的拐角处插入过渡样条曲线[11].如Zhang等[12]提出通过在拐角处插入三次参数样条曲线的方法，提高拐角处的速度，从而提高加工效率.但由于插补点在直线段和过渡样条曲线段上的循环出现，使得插补步长不一致，从而导致加工速度的波动，如果指令点越密集，速度的波动就会越频繁.另一种则是通过插值或者逼近的方式将离散的指令点拟合成平滑的加工路径[13].例如Zhang[14]提出利用二次Bézier曲线将连续微小线段拟合成光滑的样条曲线，通过样条插补来实现对自由曲线的高速高精加工.虽然拟合方式能较好地逼近原设计曲线，但由于拟合曲线较复杂，无法精确计算插补步长对应的插补参数，导致加工速度波动，降低了加工精度.

椭圆弧是一种不能被直接精确快速插补的图形.在离散小线段的形式下采用样条插补方式，可以避免复杂的椭圆弧长计算，但会存在计算精度和效率的问题.所以，基于以上研究和分析，本文提出了一种适用于高速高精加工的椭圆弧平滑压缩插补算法.

2 椭圆弧平滑压缩插补算法

本文算法的基本思想是：从由离散小线段构成的加工路径中识别出连续微小线段加工区域，然后寻找到型值点，通过控制拟合误差，将指令点指定的折线加工路径转化为平滑的二次有理Bézier曲线加工路径.然后，根据样条曲线的特征识别出椭圆弧段，将有理形式转换为几何形式，同时合并相邻的椭圆弧段.最后，通过在几何形式的椭圆弧上进行插补，实现高速高精加工.

该算法采用预处理方式，其由7个部分组成：加工区域的识别、型值点的选取、型值点的拟合、椭圆弧的识别、几何形式的转换、椭圆弧段的合并和椭圆弧的插补，其中加工区域的识别主要是从加工路径中寻找可以拟合的连续区；型值点的选取则是寻找连续区中的曲率极值点和拐点；型值点的拟合用于将连续区内的折线加工路径转化成光滑的样条曲线；椭圆弧的识别用于识别样条曲线中的椭圆弧段；几何形式的转换是将椭圆弧段的有理表示形式转换为几何表示形式；椭圆弧段的合并用于合并已识别出的相邻椭圆弧段；椭圆弧的插补用于对几何形式的椭圆弧进行插补计算.

图1 连续微小线段验证Fig.1 Continuous micro-line verification

2.1 加工区域的识别

平滑压缩算法并非针对所有的G01段，只有连续的微小线段才适合于拟合为Bézier曲线.以双弓高误差为判断依据，从由大量微小线段构成的加工路径中识别出连续微小线段加工区域.如图1所示Pi-1、Pi和Pi+1为顺序的三个指令点l1、l2，为小线段的段长，θ为小线段间的拐角，双弓高误差可由式(1)和(2)得到.

δ1=R(1-cosφ1)

(1)

δ2=R(1-cosφ2)=R(1-cos(π-θ-φ1))

(2)

(3)

(4)

如果δ1或δ2大于设定的最大弓高误差值δmax，则Pi为断点；若两个相邻的断点间存在指令点，那么两个断点连同它们之间的指令点就构成了一个连续微小线段加工区域.

2.2 型值点的选取

为了减少对连续微小线段加工区域的拟合次数、增加程序段的压缩量，通过如下三个步骤来选取型值点.

1)将连续加工区域的开始点和结束点，即断点，标记为型值点.

2)将局部曲率最大值点标记为型值点.

如图1所示，三个指令点的坐标为Pi-1(xi-1，yi-1)、Pi(xi，yi)和Pi+1(xi+1，yi+1)，由式(5)可以计算出离散指令点Pi的曲率值[15]，其中ΔPi-1PiPi+1为带符号的三角型面积.

(5)

(6)

假设kl为Pi左边的局部曲率最小值，kr为Pi右边的局部曲率最小值，如果满足下面两个条件，则Pi为局部曲率最大值点，如图2所示.

1)|ki|>|kl|并且|ki|>|kr|

2)|ki|-|kl|≥δf或者|ki|-|kr|≥δf，δf为设定的最大曲率差值.

图2 局部曲率最大值Fig.2 Local curvature maximum

3)将加工路径弯曲方向改变的点，即拐点，标记为型值点.

利用第(2)步计算出的指令点曲率值进行判断，如果ki-1ki>0并且kiki+1<0，则Pi为拐点.

2.3 型值点的拟合

假设连续微小线段加工区域中有n个型值点，为了达到平滑压缩加工路径的目的，使用n-1段二次有理Bézier曲线对其进行拟合.

二次有理Bézier曲线的定义为

(7)

其中，Pi=(xi，yi，zi)为控制点，wi是标量，为权因子，u为参变量，且u∈[0，1].P0P1P2构成了曲线的控制多边形，如图3所示.

图3 二次有理Bézier曲线Fig.3 Quadric rational Bézier curve

本文采用标准型的二次有理Bézier曲线，取w0=w1=1，w1>0，标准型的二次有理Bézier曲线为式(8)所示

(8)

由图3可知，当给定首末端点P0和P2，以及这两点处的切线方向T0和T2时，可以很容易地通过直线[P0，T0]和[P2，T2]的交点得到P1，再给定一点P就能唯一确定该曲线，因而确定了w1.

将所求的曲线看成是由点P0，P1和P2确定的抛物线的投影，P1为投影中心.如图3所示，将直线段[P0，P2]投影到要求的曲线上，则点P和Q为相对应的投影点.令w1=0，得到直线段L(u)=[P0，P2]，即

(9)

L(u)是点P0和P2的凸组合，因此|P0Q|和|QP2|的比值为u2：(1-u)2，从而推出

(10)

(11)

将u和P带入式(8)，得到w1，从而得到所求曲线.

(12)

假定型值点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)；

2)通过对sk求平均值得到逼近曲线的肩点坐标，即

(13)

则中间的权值为w=s/(1-s)；

3)根据P0，P1，P2和w即可确定型值点Qi与Qj之间的拟合曲线.

2.3.1 型值点处切矢的计算

由于数控加工程序中并不提供指令点处的切矢量，但可以通过型值点以及它周围的四个指令点来计算型值点处的切矢量.如图4所示，Qi为型值点，Qi-2、Qi-1、Qi+1、Qi+2为它周围的四个指令点，Qi的单位切矢表达式如下

图4 型值点处切矢的计算Fig.4 Tangent vector of reference point

(14)

Vi=(1-αi)qi+αiqi+1

(15)

qi=Qi-Qi-1

(16)

(17)

由于计算公式中未使用参数u，所求的单位切矢只能被看作切矢的方向.将下式带入式(14)～(17)，即可得到T0，T1和Tn-1，Tn.

q0=2q1-q2

(18)

q-1=2q0-q1

(19)

qn+2=2qn+1-qn

(20)

qn+1=2qn-qn-1

(21)

2.3.2 加工精度控制

虽然曲线段经过型值点Qi与Qj，但不能保证它到Qi与Qj之间的所有指令点的距离都满足最大拟合误差.因此，将Qi与Qj之间的所有指令点投影到拟合的曲线上，检查它们的偏差是否都在允许的误差范围内，即小于等于系统设定的最大轮廓误差δc，如果满足误差要求，进行下一段的曲线拟合；否则，将偏差最大的指令点设置为新的型值点，用新的型值点重新进行曲线拟合，再次检验拟合精度，重复这个过程，直至满足误差要求.

2.4 椭圆弧的识别

由于标准型的二次有理Bézier曲线只有一个权因子，其表达能力比非有理形式的Bézier曲线强，可以表达多种曲线，如图5所示，当01时，曲线为双曲线.

图5 不同w1定义的曲线Fig.5 Curves defined by w1

假设Ci(u)为加工路径中某段的标准型二次有理Bézier曲线，如果其权值，0

2.5 几何形式的转换

如果Ci(u)为加工路径中某段的标准型二次有理Bézier曲线，且其对应于一个椭圆弧段，如图6所示.

根据Ci(u)可得到该椭圆弧段的起点坐标p_start和终点坐标p_end.由Lee[16]的推导可知，该椭圆弧段的中心坐标为

图6 椭圆弧段几何信息的计算Fig.6 Calculation of geometric information of elliptical arcs

P1+ε(S+T)

(22)

S=P0-P1

(23)

T=P2-P1

(24)

(25)

(26)

其中，P0、P1和P2为Ci(u)的控制点，w1为权值.

该椭圆弧段的长、短半径分别为

(27)

(28)

假设λ2>λ1>0，且是以下二次方程的根

2δλ2-(kη+4β)λ+2(k-1)=0

(29)

(30)

(31)

β=S·T

(32)

该椭圆弧段长轴上的两个点为

q1=P1+(ε+r1x0)S+(ε+r1y0)T

(33)

q2=P1+(ε-r1x0)S+(ε-r1y0)T

(34)

根据q1、q2，可以得到长轴的斜率kl，和长轴与x轴正半轴的夹角d_kl；以长轴为x′轴，短轴为y′轴建立局部坐标系，利用q1、q2、中心坐标、起点坐标和终点坐标，求得相对于局部坐标系x′正半轴的起始角度d_start、终点角度d_end和椭圆弧的方向turn.

当连续微小线段加工区域完成椭圆弧的识别和几何形式转换后，会得到一个椭圆弧段数组Ellipse_Arcs[ ]，其中椭圆弧段的数据结构如下：

structEllipse_Arc_INF//椭圆弧段几何信息

{

doublecenter[2]； //中心坐标

doubler_l； //长半径

doubler_s； //短半径

doublep_start[2]； //起点坐标

doublep_end[2]； //终点坐标

doubled_start； //起点角度

doubled_end； //终点角度

doubleturn； //椭圆弧的方向

intkl； //长轴的斜率

doubled_kl； //长轴与x正半轴的夹角

}；

2.6 椭圆弧的合并

为了提高椭圆弧加工的效率和精度，将连续微小线段加工区域内的椭圆弧段合并，如果相邻的椭圆弧段具有相同的中心center、长半径r_l、短半径r_s、方向turn、长轴与x正半轴的夹角d_kl，并且上一椭圆弧段的终点坐标p_endi-1与下一椭圆弧段的起点坐标p_starti相等，则它们属于同一椭圆弧上连续的椭圆弧段，通过修改上一椭圆弧段的终点角度d_endi-1和终点坐标p_endi-1将它们合并为一个椭圆弧段.重复上述步骤，直至椭圆弧段不能再合并.

2.7 椭圆弧的插补

2.7.1 椭圆弧长计算

图7 区域及变量划分Fig.7 Region and variable partitions

(35)

根据高斯勒让德求积公式[17]，可将上式转换为

(36)

(37)

其它象限x或y为变量的情况，其弧长的推导方法相同.

2.7.2 弧长与变量的关系

2)根据公式(36)计算y0、y1、y2和y3到椭圆弧的起点y0的弧长，s0、s1、s2、s3.

3)将(s0，y0)、(s1，y1)、(s2，y2)和(s3，y3)代入如下三次多项式中进行插值，求得多项式的系数，即得到弧长与变量的关系.

y=a0+a1s+a2s2+a3s3，s∈[s0，s3]

(38)

其它象限x或y为变量的情况，其弧长与变量关系的推导方法相同.

2.7.3 椭圆弧实时插补

为了保证算法的实时性，对上述过程采用预处理方式.当加工程序读入数控系统后，对整个G01段进行处理，提前将加工路径中所有椭圆弧段的弧长与其参数关系进行计算并保存下来.在预处理结束后，再进行加工，直接利用计算好的弧长与参数关系实现实时插补.即在插补过程中，只需要7次浮点运算和5次加法运算，就能计算出第i个插补周期插补点坐标(xi，yi)，其具体步骤如下所示.

图8 椭圆弧平滑压缩算法流程图Fig.8 Flow chart of the proposed algorithm

假定插补周期为T，P0(xi-1，yi-1)和si-1分别为第i-1个插补周期的插补点坐标与弧长，第i个插补周期的进给速度为vi-1，则弧长增量为Δs=vi-1T，第i个插补周期的弧长为si=si-1+Δs，将si代入弧长与变量的关系中，得到对应变量xi或yi的值，再将其代入椭圆方程中，从而求得第i个插补周期的插补点坐标P(xi，yi).

2.8 算法流程

整个椭圆弧平滑压缩插补算法采用预处理方式，其流程如图8所示.

首先，读入一系列指令点，通过连续微小线段判断条件，识别出连续加工区域.在连续区域内，根据指令点的曲率，寻找型值点.然后，将相邻型值点间折线加工路径拟合为平滑的二次有理Bézier曲线加工路径，并检验每个指令点的拟合误差，将误差最大的点添加为新的型值点，重新进行拟合；如果都满足误差要求，从拟合曲线中识别出椭圆弧，并将其转换为几何形式；最后，将转换后的椭圆弧段与相邻的椭圆弧段合并.如果还有未拟合的型值点，则选取下一个型值点继续拟合，否则，计算椭圆弧段的弧长与变量的关系，然后判断是否还有未处理的指令点，如果有，重复上述操作直至最后一个指令点.

3 试验验证

为了验证该椭圆弧平滑压缩插补算法的正确性和有效性，本文使用了由沈阳高精数控智能技术股份有限公司自主研发的开放式四轴可控三轴联动数控系统GJ301，如图9所示，该系统采用工业级PC硬件平台，以及RTLinux实时操作系统.在验证过程中，为了与本文所提出的算法形成对比，将通过在小线段拐角处插入五次样条曲线的插补方法[11]称为传统算法.

图9 GJ301数控系统Fig.9 GJ301 NC System

设定系统的插补周期为2ms，x、y轴的最大加速度为A= 2000mm/s2，进给速度为F= 3000mm/min，双弓高误差为δmax= 0.01mm，最大曲率差值为δf= 0.0005mm-1，轮廓误差为δc= 0.01mm，试验中输入的初始加工路径如图10所示，该加工路径是按照0.001mm的精度将三个顺序相连的椭圆弧离散成211段小线段所生成的.

图10 加工路径Fig.10 Machining trajectory

由图11和图12的对比可以看出，传统算法在小线段的每个拐角处插入过渡曲线，平滑加工路径，提高加工速度，而本文的算法可以识别出所有的椭圆弧段并以三条曲线表示，加工路径的压缩比达到了70.3：1，同时具有较好的光滑性.

图11 传统算法拐角插入过渡曲线Fig．11 Inserttransitioncurveatthecorner图12 本文算法椭圆弧识别Fig．12 Recognizeellipsearc

加工过程中的速度曲线如图13、图14所示，传统算法和本文算法的速度曲线形状相同，但传统算法存在拐角，且加工时间分别为3.6s和3.464s.将匀速段放大后发现，传统算法存在频繁的波动，而本文算法的波动不明显.

图13 传统算法加工的速度曲线Fig.13 Speed curve of traditional algorithm

图14 本文算法加工的速度曲线Fig.14 Speed curve of proposed algorithm

计算轮廓误差，即由相邻插补点构成的小线段到原始(CAD)椭圆弧的最远距离，如图15和图16所示，传统算法和本文算法的轮廓误差均小于系统设置的最大值，但采用传统算法加工时，插补点位于小线段或过渡曲线上，产生的最大误差值个数较多，且平均误差值为0.0041mm，而本文算法生成的拟合曲线能够精确地表示椭圆弧，并且与原始(CAD)椭圆弧重合，使得插补点位于原始(CAD)椭圆弧上，其平均误差仅为5.6773x10-5mm，精度更高.对于系统中设置的双弓高误差δmax= 0.01mm和轮廓误差δc= 0.01mm，仅用于连续微小线段加工区域的识别与曲线拟合，虽然CAM编程相对于原始(CAD)椭圆弧存在0.001mm的误差，且拟合椭圆弧相对于CAM编程存在0.01mm的拟合误差，但这些误差主要是折线到曲线的最大距离，又因为拟合椭圆弧与原始(CAD)椭圆弧重合，且插补点位于拟合椭圆弧，而不是折线，所以本文算法能够在进给速度为3000mm/min的条件下，达到E-5mm量级的加工精度.

图15 传统算法加工的轮廓误差Fig.15 Contour error of traditional algorithm

图16 本文算法加工的轮廓误差Fig.16 Contour error of proposed algorithm

4 结 语

1)本文提出的椭圆弧平滑压缩插补算法采用预处理方式，不仅能够识别出数控加工程序中的椭圆弧，而且还能以椭圆弧的几何表示形式进行插补.

2)与传统的样条插补方法相比，提出的椭圆弧平滑压缩插补算法能够精确地识别并表示加工轨迹中的椭圆弧，提高了加工路径的压缩量和椭圆弧的加工质量，又因为本文用几何形式的椭圆弧进行插补，能够精确地计算出弧长对应的插补参数，降低了计算的复杂度和加工中速度的频繁波动，提高了加工效率.

3)本文提出的椭圆弧平滑压缩插补算法主要适用于在离散小线段形式下，椭圆弧部分的高速高精加工.

[1] Su Ai-lin，Zhang Chun-lei，Zhao Xue-feng，et al.A new method of ellipse interpolation-mapping method [J].Journal of Sichuan Union University，1999，3(1)：110-113.

[2] Hu Chi-bing，Cao Liang-wei，Yan Lin.An improved algorithm of elliptic interpolation [J].Journal of Lanzhou University of Technology，2007，33(6)：33-35.

[3] Yu Guang-wei，Zhao Yu-gang，Zhang Jian.Study of algorithm and path simulation by equal-error method to approach elliptical arc with beeline [J].Journal of Shandong University of Technology，2010，24(1)：90-92.

[4] Hou Zhi-qun.Setting-up and practice of equal segmental arc of ellipse curve points [J].Surveying and Mapping of Geology and Mineral Resources，2001，1(4)：14-16.

[5] Fukushima T.Precise and fast computation of a general incomplete elliptic integral of second kind by half and double argument transformations [J].Journal of Computational and Applied Mathematics，2011，235(14)：4140-4148.

[6] Boyd J P.Numerical，perturbative and chebyshev inversion of the incomplete elliptic integral of the second kind [J].Applied Mathematics and Computation，2012，218(13)：7005-7013.

[7] Fukushima T.Numerical inversion of a geVneral incomplete elliptic integral [J].Journal of Computational and Applied Mathematics，2013，237(1)：43-61.

[8] Gai Rong-li，Wang Yun-sen，Sun Yi-lan，et al.Survey on spline curve interpolation methods [J].Journal of Chinese Computer System，2012，33(12)：2744-2749.

[9] Tian Jun-feng，Lin Hu，Yao Zhuang，et al.Study on S-shape curve acceleration and deceleration control fast planning on CNC system [J].Journal of Chinese Computer System，2013，34(1)：168-172.

[10] Sun Shu-jie，Lin Hu，Zheng Liao-mo，et al.Algorithm of adaptive NURBS interpolation points calculation [J].Journal of Chinese Computer System，2014，35(4)：895-899.

[11] Zhang Xiao-hui，Yu Dong，Yang Dong-sheng，et al.Corner curve transition interpolation algorithm for high speed machining of micro-line segment [J].Journal of Mechanical Engineering，2010，46(19)：183-191.

[12] Zhang L B，You Y P，He J，et al.The transition algorithm based on parametric spline curve for high-speed machining of continuous short line segments [J].International Journal of Advanced Manufacturing Technology，2011，52(1)：245-254.

[13] Zhang Xiao-hui，Yu Dong，Hong Hai-tao，et al.Research on smooth compression interpolation algorithm in CNC machining [J].Journal of Mechanical Engineering，2011，47(5)：156-169.

[14] Zhang M，Yan W，Yuan C M，et al.Curve fitting and optimal interpolation on CNC machines based on quadratic B-splines [J].Science China Information Sciences，2011，54(7)：1407-1418.

[15] Yang Z Y，Shen L Y，Yuan C M，et al.Curve fitting and optimal interpolation for CNC machining under confined error using quadratic B-splines [J].Computer-Aided Design，2015，66 (7)：62-72.

[16] Farin G E.Geometric modeling：algorithms and new trends [M].Society for Industrial and Applied Mathematics(SIAM)，1987.

[17] Yin Chen-chen，Wu Wen-jiang.A new ellipse interpolation method based on arc length [J].Modular Machine Tool & Automatic Manufacturing Technique，2012，4(1)：1-4.

附中文参考文献：

[1] 苏爱林，张春雷，赵雪峰，等.椭圆插补的新方法-映射法[J].四川联合大学学报，1999，3(1)：110-113.

[2] 胡赤兵，操良伟，闫 琳.一种椭圆插补的改进算法[J].兰州理工大学学报，2007，33(6)：33-35.

[3] 于光伟，赵玉刚，张 键.椭圆弧等误差直线拟合的算法研究与轨迹仿真[J].山东理工大学学报，2010，24(1)：90-92.

[4] 侯至群.等弧长椭圆曲线点的放样与实践[J].地矿测绘，2001，1(4)：14-16.

[8] 盖荣丽，王允森，孙一兰，等.样条曲线插补方法综述[J].小型微型计算机系统，2012，33(12)：2744-2749.

[9] 田军锋，林 浒，姚 壮，等.数控系统S型曲线加减速快速规划研究[J].小型微型计算机系统，2013，34(1)：168-172.

[10] 孙树杰，林 浒，郑飂默，等.计算复杂度自适应的NURBS曲线插补算法[J].小型微型计算机系统，2014，35(4)：895-899.

[11] 张晓辉，于 东，杨东升，等.面向微线段高速加工的拐角曲线过渡插补算法[J].机械工程学报，2010，46(19)：183-191.

[13] 张晓辉，于 东，洪海涛，等.数控加工中的平滑压缩插补算法研究[J].机械工程学报，2011，47(5)：156-169.

[17] 殷晨晨，吴文江.基于弧长的椭圆插补新算法[J].组合机床与自动化加工技术，2012，4(1)：1-4.