基于Kriging模型插值的孔位修正策略

石循磊，张继文，刘顺涛，陈恳,*

1. 清华大学 机械工程系，北京 100084 2. 航空工业成都飞机工业(集团)有限责任公司，成都 610091

飞机普遍采用铆接或螺栓连接的形式实现零部件的装配，叠层制孔和连接任务占据了大部分装配任务[1]。现代飞机一般具有复杂的气动外形，且零件刚度较低，装配时容易产生较大的变形和位置偏移[2-3]，因此难以像普通加工机床一样，通过离线编程和在线对刀的形式一次性完成制孔作业任务。为了适应零部件的变形和位置移动，保证制孔孔位的边距准确性，在机器制孔之前通过人工方式或特殊方式制备基准孔(或称预连接孔)，再通过激光跟踪仪或者其他视觉测量系统检测这些基准孔的实际位置，将其与理论位置相比较，计算出基准孔的孔位偏差[4-5]。依据已制的基准孔位偏差，通过插值算法推算出待制孔的位置偏差，再对自动制孔程序进行校正，调整实际制孔位置，实现孔位补偿的功能。

在自动孔位寻找校正过程中，通过已制基准孔位推断待制孔孔位偏差的补偿算法显得尤为重要。然而，飞机零件的弱刚性以及复杂外形特性为孔位推算带来了多个难点：其一，装配过程中，工件的变形并非简单的刚性位置移动或旋转，通常还伴随有较大幅度的弹性扭转和弯曲[6]，因此无法基于一系列基准孔采用装配体的刚体变换形式获得待制孔孔位;其二，基准孔通常位于不同的零件装配件上，例如飞机的隔框纵梁长桁结构交替出现，基准孔也随之交替出现，若基准孔位与零件绑定，则为鉴别和记录带来繁杂的标记过程，因而极易出现错误，进而导致整个孔位校正失效。

目前大多数工程化解决方案采用线性或三次样条插值的方法校正两个或多个基准孔之间的孔位偏差。这是一种近似的解决方案，飞机的复杂曲面及变形特征很难服从这一假设，特别在一些曲率变化较大的位置，发生孔位校正错误的概率大大增加[7]。为了保证校正的局部特性，不得不增加基准孔的数量，并严重依赖工艺实验和人工经验选择基准孔位[8]。同时线性插值、三次样条插值需要选定插值节点，然而，对于三维空间位置而言，很难区分样值点的内外关系，也难以确定插值节点，而且在外插值情况下容易出现较大误差[9]。

在理论研究方面，Zhu等[10]通过基准孔构建了双线性误差平面，实现了对制孔区域的孔位进行线性插值补偿。周叔阳等[11]将模拟退火算法应用于预连接孔的排布问题，并抽象为数学模型，给出了预连接孔排布优化算法及流程。董辉跃等[12]从飞机曲面构造原理角度，提出了一种基于误差Coons曲面函数的补偿方法。王龙飞等[13]通过视觉辅助和手眼标定，计算出基准孔对应的机器人驱动坐标误差变换矩阵，进而补偿待制孔误差。王青等[14]提出一种基于孔边距偏差修正约束孔位和利用Shepard插值修正其余孔位的策略，实现了较好的孔位修正效果，但是缺少对孔位预测误差的分析以及基准孔的补制策略。本文提出一种采用Kriging模型插值算法的孔位偏差校正算法。该方法利用Kriging插值模型，依据基准孔位偏差推算待制孔偏差和预测标准差。算法依赖孔位估计的均方误差评估基准孔的数量和分布质量，进而可以在基准孔实际选择和布置时，依赖反复的迭代和计算，指导增添新的基准孔，从而提升孔位偏差估计的精度。

1 孔位补偿原理

根据基准孔位偏差推算待制孔位置偏差的原理如图1所示。

图1 由基准孔位偏差推算待制孔位置偏差的原理Fig.1 Principle of estimating position deviation of hole to be made from reference hole

已经依据理论位置及工件装配实际状况制备出N个基准孔。基准孔实际位置与CAD模型中的理论位置存在一定偏差。

(1)

设CAD模型中的一个待制孔理论位置坐标为(u*,v*,w*)T，则补偿目标是推算出输出实际制孔位置与待制孔理论位置的偏差值：

(2)

上述过程本质是分别建立了孔位偏差的3个空间坐标分量Δu、Δv、Δw与孔位理论位置[u,v,w]T之间的响应面模型。建立这种响应面关系的一种方法是通过材料力学、弹性力学等复杂的理论推导出这一关系，但难度过大，而本文方法通过多次实验的点以及一种特殊的函数关系模拟出这种关系，建立Kriging模型。

