大型龙门铺丝机综合运动学建模及参数标定

Integrated Kinematics Modeling and Parametric Calibration of Large Gantry Fiber Placement Machines

摘要

自动化设备通常是一个强耦合、非线性且具备一定柔性的复杂机械系统，这使得建立该系统精确的运动学模型具有较大的难度。为此，针对一台大型龙门自动铺丝机提出一种综合几何偏差和重力变形的运动学建模方法。基于多体系统理论建立理论运动学模型，分析模型坐标转换中存在的几何偏差，研究龙门铺丝机重力变形敏感部位，通过力学抽象和参数提取得到重力变形矩阵，将几何偏差与重力变形有机引入理论运动学模型，建立龙门铺丝机综合运动学模型。利用目标函数建立、测量点选择、参数辨识等手段实现模型中的参数标定。实验结果表明，所提方法可有效提高龙门铺丝机运动学模型的精度，采用该方法后龙门铺丝机位置误差降低了74%以上。

Abstract

Automation equipment was typically a strongly coupled， nonlinear， and flexible complex mechanical system， which made it challenging to establish an accurate kinematics model. A kinematics modeling method was proposed that integrated geometric deviations and gravity deformations for a large gantry automatic fiber placement machine. A theoretical kinematics model was developed based on multi-body system theory， the geometric deviations were analyzed in model coordinate conversion， and the gravity deformation-sensitive parts of the gantry fiber placement machines were studied. The gravity deformation matrix was obtained through mechanics abstraction and parameter extraction. The geometric deviations and gravity deformations were incorporated into the theoretical kinematics model to establish a comprehensive kinematics model of the gantry fiber placement machines. Parametric calibration was achieved through objective function establishment， measurement point selection， and parameter identification. Experimental results demonstrate that the proposed method may significantly improve the precision of the kinematics model， the position errors of the gantry fiber placement machines are reduced by over 74%.

Graphical abstract

关键词

龙门铺丝机 / 运动学建模 / 几何偏差 / 重力变形 / 参数辨识

Key words

gantry fiber placement machine / kinematics modeling / geometric deviation / gravity deformation / parameter identification

0 引言

碳纤维增强树脂基（carbon fiber reinforced plastics， CFRP）复合材料具有高比强度、高比模量、抗腐蚀性、耐疲劳性、可设计性、成形工艺多样性等优点，广泛应用于航空航天、轨道交通、风力发电等领域^［1-3］。自动铺丝技术是复合材料构件低成本、数字化、高效率、高质量先进制造技术，已成为大型复杂复合材料部件的典型制造工艺技术^［4］。

然而，在实际情况下，受到设备零件制造误差、部件安装偏差、重力变形、热变形等多因素耦合作用的影响，自动铺放设备末端实际位置和姿态与关节运动之间存在强烈的非线性关系，导致设备末端精度极难控制^［5］。其中，零件制造误差、部件安装偏差都属于几何误差，关于几何误差的研究方法通常有建模^［6-8］、测量、辨识^［9-10］、补偿^［10-12］等方面。郭世杰等^［9］提出了一种基于平面光栅的面—线机床空间几何误差测量与辨识方法，通过单轴和两轴联动在3个测量平面内分别连续测量一次直线轨迹，可快速实现三轴机床21项几何误差的准确辨识；李忠群等^［13］以一台四轴立式加工中心为研究对象，利用拉丁超立方采样和 Sobol 全局灵敏度分析方法分析了机床几何误差对机床空间精度的影响，并辨识出其中的关键误差项，以提高机床精度补偿的效率。

但是上述文献都是基于刚体假设来研究设备的运动学精度问题，忽略了力误差，由力导致的力误差经过运动链传递与放大的作用将会成为影响末端精度的重要因素。王明海等^［14］对弱刚性工业机器人开展了运动学分析、工作空间分析及刚度分析，采用参数优化后的粒子群算法来寻找机器人刚度最大的加工位姿，从而减小机械臂变形对末端精度的影响；陆历历等^［15］利用有限元软件得出一台大型数控落地镗铣加工中心滑枕挠曲变形的分布规律和滑枕行程的关系，并将其应用在滑枕的变形补偿中，有效降低了滑枕行程中的下垂量。

由于复合材料铺放需要恒温恒湿的工作环境，故自动铺丝设备的热误差可以暂时忽略不计。但几何误差和重力变形对其末端精度的影响却不可忽略，两者存在严重的强耦合关系，且通过测量仪器采集到的数据难以将两者有效地分离。若不对两者从机理上开展分析，则将会直接影响自动铺丝设备的铺放精度。

为提高一台大型龙门铺丝机的精度，本文研究了一种综合几何偏差和重力变形的运动学建模及参数标定方法。简单介绍了龙门铺丝机的结构组成；然后基于多体系统理论建立了龙门铺丝机理论运动学模型；详细分析了关于几何偏差和重力变形的运动学参数引入过程，并建立了龙门铺丝机综合运动学模型；介绍了运动学参数标定中目标函数建立、测量点选择、参数辨识等过程；最后通过实验对龙门铺丝机综合运动学模型和运动学参数标定方法进行了验证。

1 龙门铺丝机综合运动学建模

1.1 龙门铺丝机结构

本文研究的龙门自动铺丝机如图1所示，该龙门铺丝机具有6个运动自由度，主要由横梁、左右立柱、Y轴滑台、Z轴滑枕、液压平衡缸、蓄能器组、旋转轴框架、B轴、A轴、C轴和16丝束铺丝头组成控制。

从运动轴的连接方式来看，该龙门铺丝机属于典型的串联式机床，运动轴间的拓扑连接关系可以表示为X-Y-Z-B-A-C的形式。其中，X轴、Y轴、Z轴为相互正交的平动轴；B轴、A轴、C轴为旋转轴，分别绕Y、X、Z方向旋转。此外，该铺丝机的龙门架高度达到8 m，横梁长度接近12.7 m，铺丝头质量约有1000 kg，因此该龙门铺丝机具有大行程、大跨度、重载荷的结构特点。

1.2 理论运动学建模

将多种误差因素以运动学参数的形式引入运动学模型中可有效提高模型精度。在该过程中，如何从这些误差因素中抽象、分离出合适的运动学参数是关键，本小节就该问题展开具体研究。

