在数控加工中，刀具轨迹通常用G01微小线段表示，导致数控系统前瞻距离短，最终限制加工速度；同时，微小线段轨迹缺乏连续性，导致进给速度、加速度波动以及机床振动[1-2]．使用样条曲线拟合G01微小线段，用更少的数据表示原始轨迹，能实现轨迹压缩并提升轨迹连续性．B样条拟合计算复杂度与拟合的微小线段数成正比，为减小单次拟合计算量，整个刀具轨迹通常被划分为若干区间分段拟合．文献[3]使用文献[4]估算曲率的方法，对连续微小线段每个转角处的曲率进行估算，并由此判断连续微小线段的偏转情况，将偏转过大的点作为特征点划分压缩区域的边界．文献[5]根据微小线段间夹角和用相邻三点构成的圆弧与对应直线段之间的弓高误差作为连续短线段(CSB)的判断条件，并对连续短线段区间进行压缩．为了进一步提高压缩效率，通常会预处理待压缩刀位点，选取部分点作为初始压缩节点．现有节点选取方法对选取初始压缩节点具有借鉴意义．文献[6]使用一种迭代节点插入(IKI)方法，从少量节点开始，找到逼近误差高于给定允差的节点段，添加新节点．文献[7-8]将节点选取转化为凸优化问题，使p阶B样条的p-1阶导数阶跃最小，该优化方法可以同时计算节点的数目和位置，并在节点数较少时获得较低的近似误差，但计算成本更高．文献[9]使用一种基于支持向量机的机器学习方法进行节点选取，该方法的性能在很大程度上依赖于训练数据集，从而限制了其普适性．文献[10]通过混合导数、曲率、转角三种几何特征来选择节点．文献[11]基于高阶导数的特征函数积分快速无迭代地选取节点，同时建立了拟合误差与节点数的关系．文献[12]选取首末点、曲率极值点、曲率过零点中的弦高最大点作为拟合的主导点．文献[13]通过对输入数据的二阶导数进行分段线性逼近，选取三次B样条逼近的节点．文献[14]使用曲率峰值作为初始节点，通过迭代添加新节点或调整现有节点以减少曲率偏差．文献[15]使用输入数据集的曲率和弧长来选择初始节点，并通过插值该点集来构造B样条．为了减小单次B样条压缩计算量，对刀具轨迹进行分段压缩，有研究者提出分段后相邻轨迹段衔接点处连续性低的问题[16]．文献[17]提出了一种带有前瞻功能的实时NURBS插值，该方法基于连续短线段生成G2连续的刀具轨迹，但其只处理了少量的G01连续点集及NURBS曲线段，没有解决曲线衔接处的连续性问题．文献[16]在相邻B样条曲线间使用另一条B样条进行G2连续桥接，但其在构造桥接曲线时，须要对控制点位置进行大量分类讨论．文献[18-19]使用局部插值来构造G2连续的贝塞尔(Bezier)曲线段，并开发了一种局部修调的插值方案来处理曲线形状缺陷，但该方法会构造大量贝塞尔曲线，导致存储和计算负载．文献[20]用对称的B样条曲线桥接两个连续的直线段，文献[21]利用中间两控制点重合的三次Bezier曲线桥接相邻样条，但该方法仅适合桥接处曲率为0的情况，应用比较受局限．上述研究对于初始压缩节点的选取，一方面大多是在三维曲线上进行，或选取过程有迭代，效率低；另一方面，多数并未考虑相邻压缩轨迹间的G2连续问题，或改善曲线连续性时须要分情况讨论．为此，本研究推导出一种快速选取初始压缩节点，桥接改善整个刀具轨迹连续性的压缩方法．1 基于G2连续桥接的轨迹压缩算法如图1所示，本算法主要分为B样条压缩与G2连续桥接两个阶段．对于输入的G01微小线段，估算离散点曲率ki(i=1,2,…,q)，q为离散点个数，比较ki与曲率阈值kT，选取超过曲率阈值的点作为区间划分的特征断点．对于划分的待压缩区间，基于离散曲率积分选出合适的初始压缩节点，利用文献[22]中的方法，得到多段B样条曲线轨迹．10.13245/j.hust.230608.F001图1微小线段G2连续压缩算法流程图分段压缩后在相邻B样条曲线间进行桥接，针对桥接曲线的边界约束条件引入微分因子μ，保证其与相邻B样条曲线间G2连续，μ通过误差和保型约束求解．该方法无须分类讨论，求解简单．2 刀具轨迹压缩方法在此提出基于三次B样条刀具轨迹压缩方法，包括基于离散曲率选取区间划分断点，基于离散曲率积分选取初始压缩节点，最终压缩得到多段B样条曲线轨迹．2.1　基于离散曲率的区间划分方法区间划分的关键在于如何选取划分断点．通常曲率越大的轨迹段，拟合过后越容易形成更大的弦高误差；因此，本研究先估算各数据点的离散曲率模长，通过离散曲率确定曲率阈值，进而设计断点选取要求．具体分为以下三步：a. 通过三点拟合圆的方法[22]估算离散数据点的曲率的模ki；b. 根据ki计算曲率阈值kT=top({ki}1q,ρ⋅q)；c. 根据kT确定区间划分断点．具体算法内容如下．a. 计算离散曲率大小设Qi-1，Qi，Qi+1为离散点集中任意相邻三个点，使用三点拟合圆估算曲率的方法[22]，得到Qi的离散曲率大小为ki=2Qi-1Qi×QiQi+1Qi-1Qi⋅QiQi+1⋅Qi-1Qi+1．(1)由于曲率是估算且离散的，该计算方式可能导致曲率局部抖动频繁，因此可以通过滤波去除曲率噪点影响．b. 确定曲率阈值曲率阈值的选取直接影响到区间划分数量，若划分区间数过少，则单次拟合计算量仍较大，影响计算效率；若划分区间数过多，则区间划分与各区间拟合的总计算量也会成为负担．因此，曲率阈值由用户根据轨迹曲率大小分布特点设定，尽量使得较大曲率点均被视为断点，避免拟合后产生较大弦高误差．令曲率阈值为kT=top({ki}1q,ρq)，(2)式中：top({ki}1q,ρq)表示集合{ki}1q中第ρq大的值；ρ为用户设定的比例参数，显然ρ越小，区间划分数越少，反之越多．本研究案例取ρ=10%．c. 选取区间划分断点断点的选取应满足以下两个条件：a. 断点处曲率的模大于曲率阈值，即kikT．b. 断点处曲率的模为局部极大值，即ki=max({kj}i-Ri+R)，其中R为局部极值范围半径参数，在此选取R=2．2.2　基于离散曲率积分的初始压缩节点选取方法刀具轨迹通过上述方法划分区间后，须要对各子区间进行压缩．一组合适的初始压缩节点可以减少后续拟合迭代次数，提高压缩效率．显然，希望初始压缩节点的稀疏程度与曲率大小成正比，即一段轨迹所需的节点数与其曲率模长的累计值正相关．设第i段直线段的弦长为Δsi=QiQi-1，为了简化积分运算，以弦长代弧长，定义离散曲率积分G(j)为ki关于刀位点弧长的积分，n为刀位点数目，G(j)=∫-∞jkds≈∑i=1j(ki-1+ki2Δsi) (1≤j≤n)．(3)当曲率为定值，即轨迹为圆弧时，显然希望选取的初始压缩节点等弧长均匀分布，即节点数f(j)与G(j)值为线性关系．显然，节点数满足f(j)≥1且为整数，因此取线性关系常数项为h=1．又设该线性关系一阶因子为d且d≥0，建立前j个刀位点所需的初始压缩节点数f(j)与G(j)的线性关系为f(j)=dG(j)+1．(4)假定目标压缩比为P，且令所选取的初始压缩节点数恰好满足目标压缩比，其中压缩比为刀位点数目与压缩轨迹节点数的比值，则有P=n/f(n)=n/(dG(n)+1)．(5)由式(5)可解得d=n/[PG(n)],d≥0．(6)将式(3)、(6)代入式(4)，得f(j)=(nG(j))/[PG(n)]+1=   nP∑i=1j(ki-1+ki)Δsi  /∑i=1n(ki-1+ki)Δsi+1                           (1≤j≤n)．(7)由式(7)，计算压缩所有刀位点所需的初始压缩节点数，再对离散曲率积分均匀采样，将采样点映射回原轨迹，得到初始压缩节点．以图2刀具原始轨迹为例展示算法过程，设定目标压缩比为5．10.13245/j.hust.230608.F002图2刀具原始轨迹a. 刀具原始轨迹如图2所示．b. 计算初始压缩节点在离散曲率积分上的分布，计算过程如图3所示，其中s为刀位点弧长．10.13245/j.hust.230608.F003图3初始压缩节点选取的计算过程c. 将初始压缩节点映射回原轨迹，其位置分布如图4所示．10.13245/j.hust.230608.F004图4原轨迹中初始压缩节点位置分布由图4可见：初始压缩节点在曲率较大处分布密集，在曲率较小处分布稀疏．2.3　基于B样条曲线的刀具轨迹压缩方法区间划分完成后，使用三次B样条曲线对各子区间进行压缩[23]，算法流程分为以下五步：a. 使用上述方法选取初始压缩节点序列{Qa,Qa+1,…，Qb}，新节点在迭代过程中添加，以此来尽可能减少节点数量，提高压缩比；b. B样条曲线的节点矢量使用文献[22]中的规范积累弦长参数化方法确定；c. B样条曲线的首末端矢量使用Boehm法[24]进行构造；d. 反算控制点矩阵方程[23]求解B样条控制点；e. 迭代压缩过程中，计算当前数据点与拟合B样条间的最小距离作为误差与给定允差相比较，判断是否添加当前点到新的压缩序列中，迭代进行步骤a至e，直至所有点满足误差要求．同理完成对所有子区间的压缩，得到由多段B样条曲线表示的刀具轨迹．3 G2连续桥接方法刀具轨迹压缩完成后，根据B样条曲线性质，各子区间内部G2连续，但相邻区间衔接点处连续性还未考虑．本研究使用五次六控制点的贝塞尔曲线，截断相邻两B样条曲线轨迹，使得截断点处达到G2连续，最终减小单轴加速度波动．该算法流程分为以下五步：a. 截断参与桥接的B样条曲线，计算其在截断点处的一阶、二阶导矢；b. 通过保证桥接曲线在截断点处的C2连续，确定桥接曲线六控制点的唯一解；c. 将截断点处的一阶、二阶导矢乘以微分因子μ(μ≠0)，增加桥接曲线自由度，同时保证截断点处仍然满足G2连续；d. 基于桥接曲线与压缩后的B样条曲线间的误差，构造目标优化函数；e. 附加保型约束，优化目标函数，求解μ(μ≠0)，最终求得桥接曲线．3.1　桥接曲线截断点选取方法为了简化计算，取相邻B样条曲线衔接点前后最近的两个节点作为桥接曲线截断点，如图5所示．10.13245/j.hust.230608.F005图5桥接前后曲线示意图C1,C2为桥接前相邻两条B样条曲线，C为桥接曲线，Q为两B样条衔接点，Q0,Q1为B样条曲线上离衔接点最近的两节点，即为桥接曲线首末端点，也即首末控制点b0,b5．3.2　基于微分因子μ的边界G2连续约束方法设t0,t1,d0,d1分别为原轨迹节点Q0,Q1处的一阶、二阶导矢，P0',P1',P0″,P1″分别为桥接曲线首末端点一阶、二阶导矢．根据五次贝塞尔曲线性质，可得P0',P1',P0″,P1″与桥接曲线控制点的关系为：P0'=5(b1-b0);P1'=5(b5-b4);P0″=20(b2-2b1+b0);P1″=20(b5-2b4+b3), (8)式中bi(i=0,1,…,5)为控制点．为保证桥接曲线端点处C2连续，应当满足以下条件：b0=Q0,  P0'=t0,  P0″=d0;b5=Q1,  P1'=t1,  P1″=d1. (9)由式(8)和(9)求得桥接曲线控制点表达式为b0b1b2b3b4b5=Q0Q0+t05Q0+d020+2t05Q1+d120-2t15Q1-t15Q1. (10)根据式(10)，六个控制点由首末端点位置、一阶、二阶导矢唯一确定，无法再施加误差与保型约束，因此引入微分因子μ(μ≠0)，增加曲线自由度同时保证桥接曲线与相邻B样条曲线间G2连续．为此证明以下性质：在对上述曲线截断点处的C2连续边界约束条件施加微分因子μ(μ≠0)后，截断点处仍满足G2连续．设一截断点处的一阶、二阶导矢P',P″，则该点曲率为k=P'×P″/P'3．(11)若一阶导矢乘以微分因子μ(μ≠0)，其方向不变，为了保证截断点前后曲率不变，根据式(11)，则有下式成立：P˜'=μP';P˜″=μ2P″;μ≠0, (12)式中P˜',P˜″分别为新的一阶、二阶导矢．微分因子的引入增加了桥接曲线自由度，给后续优化留有空间；同时，在更新边界条件并修调桥接曲线控制点b1,b2,b3,b4后，曲线端点处一阶导矢方向保持不变，曲率保持不变，即桥接曲线端点处达到G2连续．根据式(10)和(12)，修调后桥接曲线控制点b˜1,b˜2,b˜3,b˜4为：b˜1=Q0+μ1t05;b˜2=Q0+μ12d020+2μ1t05;b˜3=Q1+μ22d120-2μ2t15;b˜4=Q1-μ2t15, (13)式中μ1,μ2(μ1,μ2≠0)分别为两截断点处的微分因子．3.3　基于误差和保型约束的微分因子求解方法已经得到桥接曲线控制点关于微分因子的表达式，在此建立关于桥接误差的目标优化函数，从而求解μ1,μ2(μ1,μ2≠0)．先定义桥接误差为桥接曲线与原B样条曲线间的误差．由于桥接样条较短，过渡较为平缓，作如下近似简化：取桥接曲线C(0.5)处的数据点与衔接点Q之间的距离Dmax作为桥接误差．如图5所示，C1,C2为相邻两B样条曲线，C为桥接曲线，Q0,Q1为截断点，Q为衔接点．由五次贝塞尔曲线性质，可得C(0.5)=b0+5b˜1+10b˜2+10b˜3+5b˜4+b532, (14)则误差表示为Dmax=C(0.5)-Q=(b0+5b˜1+10b˜2+10b˜3+5b˜4+b5)/32-Q．(15)为了便于求解，构造目标优优化函数F为Dmax与允差δ的差值，并整理为如下形式F=32C(0.5)-32Q2-1 024δ2=(b0+5b˜1+10b˜2+10b˜3+5b˜4+b5-32Q)2-1 024δ2. (16)将式(13)代入式(16)，整理可得F=[d0μ12/2+d1μ22/2+5t0μ1-5t1μ2+(16Q0+16Q1-32Q)]2-1 024δ2, (17)式中Q0,Q1,t0,t1,d0,d1,δ均为已知量．再增加保型约束：b1b˜2Q0Q1≥0;b˜2b˜3Q0Q1≥0;b˜3b4Q0Q1≥0. (18)代入式(13)，整理可得：(d0μ12/20+t0μ1/5)(Q1-Q0)≥0;    (-d0μ12/20+d1μ22/20-2t0μ1/5-2t1μ2/5+(Q1-Q0))(Q1-Q0)≥0;(-d1μ22/20+t1μ2/5)(Q1-Q0)≥0. (19)最终，该问题转换为在式(19)的约束下，式(17)的最优化问题，使用Matlab优化工具求得微分因子μ1,μ2(μ1,μ2≠0)，进而求得桥接曲线方程．4 仿真验证4.1　轨迹压缩算法验证结合实际加工要求和轨迹压缩算法的性能，给定压缩允差0.01 mm．为了定量评估压缩曲线，建立如下评价标准．假设原始刀位点集为{Q1，Q2,…,Qi,…,Qn}，n为刀位点数目，计算当前刀位点Qi到压缩曲线Cj(u)的最近距离minu∈[0,1]disi(u)=Qi-Cj(u)作为压缩误差，应小于给定压缩允差值，j为当前压缩B样条曲线序号．压缩效率使用压缩算法耗时进行评价，通常数控系统插补周期为1 ms，因此认为平均每个刀位点压缩耗时小于1 ms即为高效．压缩比为原始轨迹刀位点数与压缩后B样条曲线节点数的比值，压缩后减少了轨迹程序信息量，增大了数控系统前瞻距离．基于以上评价标准和压缩允差，对给定S型轨迹进行压缩如图6所示，记录压缩结果如下：平均压缩误差为1.875 77 μm，最大压缩误差为9.938 67 μm，压缩比为43.24，压缩耗时为0.764 s．对给定轨迹进行压缩的结果表明压缩误差满足给定允差要求；同时，对于该4 670个刀位点的轨迹压缩比达到43.24，减少了程序信息量；压缩耗时仅0.764 s，平均单个刀位点耗时仅0.164 ms，说明对于该轨迹压缩算法效率高．10.13245/j.hust.230608.F006图6G01轨迹压缩结果4.2　G2连续桥接算法验证对桥接前后的轨迹进行模拟插补实验，对比桥接前后刀具轨迹曲率，重点观察桥接前B样条曲线衔接点与桥接后截断点处的曲率变化，以验证G2连续桥接算法的正确性．给定刀尖进给速度为1 000 mm/min，观察单轴加速度波动，以验证桥接算法有效性．对图6刀具轨迹衔接点进行桥接，图7(a)所示为其中相邻两B样条轨迹曲线及其桥接轨迹图，图7(b)为桥接轨迹放大图，可见桥接曲线过渡平滑，无打圈打折现象．10.13245/j.hust.230608.F007图7桥接算法实施后的刀具轨迹对比桥接前后刀具轨迹曲率如图8所示，图中：k为离散曲率，ni为离散点数．在桥接前衔接点处曲率误差大小为10-3数量级，而桥接算法实施后，两截断点处的曲率误差大小数量级仅为10-6和10-15，考虑到计算误差，可以认为在桥接前衔接点处曲率不连续，桥接后截断点处曲率连续．10.13245/j.hust.230608.F008图8桥接前后刀具轨迹曲率图以z轴加速度为例，绘制其单轴加速度如图9所示．桥接前刀具轨迹衔接点处加速度均出现不同10.13245/j.hust.230608.F009图9相邻两段B样条轨迹桥接前后单轴加速度程度的波动，桥接后桥接曲线截断点附近加速度过渡平缓，波动均减小；桥接前加速度波动幅值为1.1 mm/s2，桥接后加速度波动幅值为0.7 mm/s2，减小约36%．因此，G2连续桥接算法减小了单轴加速度波动．如图10所示，为图6刀具轨迹在桥接算法实施前后的加速度对比图，显然，其中B样条曲线衔接点处加速度波动明显减小．10.13245/j.hust.230608.F010图10整条刀具轨迹桥接前后单轴加速度对比5 结语本研究提出了一种压缩与桥接相结合的刀具轨迹平滑方法，一方面利用轨迹的离散曲率积分选取初始压缩节点，考虑了轨迹的曲率特点，且避免迭代求解；另一方面引入带有微分因子μ的边界约束条件，提出了一种计算简洁的G2连续桥接方法．在此基础上，进行了模拟插补实验，得到以下结论．a. 提出的压缩算法极大地减少了原始G01程序信息量，原始程序节点数相比压缩后的节点数即压缩比达到43.24．该算法迭代快、效率高，对于具有4 670个刀位点的轨迹，压缩耗时仅0.764 s，单个刀位点平均耗时仅0.164 ms．b. 针对分段压缩轨迹衔接点处的连续性问题，提出的G2连续桥接算法改善了衔接点处的连续性，减小了单轴加速度波动，桥接后原先轨迹衔接点及桥接轨迹截断点处均未出现加速度波动．