2 Kriging模型建立

Kriging模型属于响应面模型，其本质是一种数据拟合的过程，将函数视作一种随机过程，通过有限的数据拟合出这一随机过程的概率分布规律，从而得到函数，计算出在空间坐标变化时函数数值的变化规律，进而计算出未知点的数值。Kriging模型具有良好的非线性函数拟合特性，不仅提供了未知点的预测值，还提供了预测标准差[15]。

Kriging模型的基本表达式为

y(x)=fT(x)β+μ(x)

(3)

式中：f(x)为q维基函数向量；β为q维系数向量；μ(x)为随机变量，每给出一个采样点(或称自变量)x，在f(x)、β和μ(x)已知的条件下，即可计算出其响应函数值y(x)[16]。β和μ(x)未知时，当具有一系列采样点x(1),x(2),…,x(N)及对应的y(1),y(2),…,y(N)，即可估计出其参数，进而求得其他位置x*的函数值y(x*)。

如果将μ(x)当做满足正态分布、与位置无关的随机变量μ，则可以通过多次实验消除。然而，在Kriging模型内，μ(x)是一个依赖于x的随机变量，且两个样值点x(i)、x(j)距离越近，其相关性越强。针对全局范围内的采样点x，满足数学期望E[μ(x)]=0，同时，对于任意两个采样点x(i)、x(j)，满足协方差关系：

(4)

式中：σ2为采样点对应的随机变量的方差；R(θ,p,x(i),x(j))为关联函数，体现了x(i)、x(j)两个采样点间的相关系数：

(5)

式中：θh为x向量中h分量的相关强度；ph为x向量中h分量的相关形式[17]。与拟合回归不同，Kriging模型主要通过μ(x)描述y(x)的响应值，并通常将fT(x)设置为简单的1T常数。该模型假设μ(x)和y(x)仅与状态x有关，因而不存在多项式拟合中的拟合误差。

在fT(x)=1T和p预先指定的条件下，当给出多组采样点时，通过最大似然估计法，估算参数{β,σ2,θ}，从而得到模型估计结果：

(6)

式(6)显示Kriging模型估计的结果是一系列关联函数的线性组合所表达的连续可微函数式：

(7)

式中：ai为最大似然假设所估计的N个常数[19]。

依据R(θ,x*,x(i))的定义(设ph=2时)：

(8)

3 孔位偏差估计及改进策略

针对孔位偏差问题，距离待制孔距离越大的基准孔位偏差对当前待制孔偏差的影响越小，而距离待制孔距离越小的基准孔孔位偏差则对当前待制孔的影响较大，这种特性能够有效描述待制孔与基准孔之间的偏差关联关系。如图2所示。

根据Kriging模型的定义及参数估计方法，当x*=x(i)时，即取已知实验点时，由响应面函数计算可得：

(9)

式(9)表明由式(6)计算出的基准孔x(i)处的偏差值仍然是实测值y(i)，且其估计偏差为0。因此，Kriging模型所建立的是对已知实验点的一种特殊的插值，而相比于线性插值或样条插值，孔位偏差不再限于内插值，而是可以实现一定程度的外插值。同时对于一个待制孔的误差，其不仅取决于两端的基准孔，还取决于与其临近的若干个孔的偏差，这与一般装配的误差特性一致。

图2 待制孔与基准孔之间的偏差关联关系Fig.2 Deviation correlation between hole to be made and reference hole

(10)

(11)

因此，采用Kriging模型插值推算孔位偏差，能够依赖孔位估计的均方误差评估基准孔的数量和分布质量，进而可以在基准孔实际选择和布置时，依赖反复的迭代和计算，指导增添新的基准孔，从而提升孔位偏差估计的精度。此外，可以通过式(7)采用一种连续可微的函数形式模拟孔位偏差，从而绘制出孔位误差分布图，有利于装配质量以及基准孔制备的迭代改进。算法的总体结构如图3所示。

图3 算法总体结构Fig.3 Overall structure of algorithm

4 实验与分析

飞机孔位偏差来源可分为待制孔零件的变形和平移旋转两种，实验中分别验证两种偏差来源形式下本文算法的有效性。