龙门铺丝机设有三个相互垂直的平动轴，针对其结构构型和运动形式，采用多体系统理论完成各坐标系的建立，如图2a所示。所建的坐标系包括：基坐标系\{S_{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}})，运动轴坐标系\{S_X\}、\{S_Y\}、\{S_Z\}、\{S_B\}、\{S_A\}、\{S_C\}，铺丝头末端坐标系\{S_{\mathrm{H}}\}和模具坐标系\{S_{\mathrm{M}}\}。其中，基坐标系是三维空间内的固定坐标系；运动轴坐标系固连在各运动轴上，用于描述各轴的运动；铺丝头末端坐标系放置在铺丝头末端处；模具坐标系用于表征模具在三维空间内的位置和姿态。

根据上述各坐标系的定义（图2b），利用齐次变换矩阵描述相邻坐标系的变换，可建立龙门铺丝机理论运动学模型，其表达式如下：

{}_{\mathrm{H}}{}^{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}}\mathit{T}(\mathit{q})={\mathit{T}}_X(x){\mathit{T}}_Y(y){\mathit{T}}_Z(z)\cdot

{\mathit{T}}_B(\beta ){\mathit{T}}_A(\alpha ){\mathit{T}}_C(\gamma ){\mathit{T}}_{\mathrm{H}}（1）

其中，{}_{\mathrm{H}}{}^{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}}\mathit{T}(\mathit{q})表示坐标系\left\{S_{\mathrm{H}}\right\}在坐标系\left\{S_{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}}\right\}下的位置和姿态，\mathit{q}=(x,y,z,\alpha ,\beta ,\gamma )为各轴运动量向量，{\mathit{T}}_X、{\mathit{T}}_Y、{\mathit{T}}_Z、{\mathit{T}}_B、{\mathit{T}}_A、{\mathit{T}}_C、{\mathit{T}}_{\mathrm{H}}均为齐次变换矩阵，分别描述各轴运动量和铺丝头相对位置，具体可表示为

{\mathit{T}}_X(x)=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& x\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]{\mathit{T}}_Y(y)=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& y\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]

{\mathit{T}}_Z(z)=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& z\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]

{\mathit{T}}_B(\beta )=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}\beta & \mathrm{0}& \mathrm{s}\mathrm{i}\mathrm{n}\beta & \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}\beta & \mathrm{0}& \mathrm{c}\mathrm{o}\mathrm{s}\beta & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]

{\mathit{T}}_A(\alpha )=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{c}\mathrm{o}\mathrm{s}\alpha & -\mathrm{s}\mathrm{i}\mathrm{n}\alpha & \mathrm{0}\\ \mathrm{0}& \mathrm{s}\mathrm{i}\mathrm{n}\alpha & \mathrm{c}\mathrm{o}\mathrm{s}\alpha & -L_{\mathrm{a}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]

{\mathit{T}}_C(\gamma )=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}\gamma & -\mathrm{s}\mathrm{i}\mathrm{n}\gamma & \mathrm{0}& \mathrm{0}\\ \mathrm{s}\mathrm{i}\mathrm{n}\gamma & \mathrm{c}\mathrm{o}\mathrm{s}\gamma & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]

{\mathit{T}}_{\mathrm{H}}=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& -L_{\mathrm{h}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]

式中：L_{\mathrm{a}}为坐标系\left\{S_A\right\}与\left\{S_B\right\}之间的偏置；L_{\mathrm{h}}为坐标系\left\{S_{\mathrm{H}}\right\}与\left\{S_C\right\}之间的偏置；L_{\mathrm{a}}、L_{\mathrm{h}}两个参数的理论值可从龙门铺丝机数模中获取。

1.3 几何偏差分析

几何偏差是指零部件在制造、安装过程中引入的结构尺寸偏差、位置误差、转角误差等，下文对各坐标转换中存在的几何偏差进行分析。

1）基坐标系。在实际操作中，基坐标系通常采用激光跟踪仪等测量仪器建立。受测量误差等因素影响，实际建立的基坐标系与理论基坐标系会存在微小偏差。引入实际基坐标系\left\{S_{\mathrm{A}\mathrm{c}\mathrm{t}}\right\}，该坐标系与理论基坐标系的偏差可用6个参数描述，分别为沿X、Y、Z方向的位置误差和绕X、Y、Z方向的转角误差，可表示为(x_{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}},y_{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}},z_{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}},{\alpha }_{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}},{\beta }_{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}},{\gamma }_{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}})。

2）各运动轴局部坐标系。已知每个运动轴在三维空间中有6个自由度，经过一个空间运动后运动轴会在6个自由度方向上产生偏差（即3个位置偏差和3个角度偏差）。但根据参数冗余性要求，每个运动轴对应的6个参数中存在冗余项，因此需要对各运动轴6个参数进行分析，以确定最终的运动学参数项。

3）平动轴坐标系。平动轴坐标系所处空间位置不影响平移运动的描述，且绕该轴的转角误差不影响该轴的最终位置，因此由平动轴引入的运动学参数仅为与该轴不平行的两项转角误差。考虑到X轴坐标系与基坐标系重合，可以认为X轴运动不引入运动学参数，而由Y轴运动引入的运动学参数为绕X、Z方向的转角误差({\theta }_Y^{(X)},{\theta }_Y^{(Z)})，由Z轴运动引入的运动学参数为绕X、Y方向的转角误差({\theta }_Z^{(X)},{\theta }_Z^{(Y)})。

4）旋转轴坐标系。旋转轴坐标系建立在该轴旋转轴线上，沿该轴轴线方向的位置变化不影响对其旋转运动的描述，且该轴旋转方向不引入转角误差。考虑到B轴坐标系原点与上一级运动轴Z轴坐标系原点重合，可以认为B轴运动只引入绕X、Z方向的转角误差({\theta }_B^{(X)},{\theta }_B^{(Z)})；A轴和B轴坐标系原点均处于C轴旋转轴线上，可以认为A轴运动引入了Z方向的位置误差和绕Y、Z方向的转角误差(z_A,{\theta }_A^{(Y)},{\theta }_A^{(Z)})；C轴运动则引入了沿X、Y方向的位置误差和绕X、Y方向的转角误差(x_C,y_C,{\theta }_C^{(X)},{\theta }_C^{(Y)})。

根据表1中各坐标系对应的齐次变换矩阵，可以建立结合几何偏差的龙门铺丝机运动学模型，其表达式如下：

{}_{\mathrm{H}}{}^{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}}\mathit{T}(\mathit{q},{\mathit{p}}_{\mathrm{d}})={}_{\mathrm{A}\mathrm{c}\mathrm{t}}{}^{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}}\mathit{T}{}_X{}^{\mathrm{A}\mathrm{c}\mathrm{t}}\mathit{T}{}_Y{}^X\mathit{T}{}_Z{}^Y\mathit{T}{}_B{}^Z\mathit{T}{}_A{}^B\mathit{T}{}_C{}^A\mathit{T}{}_{\mathrm{H}}{}^C\mathit{T}（2）

