目前，多自由度工业机器人已经广泛应用于自动化生产线、焊接制造及物料搬运等重复性较高的应用场景；同时，工业机器人也正逐步应用于航空航天、军工、汽车制造等领域代替人工制孔．但由于机器人本体刚性较弱，仅为传统数控机床的1/15[1]，在受到大钻削力的作用下，系统的刚性不足产生不稳定振动，从而产生较大的加工孔毛刺，影响制孔质量．而去除毛刺过程相对繁琐，同时去除毛刺容易造成安装孔部位材料疲劳失效，会带来安全隐患．因此，须要从机器人制孔毛刺产生机理入手，结合数据方法对机器人制孔毛刺高度进行预测，从而最终实现制孔毛刺的闭环控制．在钻削毛刺生成机理方面，众多研究者已对其做了大量的研究．文献[2]通过实验建立钻削轴向力经验公式，然后采用有限元方法计算单层板在理论钻削力下的变形量，进而预测层间毛刺生成初始位置，最后运用能量法建立叠层钻削层间毛刺高度理论模型；文献[3]从各部分所做的功开始，包括轴向力功、材料变形功和总功，同时考虑了声软化效应，最后根据能量守恒法计算毛刺高度；文献[4]建立了一种新的分析模型，用于描述延性金属加工中刀具几何形状、操作参数和工件材料性能的函数，以描述进口毛刺的形成．但以上研究大多基于传统机床，未考虑机器人制孔系统相关特性，同时单纯的数学解析建模无法涵盖所有影响毛刺生成因素，且一些参数的选择较为理想化，预测准确率难以满足要求．与此同时，以实时加工过程数据驱动的黑箱预测模型目前也已广泛应用于各种质量预测场景，该种预测方式主要以神经网络为基础．文献[5]以工艺参数与主轴电流信号为网络输入，以制孔出口毛刺高度为预测目标，从而实现毛刺高度的预测；文献[6]提出了一种基于信息融合的刀具剩余寿命在线预测方法，在线预测刀具剩余寿命；文献[7]基于改进CNN-BiLSTM模型，通过提取的切削力、振动和声音信号的空间和时序特征实现刀具磨损状态监测；文献[8]基于提取的振动信号特征及混合神经网络模型实现对机械设备健康度的预测．然而，基于数据驱动的黑箱模型也存在可解释性较差、泛化性不强等缺点．因此，在部分复杂系统状态研究领域，机理与数据融合的预测方法融合了白箱模型与黑箱模型二者优点，成为研究重点．文献[9]通过建立的螺栓松脱过程机理模型及基于高斯过程回归的螺栓松脱特性预测模型，对螺栓连接松脱进行预测，与传统模型相比具有更好的预测效果．文献[10]提出了一种基于知识和数据融合的设备故障诊断方法，对比传统分类模型准确率表现突出，并实现了机理知识与数据驱动相融合的设备故障领域图谱构建．文献[11]提出了数字孪生驱动下的弱刚性钻削毛刺控制方法，实现对钻削加工过程的实时优化调控．由于制孔毛刺的生成机理复杂，单纯通过机理模型进行预测准确率难以保证，因此本研究提出一种机理与数据融合的机器人制孔毛刺预测方法，通过对实时数据的采集，采用机器学习对原始物理模型进行动态更新优化，提高制孔毛刺预测的准确程度．本研究在建立的静态物理模型基础上引入动态更新模型，对物理模型结果进行更新，得到具备自我演化的融合预测模型，并在此基础上开展相应的机器人制孔实验进行验证．1 出口毛刺高度物理模型建立随着钻削过程中主轴的进给，钻头不断切削金属材料，当接近出口部位时，孔中心部位的下表面材料发生屈服，并由于材料断裂类型的不同产生不同种类的毛刺．图1是钻削后期的两种主要情形，由于材料底部加工硬化程度的不同，分为钻头刺穿工件材料与钻头未刺穿工件两种情形．10.13245/j.hust.220616.F001图1钻削后期毛刺生成的两种情形由于在切削力以及系统弱刚性的影响下，制孔系统极易产生颤振，从而改变钻头钻削的实际钻削厚度，并进一步影响该条件下钻削毛刺的生成机理；因此，基于钻削末期工件的振动模型，并在其基础上联立工件材料特性建立出口毛刺高度预测模型．由于未钻削部分工件阻尼较小，因此假设振动模型为无阻尼模型mx¨+kx=Fsinωt，(1)式中：m为未切削部分材料质量；x为系统振动位移；k为末端刚度；F为激振力；ω为振动频率；t为时间．基于式(1)，通过建立系统刚度、钻削力、振动模型，最终实现毛刺高度建模．1.1　机器人刚度建模机器人 D-H建模与D-H参数分别见图2和表1．基于D-H约定对机器人进行运动学建模并得到机器人雅可比矩阵，通过对机器人关节刚度进行辨识，拟合得到机器人的关节刚度矩阵，在已知机器人关节刚度的基础上，可以得知机器人在任意位姿下的操作刚度矩阵K(q)，将操作刚度矩阵分解为四个对称部分，即K(q)=KfdKfδKmdKmδ，(2)式中：Kfd为力-线位移刚度矩阵；Kfδ为力-角位移刚度矩阵；Kmd为转矩-线位移刚度矩阵；Kmδ为转矩-角位移刚度矩阵．参照文献[12]，得到机器人末端在所受力矢量方向上的刚度为kf=1/(eTKfd-1e)，(3)式中e为机器人末端所受力矢量的方向向量．根据力矢量方向与两投影方向的夹角可以得到机器人制孔系统中钻头进给方向的刚度k2与垂直进给方向的刚度k1．10.13245/j.hust.220616.F002图2机器人D-H建模10.13245/j.hust.220616.T001表1机器人D-H参数连杆di/mmai-1/mmαi-1/(°)θi/(°)154000020320-90-90308700041 020235-90050090062000-9001.2　钻削力建模钻削过程中钻削力的产生情况如图3所示，图3(b)为主切削刃微元局部放大示意图． 钻削力按照形成的部位可分为主切削刃、副切削刃及横刃这三部分切削分力的合力．10.13245/j.hust.220616.F003图3钻头受力情况其中主切削刃部分的切削力由工件被钻削时产生的弹性抗力和钻头与工件摩擦所产生的摩擦力组成．基于微元法建立单元体弹性抗力dFe关系式，并对其积分可得到主切削刃所受压力Fp，dFe=εhdEtanθdxdL;Fp=∫0L∫0fa2εhdExtanθcosθdxdL, (4)式中：ε为实验系数；hd为单位切削厚度；E为工件弹性模量；θ为流屑角；L为主切削刃长度；fa为刃带宽度．将弹性抗力分解为垂直于主切削刃所在平面的压力分力及沿平面的分力，并将上述空间力系建模投影到正交平面内，得到压力投影Fp'、摩擦力投影Ff、沿斜面分力投影Fn分别为Fp'=FpsinΦsinγ0;Ff=FpffsinΦcosγ0;Fn=FptanθsinΦcosγ0, (5)式中：Φ为投影角度余角；ff为动摩擦系数；γ0为单元钻削点处的前角．则主切削刃上的钻削力Fm=2hdε[(ff-tanθ)cosγ0-sinγ0]·∫0L∫0faExtanθdxcosθdL. (6)横刃所受的钻削力主要由横刃切入工件受力与横刃切除材料所受塑性变形抗力组成，第一部分力Fc1与工件硬度相关[13]，Fc1=0.002 2HBd2，(7)式中：HB为材料硬度；d为钻头直径．第二部分力Fc2=2∫0lσ1δfsinαdr，(8)式中：l为横刃长度；σ1为工件材料剪切强度；δf为横刃工作宽度；α为摩擦力与水平方向夹角．则横刃所受总钻削力Fc为Fc1与Fc2之和.各切削刃部分的切削力按照能量占比分析，主切削刃占总切削力的40%，横刃占57%，副切削刃占3%[14]，因此钻削过程中钻削力的合力为Fr=Fm+Fc+Fm40%+Fc57%3%2．(9)1.3　系统振动建模本研究定义在钻削过程末期存在的瞬时切削厚度降低的时期为不稳定期，将轴向力分解为挤压力与切削力，同时建立颤振模型与强迫振动模型加以叠加．由于在钻削过程中垂直于进给方向上受到动态钻削力，进给方向受到静态钻削力，因此可将钻削过程简化为图4所示二自由度问题．10.13245/j.hust.220616.F004图4简化振动模型当钻削颤振发生时，钻头每旋转一个周期T，切削厚度便会产生一个周期性的叠加，假设两个方向上的切削厚度的微变量分别为∆x和∆z，则钻削的瞬时切削厚度为Δh=Δxsin(90°-pp/2)+Δzcos(90°-pp/2)，(10)式中pp为钻头锋角．由于两个方向上的瞬态切削力与瞬时切削厚度成正比，因此两个方向上的瞬态切削力为：Fx=λ1Δh;Fz=λ2Δh. (11)为了突出刚度耦合的效应，削弱位移反馈带来的负阻尼效应，本研究在颤振模型中忽略阻尼的影响．结合上式建立颤振的解析模型：mx¨+k1x=-Fx;mz¨+k2z=-Fz;x=A1ept;z=A2ept. (12)将位移、加速度代入式(12)进行化简后得：(mp2+k11-λ11)A1+(k12-λ12)A2=0;(k22-λ21)A1+(mp2+k21-λ22)A2=0. (13)若方程组(13)存在非零解，其系数矩阵必等于0，根据求值公式得到p的解，鉴于自激振动的条件，p的解为虚数；同时根据稳态情况下的加速度信号求解A，最终求得x和z方向上的位移方程x1和z1：x1=x1+x2;z1=z1+z2, (14)式中x1和z1等号右侧为方程的两个解.强迫振动的激振力为挤压力，包括孔壁的挤压力及轴向的挤压力．在振动基本方程右端乘以一个频率为振动频率的正弦函数，并两边同时除以m，得x¨+2ζωnx˙+ωn2x=xstωn2sinωt;xst=F/k, (15)式中：ζ为阻尼比；ωn为系统自然频率；xst为系统在力F作用下的静变形．式(15)对应的齐次方程通解和特解分别为xtong=exp(-ξωnt)(c1cosωdt+c2sinωdt);xte=Asin(ωt-φ), (16)式中：ωd为阻尼频率；A为振幅；φ为相位角．将方程解带入运动方程，若方程在任意时间t下的值均相等，则方程中关于项sin ωt和cos ωt前的系数相等，从而可以求得振幅A和初始相位角φ0：A=xst/1-(ω/ωn)2+4ξ2(ω/ωn)2;φ0=arctan{2ξ(ω/ωn)/[1-(ω/ωn)2]}. (17)代入初始条件之后即可计算出x和z两个方向上的振动值x2和z2：x2=xtong+xte;z2=ztong+zte. (18)综上所述，整个系统在两个方向上的最终振动位移方程为：x=x1+x2;z=z1+z2. (19)1.4　毛刺高度建模钻削后期分为两种情况，通过钻削末期传感器信号变化平稳程度判断：第一种情况为钻头刺穿工件材料，假设未参与毛刺形成部分的材料最终均被切除，参与毛刺形成的部分保留下来并最终转化为毛刺，此时参与毛刺形成部分的厚度为最大未切削厚度h1；第二种情况为钻头未刺穿工件材料，此时刀具使未切削部分发生挤压变形，在边缘部分发生断裂，确定毛刺的高度就是确定产生初始断裂纹的位置．激振力F为挤压力的轴向分量，大小为F=fR2/(2h)-2R(kp+2Zn/h)(1-e(G/Z)t);wp=f/2, (20)式中：wp为单个切削刃刀具旋转一周切削的刀具材料厚度；h为未切削部分材料厚度；G为剪切弹性模量；Z为黏性系数．定义未切削部分材料刚度k=c(h/d)，(21)式中c为系数．联立强迫振动模型，得到稳定期最大未切削厚度h0=(F0sinωt)/[ρπr2z¨1+(c/d)z1]，(22)式中ρ为工件材料的密度．假设h1与h0之间的关系为h1=Ah02+Bh0+C．(23)在钻削后期的第一种情况，由于径向振动，钻头在孔中存在一个等效间隙δs，其对毛刺产生有一定影响，假设其等于径向振动的振幅，并假设在Kδs处产生裂纹，则毛刺高度为  H=khz/x{(r-Kδs)/tan(pp/2)+[h1-(r-Kδs)/tan(pp/2)]/sin(pp/2)}(1-ψ), (24)式中ψ为材料断面切削率．第二种情况下，毛刺的形成是由于塑性变形后的弹性断裂，假设毛刺形成的位置在塑形生成的一半位置，此时工件材料的突出高度为H1=keh1H．(25)假设破坏应变为εf，其与δs和H的关系为εf=Kδsεf2+2εf/(keh1)．(26)考虑材料属性断面收缩率，同时轴向振动的剧烈程度也会影响钻削状态，因此第二种情况的钻削毛刺高度为H=K|z|[δsεf2+2εf/(2keh1)](1-ψ)．(27)2 数据驱动模型建立由于机加工系统较为复杂，各个系统参数的变化均可能会对加工质量产生影响，因此当进行物理建模时，难以将所有加工参数考虑到模型当中，且可能会产生较大误差．机器人制孔系统在钻削过程中会产生大量的加工状态参数，这些参数会在一定程度上反映加工状态，为加工质量、刀具磨损等提供一定的参考价值．2.1　机器人制孔信息物理系统构建本研究设计的机器人制孔信息物理系统体系架构如图5所示，其中钻头相关信息包括钻头的材料、几何参数(直径、锋角、刃倾角等)；工件相关信息包括工件的材料、几何尺寸；加工参数包括机器人关节角、末端主轴转速及进给量．末端执行器工作状态的读取控制通过串口与上位机通信；机器人本体的工作状态读取及控制由Ethernet通信实现，最终实现制孔系统的控制及状态监控．在上述模型的基础上，建立基于数据驱动的动态更新模型．虽然建立了制孔毛刺预测的机理模型，但在实际钻削过程中可能会受到刀具磨损等其他加工环境因素的影响，因此须要采集实时加工数据对机理模型进行修正．10.13245/j.hust.220616.F005图5机器人制孔信息物理系统体系架构2.2　动态更新模型构建动态更新模型的工作流程如下：先基于静态物理模型中的机理公式对毛刺状态进行初步评估，基于毛刺高度分类表对机理层面的预测得出分类结果，并与实际测得的毛刺高度分类结果进行比对，实际分类等级与机理预测分类等级的差值作为标签．动态更新模型以钻削阶段的力信号、加工参数作为输入，判断基于机理预测的毛刺分类等级是否须要调高、不变或者降低，以此对神经网络模型进行训练，从而最终的预测结果为数据模型对机理模型的更新迭代，即机理模型预测结果y1与动态更新模型结果y2之和．制孔过程的状态参数为时序数据，因此须要采用时序信息提取能力较好的模型．本研究动态更新模型神经网络部分采用注意力机制+双向长短时记忆神经网络+支持向量机(attention+Bi-LSTM+SVM)对输入数据进行训练．其中注意力机制(attention)可以对模型中不同的关注部分赋予不同的权重，从而优化模型以做出更好的判断；长短时记忆(long short-term memory，LSTM)神经网络是循环神经网络的一个变种，可以解决传统循环神经网络(recurrent neural network，RNN)对数据输入的长期依赖问题，而Bi-LSTM是双向循环结构，相较于LSTM增加了从未来到过去的数据流向，可以更好地发掘数据的时序特征，其结构如图6～7所示；支持向量机(support vector machine，SVM)适合处理小样本、非线性等问题，同时由于其采用结构最小化，因此可以得到全局最10.13245/j.hust.220616.F006图6双向长短时记忆结构网络10.13245/j.hust.220616.F007图7长短时记忆单元结构优解．由于经过降噪之后的原始数据信息容量巨大，不能直接用于神经网络的训练，须要将信号进行特征提取、归一化等操作之后才能反映毛刺状态，因此该动态更新模型的输入为当前制孔条件下的主轴转速、进给量及三向钻削力与轴向转矩的时域特征值，输出为实际毛刺高度分类等级与机理模型预测分类等级之间的差值，即是否须要提高或降低毛刺预测分类等级．动态更新模型结构如图8所示，表2是动态更新模型的预测输出分类情况．10.13245/j.hust.220616.F008图8动态更新模型结构10.13245/j.hust.220616.T002表2模型预测输出标签序号预测结果误差标签1降低预测等级-12预测等级不变03提高预测等级1神经网络的输入信号结构为120×18，120为经过时域分析的120个时间截断，18则为4×4+2，即16种时域特征参数加上主轴转速与进给量．输出的结果为三分类结果，即提升毛刺高度分类等级、保持分类等级不变及降低毛刺分类预测等级．3 机器人制孔实验3.1　关节刚度辨识实验机器人关节刚度主要由伺服电机、减速机、齿轮、传动轴等组成，可以通过传动部件的结构、传动比等对关节刚度进行计算．然而在上述关节组成部分具体参数未知的情况下，须要采用间接方式对关节刚度进行辨识．基于文献[15]的方法，通过在末端施加不同的负载，由末端上安装的激光跟踪仪监测施加负载前后的末端位置变化，基于雅可比矩阵及最小二乘法处理数据，从而得到机器人的关节刚度．具体实验步骤如下：a. 通过示教器将机器人调整至合适的位姿，将调整好的姿态作为初始姿态，记录初始状态下的靶球位置、力传感器数值；b. 向机器人末端施加不同质量的负载(5 kg，10 kg)，记录施加负载之后的靶球位置及力传感器数值；c. 完成步骤a中姿态测量之后通过示教器变换另一种姿态，重复前两个步骤直至测量完成；d. 整理实验仪器和试验结果，并对结果进行分析．在经过试验之后，可以得到不同位姿下施加载荷前后的三个方向上的位移∆D及三向力∆F，采用最小二乘法拟合得到机器人的关节刚度Kq，      ΔD=[dx,dy,dz]T=[J(q)Kq-1JT(q)F]3×n=[J(q)diag(JT(q)F)Cq']3×n,式中：J(q)为机器人雅可比矩阵；Cq'=(1/k1，1/k2，1/k3，1/k4，1/k5，1/k6)，ki为各关节刚度；3×n表示取各矩阵的前3行．最终结合建立的机器人D-H模型得到机器人的关节刚度矩阵为Kq=diag(1.71×1010，1.7×1010，2.95×109，7.23×108，3.74×106，1.91×107)．3.2　制孔实验本研究将传感器安装在制孔系统合适位置采集信号，同时记录不同钻削条件下的转速和进给量．机器人制孔系统的组成主要包括：机器人本体，型号为安川MOTOMAN-GP110，负载110 kg；末端执行器，主要包含电主轴、进给伺服电机等；六维力传感器，型号为GGT92004BAO，量程为Fx/Fy为1 000 N，Fz为2 000 N，Tx/Ty/Tz为100 N·m；加速度传感器，采集器型号为东方所INV3062T．通过该平台进行直径为6 mm的钻孔实验，加工材料为200 mm×100 mm×10 mm的奥氏体不锈钢，牌号为1Cr18Ni9Ti．通过六维力传感器采集力和转矩，通过加速度传感器采集振动信号，利用以上数据进行毛刺高度预测模型的搭建．钻孔加工时机器人位姿六个关节角分别为-0.45°,7.30°,-45.29°,0.57°,52.99°,-2.05°，部分实验参数见表3，以主轴转速与进给量为因子进行27组全因素水平实验，每组实验进行3次．10.13245/j.hust.220616.T003表3实验参数实验序号主轴转速/(r∙min-1)进给量/(mm∙s-1)11 8000.2021 8000.2531 8000.30⋮⋮⋮252 4000.20262 4000.25272 4000.30毛刺高度的测量通过基恩士VHX-950F数码显微系统进行，基于测量仪器测得结果如图9所示．毛刺高度测量方法为：利用显微系统对毛刺进行3D形貌扫描建模，以工件上表面作为基准面，以多次测得毛刺最高点处均值作为毛刺高度的实验数据．10.13245/j.hust.220616.F009图9毛刺高度测量实例我国冲压件毛刺高度采用GB/T33217—2016标准，国际上也有许多国际组织对毛刺边缘的平整度提出了相关规定，由于毛刺高度评价分类评级种类较多，本研究根据上述标准针对弱刚性制孔系统毛刺高度分类评级，以毛刺高度占工件总高度的百分比作为界限对毛刺进行分类，具体分类情况见表4．10.13245/j.hust.220616.T004表4毛刺高度分类毛刺分类毛刺高度占比评级介绍第一级5%可视为无须去毛刺操作第二级5%~7%须去毛刺的合格孔第三级7%须去毛刺并评价孔是否合格由于初始采集到的数据量十分庞大，且并非整个时间线上采集的数据对训练都是有用的，因此须要对采集到的数据进行特征提取，并进行归一化预处理．在机器人制孔系统中，选取与制孔质量相关性较大的几种信号采集输入，包括x，y，z三个方向上的力信号和轴向转矩信号，力和转矩与钻削质量直接关联，过大的钻削力会引起工件、刀具及整个制孔系统发生应变，产生不利影响．因此，更新模型的输入最终选择为当前加工状态下的主轴转速、进给量及三向力和转矩的时域特征值．其中力与转矩信号提取其均值、均方根值、方差及峰值这4种时域特征，总共4×4=16种时域特征参数．将经过处理的数据与其对应的标签输入神经网络中进行训练，从而实现机理与数据融合的毛刺高度分类预测．4 实验验证通过搭建的机器人制孔系统控制软件，在机器人平台进行实验与仿真验证．首先基于建立的制孔毛刺高度预测静态物理模型对毛刺高度进行计算，并根据表4进行毛刺高度评级分类．其次将处理过的实验采集数据输入所建立的动态更新神经网络模型，输出结果为实际毛刺高度评级与物理模型预测毛刺高度评级之差．基于物理模型计算结果与其对应误差标签举例见表5．10.13245/j.hust.220616.T005表5基于静态物理模型毛刺高度分类举例实验序号基于物理模型毛刺高度结果/mm基于物理模型毛刺高度分类实际测得毛刺高度分类预测误差标签10.133第一级第一级不变20.414第一级第一级不变30.441第一级第一级不变40.506第二级第二级不变50.566第二级第三级提高60630第二级第一级降低70.689第二级第二级不变80.658第二级第一级降低90.390第一级第一级不变将实验获得的加工信号时序数据集、实验中测量的刀具磨损状态数据组成带有标签的数据集，并对其进行分类，剩余的数据作为测试集．数据集的标签为每组数据实际测量出的毛刺高度分类等级与基于静态物理模型得出的毛刺高度分类等级的误差，最终得到模型训练损失-精度情况如图10所示．以单纯采用物理模型与数据模型作为对照组，经过对比可以发现机理与数据融合预测模型准确率最高，为94.70%，相比于纯物理模型与单纯数据模型更高，单纯采用物理模型及单纯采用数据模型预测准确率分别为67.00%和78.95%．而由于钻削过程中的复杂程度，纯物理模型的预测分类效果较差，从而验证了该预测方法的可行性．10.13245/j.hust.220616.F010图10模型训练结果5 结语本研究针对机器人制孔系统钻削毛刺高度状态预测问题，建立了机理数据融合的毛刺高度预测模型．a. 针对机器人制孔系统相较于传统机床刚性较弱这一特点，探究了机器人制孔系统颤振机理，并基于钻削力和振动理论模型建立了机器人制孔出口毛刺高度预测机理模型．b. 以毛刺高度预测机理模型为基础，建立了制孔毛刺高度分类动态更新模型，机理与数据融合预测的方法在保留机理模型可解释性的同时提高了模型的泛化性和准确性．本研究通过自行搭建的机器人制孔控制原型系统，以机器人制孔系统钻削过程的毛刺生成机理模型建立行为衍射模型，并建立以加工采集数据为输入的动态更新模型，从而建立钻削毛刺生成的虚拟仿真系统，对机器人制孔系统进行监控．通过原型系统对所建立的预测方法进行验证，经过实验数据的对比发现该模型预测准确率达到94.74%，表明该预测方法具有一定的实际应用性能，同时可以对该种弱刚性类型系统进行监测，为提升加工质量奠定基础．