针对零件变形带来的孔位偏差，采用有限元仿真实验来验证算法的有效性。利用有限元软件Abaqus建模，如图4所示，将蒙皮简化为边长约4.8 m、厚度3 mm的梯形曲面工件，连接于骨架上，并施加载荷。坐标系y轴定义为与曲面蒙皮3个相距较远顶点组成的平面垂直的方向，即y轴近似与蒙皮表面垂直，而x和z轴在该平面内。

图4 有限元仿真载荷示意图Fig.4 Load diagram of finite element simulation

有限元仿真结果如图5所示，图中标出了25个基准孔，在蒙皮上另取50个点作为待制孔的位置。将有限元仿真结果中基准孔的位置作为实际孔位，可得基准孔的位置偏差(位移)Δu、Δv、Δw作为Kriging模型的因变量，以基准孔理论位置为自变量，可建立3个Kriging模型。利用建立好的模型估计待制孔的应制孔位置，结果如图6所示。

图6中，eu、ev、ew是Kriging模型估计的孔位偏差值(位移)的误差，定义为

(12)

从图6中可以看出，变形导致的x和z方向位移量较小，孔位估计误差均小于0.1 mm。y方向位移达到10 mm量级，可认为是实际装配时壁板的极限变形，因此主要分析y方向实验结果。

针对零件整体平移及旋转带来的孔位偏差，偏移前后孔位坐标有Q=R0Q0+T0

(13)

式中：Q0为孔位理论坐标；Q为偏移后坐标；R0为旋转矩阵；T0为平移变换矩阵，因此孔位偏差量为Q-Q0。

在实际装配过程中，蒙皮壁板放置在飞机骨架上，受骨架约束，可认为壁板仅绕图8中y轴旋转。验证实验中取旋转角度为1.15°，可得旋转矩阵：

图5 有限元仿真位移结果云图Fig.5 Displacement nephogram of finite element simulation

(14)

同样采用图5中的25个基准孔建立Kriging模型，50个待制孔的x方向和z方向最大位移均达10 mm量级，而孔位估计误差如图9所示，均小于0.05 mm。因此对于零件整体旋转平移导致的孔位偏差，该算法能很好地估计偏移后的位置。

利用上述方法在某型飞机数字化装配现场实验的结果见表1，实物实验中基准孔的偏移量均小于1 mm,该方法估计出的待制孔的偏差量与测量值相差很小。

考虑该算法的计算效率，实验使用的计算机CPU为i7-8700k，利用R语言编程，由25个基准孔建立Kriging模型并估计50个待制孔位所需运算时间为0.214 5 s，50个基准孔估计100个待制孔耗时0.306 2 s，90个基准孔估计150个待制孔耗时0.480 7 s，而实际生产中孔位估计并无实时性要求，所以计算效率完全能达到要求。与线性插值等传统方法相比较，该算法除了能给出待制孔的预测标准误差外，另一个显著优势体现在曲面区域外插值的情况。例如在上述有限元仿真实验中，一个待制孔的坐标为(-1 592.7, 75.251, -1 008.3)，该待制孔附近的两个基准孔坐标为(-1 692.1, 66.179, -1 334.5)和(-1 690.7, 66.312, -814.05)，以这两个点为插值节点，线性插值得待制孔的y方向偏差值为154.78 mm，而仿真结果为2.990 3 mm，本文算法估计结果为2.847 1 mm，此时线性插值出现很大误差，而该算法仍有较好效果。

图6 变形导致孔位偏差的实验结果Fig.6 Experimental results of hole displacement caused by deformation

图7 增添基准孔后实验结果Fig.7 Experimental results after adding reference holes

图8 实验件坐标系定义Fig.8 Definition of experimental coordinate system

多组实验表明本文方法能较准确、高效地通过基准孔信息估计待制孔的位移，并且给出有实际参考意义的待制孔位置的预测标准误差，用以指导排布或增添基准孔，最后通过在合适位置增添基准孔能有效提高待制孔位移估计的准确度。

图9 平移旋转导致孔位偏差的实验结果Fig.9 Experimental results of hole displacement caused by translation and rotation

表1 实物实验待制孔位估计结果

5 结 论

1) 提出了一种基于基准孔的孔位校正方法，该方法利用Kriging模型，依据基准孔位偏差推算出待制孔偏差和预测标准误差。

2)相比于线性或三次样条插值，该方法适用范围更大，且能给出孔位估计的标准误差，以评估基准孔的数量和分布质量，进而指导增添新的基准孔。

3)利用有限元仿真实验，验证了在零件变形及平移旋转导致孔位偏差情况下本文方法的有效性，在合理增添基准孔后，孔位估计的误差小于0.3 mm，满足航空装配要求。