其中，{\mathit{p}}_{\mathrm{d}}为由几何偏差引入的22个运动学参数组成的向量，{}_{\mathrm{A}\mathrm{c}\mathrm{t}}{}^{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}}\mathit{T}为理论基坐标系\left\{S_{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}}\right\}到实际基坐标系\left\{S_{\mathrm{A}\mathrm{c}\mathrm{t}}\right\}的变换矩阵，{}_X{}^{\mathrm{A}\mathrm{c}\mathrm{t}}\mathit{T}为\left\{S_{\mathrm{A}\mathrm{c}\mathrm{t}}\right\}到\left\{S_X\right\}的变换矩阵，{}_Y{}^X\mathit{T}为\left\{S_X\right\}到\left\{S_Y\right\}的变换矩阵，其他同理。

1.4 结构重力变形研究

相较于普通机床，龙门铺丝机的工作行程与结构尺寸大幅增大，龙门架高8 m、横梁跨度近12.7 m、Z轴滑枕行程超2 m；龙门铺丝机的末端执行器铺丝头因集成了大量预浸料丝束控制机构，使得铺丝头质量近1000 kg，含Z轴滑枕、旋转轴框架、各运动轴及铺丝头的垂直运动组件，总质量更达10 t。

重力载荷下，龙门铺丝机易发生重力变形的结构为旋转轴框架、Y轴滑台及龙门架横梁，分布如图3所示，原因如下：①旋转轴框架为铝质，弹性模量小，受B轴、A轴和铺丝头偏置重力易变形；②Y轴滑台为龙门架与垂直运动组件的连接部位，双侧导轨滑块的连接形式刚度较差；③龙门架横梁跨度长、自重大、载荷复杂，且承受垂直运动组件重量，极易弯曲变形。

1.4.1 旋转轴框架分析

旋转轴的旋转运动是龙门铺丝机重心位置变化的主要来源，也是旋转轴框架结构变形变化的主要因素。考虑到旋转轴框架结构变形的复杂性，本文在分析其重力变形时对其进行了力学简化处理，在旋转轴框架与Z轴滑枕的连接处P引入3个方向的线性弹簧和绕3个方向的扭簧来反映该结构变形对末端精度的影响，其中线性弹簧弹性系数用(k_X^{(P)},k_Y^{(P)},k_Z^{(P)})表示，扭簧弹性系数用（r_X^{(P)},r_Y^{(P)},r_Z^{(P)}）表示，如图4所示。

图4中l_{PB}为P点到B轴坐标系原点的距离，l_{BA}为B轴坐标系原点到A轴坐标系原点的距离。对旋转轴框架进行力学简化后，基于力学原理可计算出B轴、A轴在不同位置时P点处沿3个方向的力：

\left.\begin{array}{l}{F}_P^{\left(X\right)}=\mathrm{0}\\ F_P^{\left(Y\right)}=\mathrm{0}\\ F_P^{\left(Z\right)}=g(m_A+m_B+m_{\mathrm{r}\mathrm{f}})\end{array}\right\}（3）

以及绕3个方向的力矩：

\left.\begin{array}{l}{N}_P^{\left(X\right)}=g(m_B{y}_B+m_A{y}_A\mathrm{c}\mathrm{o}\mathrm{s}\alpha -m_A{z}_A\mathrm{s}\mathrm{i}\mathrm{n}\alpha )\\ \begin{array}{l}{N}_P^{\left(Y\right)}=-g[(m_A{x}_B+m_B{x}_B)\mathrm{c}\mathrm{o}\mathrm{s}\beta +(l_{BA}{m}_A+m_B{z}_B)·\\ \mathrm{s}\mathrm{i}\mathrm{n}\beta +m_A{z}_A\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta +m_A{y}_A\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta ]\end{array}\\ N_P^{\left(Z\right)}=\mathrm{0}\end{array}\right\}（4）

式中：g为重力加速度；m_A、m_B、m_{\mathrm{r}\mathrm{f}}分别为A轴、B轴和旋转轴框架的质量；\alpha、\beta分别为B轴、A轴的运动量；(x_A,y_A,z_A)为A轴重心相对于A轴局部坐标系的位置；(x_B,y_B,z_B)为B轴重心相对于B轴局部坐标系的位置。

基于上述计算得到的力和力矩，可计算出旋转轴框架的静力变形：

\left.\begin{array}{l}\mathrm{\Delta }{z}_P=F_P^{\left(Z\right)}/k_Z^{\left(P\right)}\\ \mathrm{\Delta }{\theta }_P^{\left(X\right)}=N_P^{\left(X\right)}/r_X^{\left(P\right)}\\ \mathrm{\Delta }{\theta }_P^{\left(Y\right)}=N_P^{\left(Y\right)}/r_Y^{\left(P\right)}\end{array}\right\}（5）

式中：\mathrm{\Delta }{z}_P为P点在Z方向的重力变形；\mathrm{\Delta }{\theta }_P^{\left(X\right)}为P点处重力引起的绕X方向角度偏差；\mathrm{\Delta }{\theta }_P^{\left(Y\right)}为P点处重力引起的绕Y方向角度偏差。

1.4.2 Y轴滑台分析

对Y轴滑台连接部位Q进行受力分析时，同样将其进行力学简化处理，引入3个方向的线性弹簧（线性弹簧弹性系数为(k_X^{(Q)},k_Y^{(Q)},k_Z^{(Q)})）和绕3个方向的扭簧（扭簧弹性系数为(r_X^{(Q)},r_Y^{(Q)},r_Z^{(Q)})）来模拟Q点处的结构变形，如图5所示，其中l_{PQ}^{\left(X\right)}为Q点与P点在X方向上的距离，l_{PQ}^{\left(Z\right)}为龙门铺丝机Z轴处于零位时Q点与P点在Z方向上的距离。

连接部位Q点处不仅受到Y轴滑台、Z轴滑枕的重力，还在连接点P处受到旋转轴框架的作用反力。根据静力学原理，可计算出Q点处的受力情况：

\left.\begin{array}{l}{F}_Q^{\left(Z\right)}=F_P^{\left(Z\right)}+g\left(m_Z+m_Y\right)\\ N_Q^{\left(X\right)}=N_P^{\left(X\right)}\\ N_Q^{\left(Y\right)}=N_P^{\left(Y\right)}-m_{Z}g{l}_{PQ}^{\left(X\right)}\end{array}\right\}（6）

其中，m_Z、m_Y分别为Z轴滑枕和Y轴滑台的质量。则可计算得到Y轴滑台处的静力变形：

\left.\begin{array}{l}\mathrm{\Delta }{z}_Q=F_Q^{\left(Z\right)}/k_Z^{\left(Q\right)}\\ \mathrm{\Delta }{\theta }_Q^{\left(X\right)}=N_Q^{\left(X\right)}/r_X^{\left(Q\right)}\\ \mathrm{\Delta }{\theta }_Q^{\left(Y\right)}=N_Q^{\left(Y\right)}/r_Y^{\left(Q\right)}\end{array}\right\}（7）

式中：\mathrm{\Delta }{z}_Q为Q点在Z方向的重力变形；\mathrm{\Delta }{\theta }_Q^{\left(X\right)}为Q点处重力引起的绕X方向角度偏差；\mathrm{\Delta }{\theta }_Q^{\left(Y\right)}为Q点处重力引起的绕Y方向角度偏差。

1.4.3 龙门架横梁分析

为分析横梁自重和其他部件质量对横梁结构的影响，本文将横梁简化为简支梁结构，其弯曲变形可视为均布载荷q_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}和集中载荷F_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}作用叠加的结果，如图6所示，其中均布载荷模拟横梁自重，集中载荷模拟其他部件的重力载荷。

当施加载荷为均布载荷和集中载荷时，梁的挠曲线方程分别为

\left.\begin{array}{l}{w}_{\mathrm{u}}(x)=-\frac{q_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}x}{\mathrm{24}EI}(l^{\mathrm{3}}-\mathrm{2}l{x}^{\mathrm{2}}+x^{\mathrm{3}})\mathrm{0}\le x\le l\\ w_{\mathrm{c}}(x)=-\frac{F_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}bx}{\mathrm{6}EIl}(l^{\mathrm{2}}-x^{\mathrm{2}}-b^{\mathrm{2}})\mathrm{0}\le x\le a\end{array}\right\}（8）

式中：w_{\mathrm{u}}为均布载荷下梁的挠曲线方程；w_{\mathrm{c}}为集中载荷下梁的挠曲线方程；l为梁的长度；E为梁的弹性模量；I为梁的惯性矩；a、b分别为集中载荷到两侧支撑点的距离，如图6所示。

求挠曲线转角\theta时，由于梁的挠度一般远小于跨度，挠曲线是一条平坦的曲线，转角\theta也是一个很小的角度，因此该角度可表示为

\theta =\mathrm{t}\mathrm{a}\mathrm{n}\theta =\mathrm{d}w/\mathrm{d}x=w^{\text{'}}(x)（9）

其中，w表示式（8）中的两个挠曲线方程，则两种载荷下的梁挠曲线转角方程为

\left.\begin{array}{l}{\theta }_{\mathrm{u}}(x)=-\frac{q_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}}{\mathrm{24}EI}(l^{\mathrm{3}}-\mathrm{6}l{x}^{\mathrm{2}}+\mathrm{4}{x}^{\mathrm{3}})\mathrm{0}\le x\le l\\ {\theta }_{\mathrm{c}}(x)=-\frac{F_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}b}{\mathrm{6}EIl}(l^{\mathrm{2}}-\mathrm{3}{x}^{\mathrm{2}}-b^{\mathrm{2}})\mathrm{0}\le x\le a\end{array}\right\}（10）

根据式（8）和式（10），基于简支梁小变形叠加原理可以计算出Z轴滑枕处（即x=a时）横梁的挠度变形和转角变形：

\left.\begin{array}{l}{w}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}(a)=w_{\mathrm{u}}(a)+w_{\mathrm{c}}(a)\\ {\theta }_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}(a)={\theta }_{\mathrm{u}}(a)+{\theta }_{\mathrm{c}}(a)\end{array}\right\}（11）

此时，再令

\left.\begin{array}{l}a=\frac{l}{\mathrm{2}}-y\\ b=\frac{l}{\mathrm{2}}+y\end{array}\right\}（12）

可将式（11）中挠度变形和转角变形转换成关于Y轴滑台位置y的函数，则有

\left.\begin{array}{l}{w}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}(y)=-[(\mathrm{128}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}+\mathrm{16}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}l)y^{\mathrm{4}}-(\mathrm{64}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{2}}+\\ \mathrm{24}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{3}})y^{\mathrm{2}}+(\mathrm{5}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{5}}+\mathrm{8}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{4}})]/(\mathrm{384}EIl)\\ {\theta }_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}(y)=[(\mathrm{4}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}l+\mathrm{16}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}})y^{\mathrm{3}}-(\mathrm{3}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{3}}+\\ \mathrm{4}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{2}})y]/(\mathrm{24}EIl)\end{array}\right\}（12）

由此得到横梁在均布载荷和集中载荷下发生的弯曲变形。

1.4.4 综合运动学建模

为方便描述各处结构变形对末端精度的影响，根据小变形假设将各变形统一转换到基坐标系下，其表达式如下：

\left.\begin{array}{l}\mathrm{\Delta }\alpha (y,z,\beta ,\alpha )=\mathrm{\Delta }{\theta }_Q^{(X)}+\mathrm{\Delta }{\theta }_P^{(X)}+{\theta }_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}\\ \mathrm{\Delta }\beta (y,z,\beta ,\alpha )=\mathrm{\Delta }{\theta }_Q^{(Y)}+\mathrm{\Delta }{\theta }_P^{(Y)}\\ \mathrm{\Delta }\gamma (y,z,\beta ,\alpha )=\mathrm{0}\\ \mathrm{\Delta }x(y,z,\beta ,\alpha )=\mathrm{\Delta }\beta (l_{PQ}^{(Z)}-z+l_{PB})-\mathrm{\Delta }{\theta }_P^{(Y)}(l_{PQ}^{(Z)}-z)\\ \mathrm{\Delta }y(y,z,\beta ,\alpha )=\mathrm{\Delta }\alpha (l_{PQ}^{(Z)}-z+l_{PB})-\mathrm{\Delta }{\theta }_P^{(X)}(l_{PQ}^{(Z)}-z)\\ \mathrm{\Delta }z(y,z,\beta ,\alpha )=-\mathrm{\Delta }{z}_Q-\mathrm{\Delta }{z}_P+w_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}\end{array}\right\}（14）

式（14）虽然统一了各处变形的表征，但基坐标系是在各轴处于零位时建立的，而建立该坐标系时龙门铺丝机的重力变形已经存在，因此重力变形随各轴位置的变化量才是真正影响末端精度的最终因素。重力变形随各轴位置的变化量可表示为

\left.\begin{array}{l}\mathrm{\Delta }{\alpha }^{\text{'}}=\mathrm{\Delta }\alpha (y,z,\beta ,\alpha )-\mathrm{\Delta }\alpha (\mathrm{0,0},\mathrm{0,0})\\ \mathrm{\Delta }{\beta }^{\text{'}}=\mathrm{\Delta }\beta (y,z,\beta ,\alpha )-\mathrm{\Delta }\beta (\mathrm{0,0},\mathrm{0,0})\\ \mathrm{\Delta }{\gamma }^{\text{'}}=\mathrm{\Delta }\gamma (y,z,\beta ,\alpha )-\mathrm{\Delta }\gamma (\mathrm{0,0},\mathrm{0,0})\\ \mathrm{\Delta }{x}^{\text{'}}=\mathrm{\Delta }x(y,z,\beta ,\alpha )-\mathrm{\Delta }x(\mathrm{0,0},\mathrm{0,0})\\ \mathrm{\Delta }{y}^{\text{'}}=\mathrm{\Delta }y(y,z,\beta ,\alpha )-\mathrm{\Delta }y(\mathrm{0,0},\mathrm{0,0})\\ \mathrm{\Delta }{z}^{\text{'}}=\mathrm{\Delta }z(y,z,\beta ,\alpha )-\mathrm{\Delta }z(\mathrm{0,0},\mathrm{0,0})\end{array}\right\}（15）

于是，采用齐次变换矩阵表征各变形变化量，可得到龙门铺丝机重力变形矩阵：

\mathrm{\Delta }{\mathit{T}}_{\mathrm{g}}={\mathit{T}}_X(\mathrm{\Delta }{x}^{\text{'}}){\mathit{T}}_Y(\mathrm{\Delta }{y}^{\text{'}}){\mathit{T}}_Z(\mathrm{\Delta }{z}^{\text{'}})\cdot

{\mathit{T}}_A(\mathrm{\Delta }{\alpha }^{\text{'}}){\mathit{T}}_B(\mathrm{\Delta }{\beta }^{\text{'}}){\mathit{T}}_C(\mathrm{\Delta }{\gamma }^{\text{'}})（16）

将重力变形矩阵\mathrm{\Delta }{\mathit{T}}_{\mathrm{g}}结合至龙门铺丝机运动学模型时，还需分析该矩阵引入的运动学参数。以重力变形分量\mathrm{\Delta }{\alpha }^{\text{'}}为例，将其表达式代入后可得

\begin{array}{c}\mathrm{\Delta }{\alpha }^{\text{'}}=[m_{A}g{y}_A(\mathrm{c}\mathrm{o}\mathrm{s}\alpha -\mathrm{1})-m_{A}g{z}_A\mathrm{s}\mathrm{i}\mathrm{n}\alpha ]/r_X^{(P)}+\\ (m_{A}g{y}_A(\mathrm{c}\mathrm{o}\mathrm{s}\alpha -\mathrm{1})-m_{A}g{z}_A\mathrm{s}\mathrm{i}\mathrm{n}\alpha )/r_X^{(Q)}+\end{array}

[(\mathrm{4}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}l+\mathrm{16}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}})y^{\mathrm{3}}-(\mathrm{3}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{3}}+\mathrm{4}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{2}})y]/(\mathrm{24}EIl)（17）

为了减少式（17）中未知参数的引入且避免参数冗余，令r_X^{(Q)}=k_{\alpha }{r}_X^{(P)}，k_{\alpha }为两个线性弹簧弹性系数r_X^{(P)}和r_X^{(Q)}的比值，则式（17）可变为

\begin{array}{c}\mathrm{\Delta }{\alpha }^{\text{'}}=(\mathrm{1}/k_{\alpha }+\mathrm{1})[m_{A}g{y}_A(\mathrm{c}\mathrm{o}\mathrm{s}\alpha -\mathrm{1})-\\ m_{A}g{z}_A\mathrm{s}\mathrm{i}\mathrm{n}\alpha ]/r_X^{(P)}+[(\mathrm{4}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}l+\mathrm{16}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}})y^{\mathrm{3}}-\end{array}

(\mathrm{3}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{3}}+\mathrm{4}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{2}})y]/(\mathrm{24}EIl)（18）

此时，可引入4个运动学参数(a_{\mathrm{1}},a_{\mathrm{2}},a_{\mathrm{3}},a_{\mathrm{4}})：

\left.\begin{array}{l}{a}_{\mathrm{1}}=(\mathrm{1}/k_{\alpha }+\mathrm{1})m_{A}g{y}_A/r_X^{(P)}\\ a_{\mathrm{2}}=-(\mathrm{1}/k_{\alpha }+\mathrm{1})m_{A}g{z}_A/r_X^{(P)}\\ a_{\mathrm{3}}=(\mathrm{4}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}l+\mathrm{16}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}})/(\mathrm{24}EIl)\\ a_{\mathrm{4}}=-(\mathrm{3}{q}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{3}}+\mathrm{4}{F}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}}{l}^{\mathrm{2}})/(\mathrm{24}EIl)\end{array}\right\}（19）

则式（18）可变为

\mathrm{\Delta }{\alpha }^{\text{'}}=a_{\mathrm{1}}(\mathrm{c}\mathrm{o}\mathrm{s}\alpha -\mathrm{1})+a_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha +a_{\mathrm{3}}{y}^{\mathrm{3}}+a_{\mathrm{4}}y（20）

式（20）揭示了各轴运动对重力变形分量\mathrm{\Delta }{\alpha }^{\text{'}}的影响规律。此外，由于式（18）中弹性系数r_X^{(P)}和弹性模量E均较大，因此引入的4个运动学参数的数值均很小。

同理，对式（15）中的其他重力变形分量进行类似的处理，最终得到各轴运动对各重力变形分量的影响规律，其表达式如下：

\left.\begin{array}{l}\mathrm{\Delta }{\alpha }^{\text{'}}=a_{\mathrm{1}}(\mathrm{c}\mathrm{o}\mathrm{s}\alpha -\mathrm{1})+a_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha +a_{\mathrm{3}}{y}^{\mathrm{3}}+a_{\mathrm{4}}y\\ \begin{array}{l}\mathrm{\Delta }{\beta }^{\text{'}}=a_{\mathrm{5}}(\mathrm{c}\mathrm{o}\mathrm{s}\beta -\mathrm{1})+a_{\mathrm{6}}\mathrm{s}\mathrm{i}\mathrm{n}\beta +a_{\mathrm{7}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta +\\ a_{\mathrm{8}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta \end{array}\\ \mathrm{\Delta }{\gamma }^{\text{'}}=\mathrm{0}\\ \begin{array}{l}\mathrm{\Delta }{x}^{\text{'}}=\mathrm{\Delta }{\beta }^{\text{'}}(l_{PQ}^{(Z)}-z+l_{PB})-(l_{PQ}^{(Z)}-z)[a_{\mathrm{9}}(\mathrm{c}\mathrm{o}\mathrm{s}\beta -\mathrm{1})+\\ a_{\mathrm{10}}\mathrm{s}\mathrm{i}\mathrm{n}\beta +a_{\mathrm{11}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta +a_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta ]\end{array}\\ \begin{array}{l}\mathrm{\Delta }{y}^{\text{'}}=\mathrm{\Delta }{\alpha }^{\text{'}}(l_{PQ}^{(Z)}-z+l_{PB})-(l_{PQ}^{(Z)}-z)[a_{\mathrm{13}}(\mathrm{c}\mathrm{o}\mathrm{s}\alpha -\mathrm{1})+\\ a_{\mathrm{14}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha ]\end{array}\\ \mathrm{\Delta }{z}^{\text{'}}=a_{\mathrm{15}}{y}^{\mathrm{4}}+a_{\mathrm{16}}{y}^{\mathrm{2}}\end{array}\right\}（21）

其中，a₅~a₁₆均为引入的运动学参数。由此可见重力变形共引入了18个运动学参数，具体参数如表2所示。将这些运动学参数代入龙门铺丝机重力变形矩阵\mathrm{\Delta }{\mathit{T}}_{\mathrm{g}}中，再将该矩阵结合运动学模型（式（2）），最终得到结合几何偏差和重力变形的龙门铺丝机综合运动学模型，其表达式如下：

{}_{\mathrm{H}}{}^{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}}\mathit{T}(\mathit{q},{\mathit{p}}_{\mathrm{d}},{\mathit{p}}_{\mathrm{g}})={}_{\mathrm{A}\mathrm{c}\mathrm{t}}{}^{\mathrm{B}\mathrm{a}\mathrm{s}\mathrm{e}}\mathit{T}\mathrm{\Delta }{\mathit{T}}_{\mathrm{g}}{}_X{}^{\mathrm{A}\mathrm{c}\mathrm{t}}\mathit{T}{}_Y{}^X\mathit{T}{}_Z{}^Y\mathit{T}{}_B{}^Z\mathit{T}{}_A{}^B\mathit{T}{}_C{}^A\mathit{T}{}_{\mathrm{H}}{}^C\mathit{T}（22）

其中，{\mathit{p}}_{\mathrm{d}}为由几何偏差引入的22个运动学参数向量，{\mathit{p}}_{\mathrm{g}}为由重力变形引入的18个运动学参数向量，可见龙门铺丝机综合运动学模型共引入40个运动学参数。

2 运动学参数标定

在所建的3个龙门铺丝机运动学模型中，引入的运动学参数虽提高了模型正确性，但增加了参数确定的复杂度。几何偏差与重力变形引入的参数无法从设备数模中直接获取，也不能通过测量仪器直接测得，必须借助运动学参数标定技术精准确定，该标定包含运动学建模、数据测量、参数辨识、模型修正四个步骤。

2.1 测量数据获取

测量数据的准确性是参数标定的前提条件。激光跟踪仪作为一种高精度三维坐标测量仪器，具有测量精度高、测量范围大、操作方便灵活等特点，本文采用激光跟踪仪获取铺丝机的实际测量数据。

鉴于激光跟踪仪只能测量三维空间中单点的位置坐标，而无法获取龙门铺丝机末端的位置和姿态，故本文设计了一种可利用激光跟踪仪测量末端位姿的测量工装，如图7所示。该工装上放置了3个激光跟踪仪靶座(P_O,P_X,P_Y)，它们的连线构成一个直角三角形。测量前将铺丝头末端压辊拆下，换上本文所设计的测量工装，此时直角三角形的两条直角边分别沿着龙门铺丝机X轴和Y轴的正方向。

测量某一位姿时，只需依次用激光跟踪仪测量这三个靶座上靶球的坐标，再利用右手螺旋定则便可快速构建用于描述末端位姿的坐标系\{S_{\mathrm{H}}\}，从而实现龙门铺丝机末端位姿的测量。

为保证辨识结果的精度，需测量足够数量的位姿点用于参数辨识和验证工作。文献［16］选取了35个测量点对21个参数进行了辨识，文献［17］选取了100个测量点对36个参数进行了辨识，均有效保证了辨识精度。本文龙门铺丝机综合运动学模型中未知运动学参数为40个，考虑运动学参数的数量和龙门铺丝机工作空间的大小，本文共随机选取了均匀分布在工作空间内的250个位姿用于龙门铺丝机参数辨识和验证过程，其中，180个位姿用于参数辨识过程，70个位姿用于验证模型及参数的准确性。

2.2 运动学参数辨识

获取得到测量数据后便可开展运动学参数辨识。常见的参数辨识算法包括高斯牛顿（Gauss-Newton，G-N）算法、鲍威尔（Powell）算法、L-M（Levenberg-Marquardt）算法等，其中G-N法十分依赖初值的准确性，且对残差比较敏感；Powell算法会随着参数维度的增加，使所需迭代次数和函数求值次数急剧增加，不适合多参数优化场景；L-M法作为G-N法的改进版，兼具收敛速度快、收敛性强、残差处理能力好、鲁棒性强等优势，因此本文选用L-M法作为运动学参数辨识的优化算法。

优化算法中目标函数决定优化方向，参数会向其值最小的方向迭代。为使综合模型计算位姿贴合实测位姿，本文兼顾位置与姿态误差，将目标函数定义如下：

min J=k_{\mathrm{P}}{\displaystyle \sum _{i=\mathrm{1}}^n}{‖{\mathit{P}}_{\mathrm{m}}^{(i)}-{\mathit{P}}_{\mathrm{p}}^{(i)}‖}_{\mathrm{2}}+k_{\mathrm{R}}{\displaystyle \sum _{i=\mathrm{1}}^n}{‖{\mathit{R}}_{\mathrm{m}}^{(i)}-{\mathit{R}}_{\mathrm{p}}^{(i)}‖}_{\mathrm{2}}（23）

式中：{\mathit{P}}_{\mathrm{m}}^{(i)}、{\mathit{R}}_{\mathrm{m}}^{(i)}分别为第i个实际末端测量位姿的位置矩阵和姿态矩阵；{\mathit{P}}_{\mathrm{p}}^{(i)}、{\mathit{R}}_{\mathrm{p}}^{(i)}分别为第i个运动学模型预测位姿的位置矩阵和姿态矩阵；n为用于辨识过程的测量数据数量，本文取n=\mathrm{180}；k_{\mathrm{P}}、k_{\mathrm{R}}分别为位置误差和姿态误差的权重，两者的确定取决于位置误差和姿态误差的数量级，位置误差的数量级在毫米级，姿态误差的数量级一般在{\mathrm{10}}^{-\mathrm{2}}，为防止参数辨识过程对其中一部分误差不敏感，目标函数的建立需要通过两个权重系统对两部分误差进行归一化处理，即k_{\mathrm{P}}=\mathrm{1}，k_{\mathrm{R}}=\mathrm{100}。

优化算法的优化变量为运动学模型中的运动学参数，选择优化变量初值时，由几何偏差引入的运动学参数初值均设置为0，而由重力变形引入的运动学参数初值见表2。此外，优化算法的收敛准则采用点距准则（收敛容差为{\mathrm{10}}^{-\mathrm{6}}）和函数下降量准则(\mathrm{收敛}\mathrm{容差}\mathrm{为}{\mathrm{10}}^{-\mathrm{6}})。

完成优化算法目标函数和算法参数设置后即可进行运动学参数辨识过程，辨识得到的运动学参数随后代入到相应运动学模型中，最终实现完整的运动学模型的建立。

3 实验验证

在大型龙门铺丝机上开展相关实验，如图8所示。本实验所采用的激光跟踪仪型号为Leica AT901，其重复测量精度和绝对测量精度分别为±5 μm/m和±10 μm/m。为保证测量数据的可靠性，在实验过程中跟踪仪不受干扰。此外，实验前还需要将龙门铺丝机铺丝头的压辊更换为测量工装。

在测量时，首先采用激光跟踪仪建立龙门铺丝机基坐标系；然后，在控制系统中输入各轴坐标，并控制龙门铺丝机各轴运动到指定位置；最后，采用激光跟踪仪依次测量铺丝头末端测量工装三个靶座的坐标，从而完成一个位姿的测量。测量的最终结果便是得到一系列龙门铺丝机各轴运动量及对应的末端实测位姿。在此基础上，采用交叉验证的思想将这些位姿随机分为两组，其中一组180个位姿用于运动学参数辨识，另一组70个位姿用于辨识结果的验证。

为验证本文建立的综合运动学模型的准确性，现对建立的三个运动学模型进行分析和对比。将理论运动学模型定义为模型Ⅰ、结合几何偏差的运动学模型定义为模型Ⅱ、结合几何偏差和重力变形的综合运动学模型定义为模型Ⅲ，并根据测量数据和辨识算法对模型Ⅱ和模型Ⅲ开展运动学参数辨识。三个模型的参数辨识结果如表3和表4所示，参数的含义及初值见表1和表2。

从表3和表4中可以看出，模型Ⅱ和模型Ⅲ在由几何偏差引入的运动学参数部分相差不大，且与理论值接近，通过该相互验证说明模型Ⅱ和模型Ⅲ的辨识结果是准确可信的；而模型Ⅲ因存在由重力变形引入的运动学参数而不完全与理论参数相同，说明这些变化的参数产生了模型修正效果。

龙门铺丝机末端为柔性压辊，位置精度对铺丝精度的影响远大于姿态精度，故本文仅以位置误差为模型精度评判标准。图9展示了三个模型经参数辨识后的位置残余误差，可见模型Ⅱ、Ⅲ的位置残差均小于模型Ⅰ的位置残差，且模型Ⅲ的残差减小程度更优，Y、Z向的位置残差尤为明显，更贴合设备实际情况。模型Ⅰ的Z向位置残差约为2 mm，远大于X、Y向的位置残差，凸显重力变形对设备精度的显著影响；而与模型Ⅱ相比，模型Ⅲ的Z向位置精度控制能力更优，误差可控制在0.2 mm内，验证了模型中重力变形引入参数的有效性。

三个运动学模型在70个验证点处的预测误差如图10所示，可以看出，预测误差和辨识误差具有类似的误差分布形式，说明本文所建立的三个运动学模型均具有较好的一致性。

从图9和图10中还可以看出，模型Ⅰ在X方向上的位置误差主要集中在零线以下，而在Z方向上的位置误差主要集中在零线以上。这是因为龙门铺丝机Z轴滑枕绕横梁的旋转变形以及横梁、旋转轴框架等在Z方向上的重力变形使得龙门铺丝机末端更趋向于沿着X正向和Z负向，因此在模型Ⅰ计算值减去实际测量值时，X方向误差偏向负值，Z方向误差偏向正值。而模型Ⅲ通过引入重力变形部分的运动学参数，能很好地表征这部分的固有偏差，使得该模型在三个方向上的位置误差能够均匀分布在零线附近，从而显著提高了模型预测精度。

最后，对三个运动学模型的残余误差和预测误差进行数据统计分析，对比列出了三者误差求绝对值后的平均值和最大值，并计算出了模型Ⅱ和模型Ⅲ相对于模型Ⅰ的误差减小百分比，如表5所示，其中\mathrm{\Delta }x、\mathrm{\Delta }y、\mathrm{\Delta }z分别为三个方向上的位置误差。

由表5可知，与模型Ⅰ相比，模型Ⅱ、Ⅲ在三个方向上的位置误差均有减小：模型Ⅱ的位置误差均值减小了56.9%~75.8%，模型Ⅲ的位置误差均值减小了74.0%~88.6%，模型Ⅲ的减小范围更大、精度更高；且模型Ⅲ的Z向误差减小程度高于其他两个方向的误差减小程度，而模型Ⅱ无此特征，由此证明综合运动学模型对重力变形的预测有效。

4 结论

1）基于多体系统理论建立龙门铺丝机理论运动学模型，结合几何偏差构建对应模型；重点研究设备重力变形机理，分析旋转轴框架、Y轴滑台、龙门架横梁的结构变形对运动学精度的影响规律，融入模型后建立综合运动学模型。

2）针对模型未知参数，研究了基于激光跟踪仪的位姿测量、测点选择等方法；通过确定L-M算法的目标函数、优化变量初值与收敛准则，实现运动学参数辨识。

3）开展模型参数辨识与验证实验，结果表明融合几何偏差和重力变形的综合模型精度显著优于仅含几何偏差的模型，可有效表征设备X、Z向的固有偏差，位置误差平均值较理论模型减小74.0%以上，多数控制0.2 mm内，且工作空间精度一致性良好，验证了所提综合建模及参数辨识方法能有效提高龙门铺丝机静态精度。

参考文献

[1] 黄亿洲， 王志瑾， 刘格菲. 碳纤维增强复合材料在航空航天领域的应用［J］. 西安航空学院学报， 2021， 39（5）： 44-51.

[2] HUANG Yizhou， WANG Zhijin， LIU Gefei. Application of Carbon Fiber Reinforced Composite in Aerospace［J］. Journal of Xi'an Aeronautical University， 2021， 39（5）： 44-51.

[3] 郝元恺， 肖加余. 高性能复合材料学［M］. 北京： 化学工业出版社， 2004.

[4] HAO Yuankai， XIAO Jiayu. High Performance Composite Materials Science［M］. Beijing： Chemical Industry Press， 2004.

[5] 张晓虎， 孟宇， 张炜. 碳纤维增强复合材料技术发展现状及趋势［J］. 纤维复合材料， 2004， 21（1）： 50-53.

[6] ZHANG Xiaohu， MENG Yu， ZHANG Wei. The State of the Art and Trend of Carbon Fiber Reinforced Composites［J］. Fiber Composites， 2004， 21（1）： 50-53.

[7] 王显峰， 严飙， 薛柯， 等. 机器人高效自动铺丝技术研究进展［J］. 航空制造技术， 2019， 59（16）： 14-20.

[8] WANG Xianfeng， YAN Biao， XUE Ke， et al. Research Progress of Robot Automatic Placement Technology with High Efficiency［J］. Aeronautical Manufacturing Technology， 2019， 59（16）： 14-20.

[9] ZHAO Dan， BI Yunbo， KE Yinglin. An Efficient Error Compensation Method for Coordinated CNC Five-axis Machine Tools［J］. International Journal of Machine Tools and Manufacture， 2017， 123： 105-115.

[10] ZHANG Zhao， JIANG Feng， LUO Ming， et al. Geometric Error Measuring， Modeling， and Compensation for CNC Machine Tools： a Review［J］. Chinese Journal of Aeronautics， 2024， 37（2）： 163-198.

[11] 谢核. 机器人加工几何误差建模及工程应用［D］. 武汉： 华中科技大学， 2019.

[12] XIE He. Geometric Error Modeling and Its Engineering Application in Robotic Machining［D］. Wuhan： Huazhong University of Science and Technology， 2019.

[13] 董泽园， 李杰， 刘辛军， 等. 数控机床两种几何误差建模方法有效性试验研究［J］. 机械工程学报， 2019， 55（5）： 137-147.

[14] DONG Zeyuan， LI Jie， LIU Xinjun， et al. Experimental Study on the Effectiveness of Two Different Geometric Error Modeling Methods for Machine Tools［J］. Journal of Mechanical Engineering， 2019， 55（5）： 137-147.

[15] 郭世杰， 梅雪松， 姜歌东. 基于平面光栅的机床几何误差测量与辨识［J］. 计算机集成制造系统， 2020， 26（8）： 2037-2049.

[16] GUO Shijie， MEI Xuesong， JIANG Gedong. Geometric Error Measurement and Identification of Machine Tool Based on Gross Grid Encoder［J］. Computer Integrated Manufacturing Systems， 2020， 26（8）： 2037-2049.

[17] ZHU Shaowei， DING Guofu， QIN Shengfeng， et al. Integrated Geometric Error Modeling， Identification and Compensation of CNC Machine Tools［J］. International Journal of Machine Tools and Manufacture， 2012， 52（1）： 24-29.

[18] 杨建军， 司林林， 林守金， 等. 齿轮加工机床几何误差补偿研究综述［J］. 机床与液压， 2024， 52（19）： 174-187.

[19] YANG Jianjun， SI Linlin， LIN Shoujin， et al. Review on Geometric Error Compensation of Gear Processing Machine Tools［J］.Machine Tool & Hydraulics，2024，52（19）： 174-187.

[20] 范晋伟， 秦池， 潘日. 基于数控指令修正的数控内圆磨床几何误差补偿［J］. 机床与液压， 2022， 50（21）： 157-162.

[21] FAN Jinwei， QIN Chi， PAN Ri. Geometric Error Compensation of CNC Internal Grinder Machine Based on CNC Command Correction［J］. Machine Tool & Hydraulics， 2022， 50（21）： 157-162.

[22] 李忠群， 刘鸿， 刘强， 等. 四轴联动机床几何误差建模分析与评估［J］. 机床与液压， 2024， 52（23）： 175-180.

[23] LI Zhongqun， LIU Hong， LIU Qiang， et al. Analysis and Evaluation of Geometric Error Modeling for Four-axis Linkage Machine Tools［J］. Machine Tool & Hydraulics， 2024， 52（23）： 175-180.

[24] 王明海， 张威， 刘香辰. 基于位姿变化的钻铆机器人刚度优化［J］. 机床与液压， 2024， 52（3）： 49-54.

[25] WANG Minghai， ZHANG Wei， LIU Xiangchen. Drilling Riveting Robot Stiffness Optimization Based on Pose Change［J］. Machine Tool & Hydraulics， 2024， 52（3）： 49-54.

[26] 陆历历， 郭旭红， 吉文正. TH6920滑枕挠曲变形的耦合补偿［J］. 机械设计与制造， 2019（1）： 261-263.

[27] LU Lili， GUO Xuhong， JI Wenzheng. Coupling Compensation of the Ram Flexure Deformation for TH6920［J］. Machinery Design & Manufacture， 2019（1）： 261-263.

[28] ZHU Weidong， MEI Biao， KE Yinglin. Kinematic Modeling and Parameter Identification of a New Circumferential Drilling Machine for Aircraft Assembly［J］. The International Journal of Advanced Manufacturing Technology， 2014， 72（5）： 1143-1158.

[29] ZHAO Dan， BI Yunbo， KE Yinglin. Kinematic Modeling and Base Frame Calibration of a Dual-machine-based Drilling and Riveting System for Aircraft Panel Assembly［J］. The International Journal of Advanced Manufacturing Technology， 2018， 94（5）： 1873-1884.

基金资助

航空科学基金(2024Z043076001)

浙江省自然科学基金(LY24E050003)